教学项目4 数值分析法模型

教学项目四数值分析法模型4.1插值法建模拉格朗日插值分段线性插值三次样条插值一、插值的定义二、插值的方法三、用Matlab解插值问题已知n+1个节点,,1,0(),(njyxjj其中jx互不相同，不妨设),10bxxxan求任一插值点)(*jxx处的插值.*y0x1xnx0y1y节点可视为由)(xgy产生，g表达式复杂，或无封闭形式，或未知。*x*y4.1.1插值法的定义构造一个(相对简单的)函数),(xfy通过全部节点,即),1,0()(njyxfjj再用)(xf计算插值，即).(**xfy0x1xnx0y1y*x*y4.1.2插值的方法已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn。求一n次多项式函数Pn(x)，使其满足：拉格朗日(Lagrange)插值niyxPiin,,1,0,)(称为拉格朗日插值基函数。)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10iniiinyxLxP0)()(解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：特殊情况两点一次(线性)插值多项式:101001011yxxxxyxxxxxL三点二次(抛物)插值多项式:2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL.,满足插值条件直接验证可知xLn采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n取10时，绘出插值结果图形.11,911)(2xxxg例分段线性插值其它,0,,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.xjxj-1xj+1x0xnxoy比分段线性插值更光滑。xyxi-1xiab在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值},1],,[),({)(1nixxxxsxSiii],[)()3),1,0()()2),1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs)1,,1()()(),()(),()(111nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40nxSxS)(,,,)4)3)2xSdcbaiiiig(x)为被插值函数。)()(limxgxSnn用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’：最邻近插值‘linear’：线性插值；‘spline’：三次样条插值；‘cubic’：立方插值。缺省时：分段线性插值。注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。例：在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计每隔1/10小时的温度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接输出数据将是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作图xlabel('Hour'),ylabel('DegreesCelsius’)xyX035791112131415Y01.21.72.02.12.01.81.21.01.6实验习题已知飞机下轮廓线上数据如下，求x每改变0.1时的y值。机翼下轮廓线4.2拟合法简介2.拟合的基本原理1.拟合问题引例拟合问题引例1温度t(0C)20.532.751.073.095.7电阻R()7658268739421032已知热敏电阻数据：求600C时的电阻R。2040608010070080090010001100设R=at+ba,b为待定系数拟合问题引例2t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg)求血药浓度随时间的变化规律c(t).作半对数坐标系(semilogy)下的图形02468100101102为待定系数kcectckt,)(0曲线拟合问题的提法已知一组（二维）数据，即平面上n个点（xi,yi)i=1,…n,寻求一个函数（曲线）y=f(x),使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。+++++++++xyy=f(x)(xi,yi)ii为点（xi,yi)与曲线y=f(x)的距离拟合与插值的关系函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。实例：下面数据是某次实验所得，希望得到X和f之间的关系？x1247912131517f1.53.96.611.715.618.819.620.621.1问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：•若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。•若要求所求曲线（面）通过所给所有数据点，就是插值问题；最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416180510152025已已已已已spline已已已已已已已0246810121416180510152025已已已已已linest已已已已已已已0246810121416180510152025已已已已已nearest已已已已已已已曲线拟合问题最常用的解法——线性最小二乘法的基本思路第一步:先选定一组函数r1(x),r2(x),…rm(x),mn,令f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中a1,a2,…am为待定系数。第二步:确定a1,a2,…am的准则（最小二乘准则）：使n个点（xi,yi)与曲线y=f(x)的距离i的平方和最小。记)2(])([])([),,(211211221iiknimkkininiiimyxrayxfaaaJ问题归结为，求a1,a2,…am使J(a1,a2,…am)最小。线性最小二乘法的求解：预备知识超定方程组：方程个数大于未知量个数的方程组)(221111212111mnyarararyarararnmnmnnmm即Ra=ynmnmnnmyyyaaarrrrrrR112111211,,其中超定方程一般是不存在解的矛盾方程组。如果有向量a使得达到最小，则称a为上述超定方程的最小二乘解。212211)(imniimiiyararar线性最小二乘法的求解定理：当RTR可逆时，超定方程组（3）存在最小二乘解，且即为方程组RTRa=RTy的解：a=(RTR)-1RTy所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。nmnmnmyyyaaaxrxrxrxrR111111,,)()()()(其中Ra=y（3）线性最小二乘拟合f(x)=a1r1(x)+…+amrm(x)中函数{r1(x),…rm(x)}的选取1.通过机理分析建立数学模型来确定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.将数据(xi,yi)i=1,…n作图，通过直观判断确定f(x)：用MATLAB解拟合问题1、线性最小二乘拟合2、非线性最小二乘拟合用MATLAB作线性最小二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,可利用已有程序:a=polyfit(x,y,m)2.对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。，用yRa\3.多项式在x处的值y可用以下命令计算：y=polyval（a，x）输入同长度的数组X，Y拟合多项式次数即要求出二次多项式:3221)(axaxaxf中的),,(321aaaA使得:最小])([1112iiiyxf例对下面一组数据作二次多项式拟合xi0.10.20.40.50.60.70.80.91yi1.9783.286.167.347.669.589.489.3011.21）输入以下命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；A=R\y'1111211121xxxxR此时解法1．用解超定方程的方法2）计算结果：Ａ=-9.810820.1293-0.03170317.01293.208108.9)(2xxxf1）输入以下命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出数据点和拟合曲线的图形2）计算结果：Ａ=-9.810820.1293-0.0317解法2．用多项式拟合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，…，xdatan），ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非线性最小二乘拟合Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.最小)),((21niiiydataxdataxFlsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量),使得输入格式为:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’