数学建模插值方法.pptx

《数学建模插值方法.pptx》由会员分享，可在线阅读，更多相关《数学建模插值方法.pptx（75页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、前言 函数是多种多样的，在科研与工程实际中有的函数表达式过于复杂而不便于计算，但又需要计算多点的函数值；有的函数甚至给不出数学式子，只能通过实验和测量得到一些离散数据（如某些点的函数值和导数值）。面对这种情况，很自然的一个想法就是构造某个简单的函数作为要考察的函数的近似。如果要求近似函数满足给定的离散数据，则称之为插值函数。实用上，我们常取结构相对比较简单的代数多项式作为插值函数,这就是所谓的代数插值。第1页/共75页 设 为给定的节点，为相应的函数值，求一个次数不超过 的多项式 ，使其满足 ，.这类问题称为插值问题。称为被插值函数，称为插值函数，称为插值节点一、问题提出插值部分第2页/共75

2、页 定理1 设 为给定的彼此互异的 个插值节点，则存在唯一的次数不超过 的多项式 ，满足条件 ,.二、存在性与唯一性第3页/共75页证明证明:设设 ,其中其中 为待定系数为待定系数.利用插值条件利用插值条件 ,我们得到一个线性代数我们得到一个线性代数方程方程组组 ,其中其中 观察发现矩阵观察发现矩阵A A是范德蒙矩阵是范德蒙矩阵,那么那么,由几代知识知道矩阵由几代知识知道矩阵A A 的行列的行列式式 为为 ,由定理中条件由定理中条件,插值结点为彼此互异的插值结点为彼此互异的,那那么行么行列式不为零列式不为零.故由故由CramerCramer法则知线性代数方程组法则知线性代数方程组 存在唯一解存

3、在唯一解.第4页/共75页三、三、Lagrange插值法插值法 （1）Lagrange插值多项式可以表示为 第5页/共75页 引入记号引入记号 ,易证易证 ,从而从而LagrangeLagrange插值多项式可表示为插值多项式可表示为 第6页/共75页（2）插值误差估计 定理2 设 在 上连续，在 内存在,节点 ，是拉格朗日插值多项式，则对任意 ,插值余项 其中 且依赖于 .第7页/共75页例2.求过点(2,0)(4,3)(6,5)(8,4)(10,1)的拉格朗日型插值多项式。解：用4次插值多项式对5个点插值 第8页/共75页 第9页/共75页于是有第10页/共75页function yi=l

4、agrcz(x,y,xi)n=length(x);m=length(xi);for s=1:m yi(s)=0;for i=1:n w(i)=1;dw(i)=1;for j=1:n if(j=i)w(i)=(xi(s)-x(j)*w(i);dw(i)=(x(i)-x(j)*dw(i);end end yi(s)=y(i)*w(i)/dw(i)+yi(s);endend第11页/共75页第12页/共75页缺点:当增加或减少插值节点时,基函数需要重新 构造,不便于实际的计算使用第13页/共75页 定义称 为 在 两点处的一阶差商.（1）差商定义四、Newton插值法二阶差商n 阶差商第14页/共7

5、5页(2)Newton插值公式插值公式 由差商定义由差商定义把以上各式由后向前代入把以上各式由后向前代入,可得可得第15页/共75页差商表 一阶一阶差商差商二阶差商二阶差商三阶差商三阶差商四阶差商四阶差商第16页/共75页例2:已知求满足以上插值条件的牛顿型插值多项式。解:1 2 3 4 0 -5 -6 3一阶差商二阶差商三阶差商 1 2 3 4 0 -5 -6 3 -5 -1 9 2 5 1第17页/共75页由上述差商表对角线上取得的值则牛顿三次插值多项式为 第18页/共75页function yi=newtcz(x,y,xi)n=length(x);m=length(xi);nt=zero

6、s(n,n);nt(:,1)=y;for i=2:n for j=i:n nt(j,i)=(nt(j,i-1)-nt(j-1,i-1)/(x(j)-x(j-(i-1);endEndfor i=1:n nt(i,i)Endfor i=1:m yi(i)=nt(1,1);for j=2:n t=1;for s=1:j-1 t=t*(xi(i)-x(s);end yi(i)=yi(i)+t*nt(j,j);endend第19页/共75页五、Hermite插值多项式给定的是节点上的函数值和导数值问题：已知求3次多项式 ，使得第20页/共75页第21页/共75页*多项式插值的问题 前面介绍了构造插值公式

7、的方法，并分析了前面介绍了构造插值公式的方法，并分析了它它们的余项。在实际应用插值函数作近似计算时，们的余项。在实际应用插值函数作近似计算时，总总希望插值公式余项希望插值公式余项 的绝对值小一些，即使的绝对值小一些，即使得得 逼近的精度好。从表达式看，似乎提高插值多项逼近的精度好。从表达式看，似乎提高插值多项式式的次数便可达到目的，但实际上并非如此。的次数便可达到目的，但实际上并非如此。第22页/共75页例如 给定函数 取其等距节点 ，构造的Lagrange插值多项式为 当 时，只能在 内收敛，而在这个区间以外是发散的。这种畸形现象 通常叫做Runge现象。如下图所示。第23页/共75页第24

8、页/共75页六、分段插值 所谓分段插值，就是将被插值函数逐段多项式化。在每个 子段上构造插值多项式，然后把它们装配在一，作为整个区间 上的插值函数，即称为分段多项式。如果函数 在每个子段上都是 次式，则称为 次式。一般（低次：k=1,2,3）第25页/共75页（1）分段线性插值的构造（k=1)易知 在每个子区间 上是一次插值多项式分段线性插值的余项其中第26页/共75页（2）分段抛物线插值（K=2）（3）分段三次 Hermite 插值(K=3)第27页/共75页（4）三次样条插值 在分段插值中，分段线性插值在节点上仅连续而不可导，分段三次埃尔米特插值有连续的一阶导数，如此光滑程度常不能满足物理

9、问题的需要，而引入的样条函数则可以同时解决这两个问题,使插值函数既是低阶分段函数,又是光滑的函数。第28页/共75页三次三次样条函数定条函数定义 给定区定区间 的一个划分的一个划分 ，如果函数如果函数 满足：足：（1）在每一小区间上是三次多项式；（2）在每个内节点上具有二阶连续导数；（3）则称 是 在该区间上关于该划分的一个三次样条函数。第29页/共75页其中四个待定系数为 ,子区间共有n n个所以要确定S(x)S(x)需要4n4n个待定系数。另一方面,要求分段三次多项式S(x)S(x)及其导数 和 在整个插值区间 a,ba,b 上连续,则要求它们在各个子区间的连接点 上连续，即满足条件 由样

10、条函数的定义可知,三次样条插值函数S(S(x x)是一个分段三次多项式,要求出S(S(x x),),在每个小区间 x xi i,x xi+1i+1 上要确定4 4个待定参数,若用S Si i(x x)表示它在第i i个子区间 x xi i,x xi+1i+1 上的表达式，则第30页/共75页（1 1）插值条件 （2 2）连接条件 式共给出了4n-24n-2个条件,而待定系数有4n4n个,因此还需要2 2个条件才能确定S(x),S(x),通常在区间端点上 各加一个条件,称为边界条件,常用边界条件有三种类型。第31页/共75页第一种类型：给定两端点 的一阶导数值：第二种类型：给定两端点f(x)f(

11、x)的二阶导数值：作为特例,称为自然边界条件。满足自然边界条件的三次样条插值函数称为自然样条插值函数。第三种类型：当 是以为 周期的函数时，则要求S(x)S(x)也是周期函数,这时边界条件应满足当 时，第32页/共75页这样，由上给定的任一种边界条件加上插值条件和连接条件，就能得出4n4n个方程，可以惟一确定4n4n个系数。从而得到三次样条插值函数S(x)S(x)在各个子区间 x xi i,x,xi+1i+1 上的表达式S(S(x xi i）(i=1,2,)(i=1,2,)。但是，这种做法当n n较大时，计算工作很大，不便于实际应用。因此我们希望找到一种简单的构造方法。第33页/共75页三次样

12、条插值函数的求法设S(x)S(x)在节点x xi i处的二阶导数为因为在子区间 x xi-1i-1,x,xi i 上 是三次多项式,所以 在此小区间上是x x的线性函数,且因为用线性插值,可知其表达式为记 ，则有第34页/共75页其中,A,Ai i,B,Bi i为积分常数,可利用插值条件 确定,即要求A Ai i,B,Bi i满足并记 ，则得连续两次积分得第35页/共75页由上讨论可知,只要确定 这n+1n+1个值,就可定出三样条插值函数S(x)S(x)。为了求出 ,利用一阶导数在子区间连接点上连续的条件 ，求导一次,得在区间 x xi-1i-1,x,xi i 上的表达式为 第36页/共75页

13、也就是在右端点x xi i上有 在左端点x xi-1i-1上有 将上式中的i-1i-1改为i,i,即得在子区间 x xi i,x,xi+1i+1 上的表达式 ,并由此得 利用 在内接点的连续性,即就可得到关于参数 的一个方程第37页/共75页上式两边同乘以 ,即得方程 若记 第38页/共75页则所得方程可简写成 即 这是一个含有n+1n+1个未知数、n-1n-1个方程的线性方程组.要完全确定 的值还需要补充两个条件,这两个条件通常根据实际问题的需要，根据插值区间 a,ba,b 的两个端点处的边界条件来补充。边界条件的种类很多，常见的有以下3 3种：第39页/共75页第一种边界条件：即已知插值区

14、间两端的一阶导数值：则可得到包含M Mi i的两个线性方程,S(x),S(x)在子区间 上的导数为由条件 得 即 同理,由条件 得 第40页/共75页即得确定 的线性方程组 其中第41页/共75页第二种边界条件:即已知插值区间两端的二阶导数值:,由于在区间端点处二阶导数 ，所以方程中实际上只包含有n-1n-1个未知数 ，从而得方程组 第42页/共75页第三种边界条件:由 与 ，可得 和 其中第43页/共75页得关于 的线性方程组。利用线性代数知识,可以证明方程组的系数矩阵都是非奇异的，因此有惟一解。第44页/共75页 用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时，其函数值在整体

15、上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数，不会发生龙格现象。因此三次样条在计算机辅助设计中有广泛的应用。第45页/共75页用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ：线性插值；线性插值；spline ：三次样条插值；三次样条插值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调

16、的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。第46页/共75页 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度依次为：度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的温小时的温度值。度值。x=1:12;y=5 8 9 15 25 29 31 30 22 25 27 24;xi=1:0.1:12;yi=interp1(x,y,xi,spline);plot(x,y,+,xi,yi,

17、r)第47页/共75页第48页/共75页 三次样条插值的Matlab实现 如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第2个三次多项式也做同样地处理。Matlab中三次样条插值也有现成的函数：y=interp1(x0,y0,x,spline)；y=spline(x0,y0,x)；pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)。其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。第49页/共75页对于三次样条插值，我们提倡使用函数csa

18、pe，csape的返回值是pp形式，要求插值点的函数值，必须调用函数ppval。pp=csape(x0,y0)：使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件，其值可为：complete 边界为一阶导数，一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。not-a-knot 非扭结条件periodic 周期条件second 边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为0,0。variational 设置边

19、界的二阶导数值为0,0。第50页/共75页对于一些特殊的边界条件，可以通过conds的一个12矩阵来表示，conds元素的取值为0，1，2。conds(i)=j的含义是给定端点i的j 阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件conds=2,1表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。第51页/共75页最小二乘法拟合 已知一批离散数据 (x xi i,y yi i),i i=0 0,1 1,.,n n，且 x x0 0 x x1 1 x xn n,寻找一个函数 f(x)，使达到最小.这个过程称为最小二乘拟合,f(x)称为拟合函数.

20、拟合部分第52页/共75页一、线性拟合 若设拟合函数f(x)=b+ax，则有令第53页/共75页即这是一个关于a a,b b的2 2元线性方程组.求解即可得到f f(x x)的表达式.第54页/共75页二、多项式拟合 有时所给数据点的分布并不一定近似地呈一条直线,这时仍用直线拟合显然是不合适的,可用多项式拟合。对于给定的一组数据 寻求次数不超过m m(mN)(mN)的多项式，来拟合所给定的数据，与线性拟合类似，使偏差的平方和为最小第55页/共75页由于 可以看作是关于 (j=0,1,2,j=0,1,2,m)m)的多元函数,故上述拟合多项式的构造问题可归结为多元函数的极值问题。令得 即有 第56

21、页/共75页这是关于系数 的线性方程组，通常称为正规方程组。可以证明，正规方程组有惟一解。第57页/共75页三、可化为线性拟合的非线性拟合 有些非线性拟合曲线可以通过适当的变量替换转化为线性曲线，从而用线性拟合进行处理，对于一个实际的曲线拟合问题，一般先按观测值在直角坐标平面上描出散点图，看一看散点的分布同哪类曲线图形接近，然后选用相接近的曲线拟合方程。再通过适当的变量替换转化为线性拟合问题，按线性拟合解出后再还原为原变量所表示的曲线拟合方程。下表列举了几类经适当变换后化为线性拟合求解的曲线拟合方程及变换关系 第58页/共75页 曲线拟合方程 变换关系 变换后线性拟合方程第59页/共75页多项

22、式曲线拟合函数：polyfit()polyfit()调用格式：p=polyfit(x,y,n)p=polyfit(x,y,n)p,s=polyfit(x,y,n)p,s=polyfit(x,y,n)说明：x,yx,y为数据点，n n为多项式阶数，返回p p为幂次从高到低的多项式系数向量p p。矩阵s s用于生成预测值的误差估计。例2 2：由离散数据x0.1.2.3.4.5.6.7.8.91y.3.511.4 1.6 1.9.6.4.81.5 2拟合出多项式。第60页/共75页x=0:.1:1;x=0:.1:1;y=.3.5 1 1.4 1.6 1.9.6.4.8 1.5 2y=.3.5 1 1

23、.4 1.6 1.9.6.4.8 1.5 2n=3;n=3;p=polyfit(x,y,n)p=polyfit(x,y,n)xi=linspace(0,1,100);xi=linspace(0,1,100);z=polyval(p,xi);z=polyval(p,xi);plot(x,y,o,xi,z,k:,x,y,b)plot(x,y,o,xi,z,k:,x,y,b)第61页/共75页第62页/共75页二维插值 xyO O第一种（网格节点）：第一种（网格节点）：第63页/共75页 注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值x y(x1,y1)(x1,y2)

24、(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。第64页/共75页 将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第65页/共75页 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性

25、插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O第66页/共75页(3)分片双三次样条插值第67页/共75页第二种（散乱节点）：第二种（散乱节点）：yx0 0第68页/共75页 要求x0,y0 x0,y0单调；x x，y y可取为矩阵，或x x取行向量，y y取为列向量，x,yx,y的值分别不能超出x0,y0 x0,y0的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearest nearest 最邻近插值linear linear 双线性插值cubic

26、 cubic 双三次插值缺省时,双线性插值第69页/共75页例：测得平板表面3*53*5网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2以平滑数据,

27、在x、y方向上每隔个单位的地方进行插值.第70页/共75页第71页/共75页再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.第72页/共75页第73页/共75页 插值函数griddata格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求cxcx取行向量，cycy取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值linear linear 双线性插值cubic cubic 双三次插值v4-Matlab提供的插值方法缺省时,双线性插值第74页/共75页感谢您的观看！第75页/共75页