第4讲 插值(简化)

后勤工程学院数学教研室MATLAB与数学实验数学教研室刘淑环第四讲插值实例：海底探测问题某公司用声纳对海底进行测试，在5×5海里的坐标点上测得海底深度的值，希望通过这些有限的数据了解更多处的海底情况。并绘出较细致的海底曲面图。在应用领域中，由有限个已知数据点，构造一个解析表达式，由此计算数据点之间的函数值，称之为插值。3实验目的实验内容2、掌握用MATLAB求解插值问题。1、了解插值的基本内容。[1]一维插值[2]二维插值[3]实验作业4拉格朗日插值分段线性插值三次样条插值一维插值一、插值的定义二、插值的方法三、用Matlab解插值问题5一维插值的定义已知n+1个节点,,1,0(),(njyxjj其中jx互不相同，不妨设),10bxxxan求任一插值点)(*jxx处的插值.*y0x1xnx0y1y节点可视为由)(xgy产生，g表达式复杂,或无封闭形式,或未知.。*x*y6构造一个(相对简单的)函数),(xfy通过全部节点,即),1,0()(njyxfjj再用)(xf计算插值，即).(**xfy0x1xnx0y1y*x*y7称为拉格朗日插值基函数。n0iiiny)x(L)x(P已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn。求一n次多项式函数Pn(x)，使其满足：Pn(xi)=yi,i=0,1,…,n.解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i拉格朗日(Lagrange)插值8拉格朗日(Lagrange)插值特别地:两点一次(线性)插值多项式:101001011yxxxxyxxxxxL三点二次(抛物)插值多项式:2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL拉格朗日插值的Matlab实现Matlab中没有现成的拉格朗日插值函数，必须编写一个M文件实现拉格朗日插值。设n个节点数据以数组x0,y0输入（注意Matlab的数组下标从1开始），m个插值点以数组x输入，输出数组y为m个插值。%编写一个名为lagrange.m的M文件：functiony=lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end1155,11)(2xxxg采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例12x0=[-5:2:5];y0=1./(1+x0.^2);x=[-5:0.1:5];y=lagrange(x0,y0,x);y1=1./(1+x.^2);plot(x,y,'-r')holdonplot(x,y1,'-b')-5-4-3-2-1012345-0.200.20.40.60.811.2插值的目的是函数逼近的一种手段。函数逼近的目的是得到一个数学问题的精确解或足够精确地解。那么，是否插值多项式的次数越高（即插值节点越多），越能够达到这个目的？下面看n=10的情况。55,11)(2xxxgx0=[-5:2:5];y0=1./(1+x0.^2);x=[-5:0.1:5];y=lagrange(x0,y0,x);y1=1./(1+x.^2);plot(x,y,'-r')holdonplot(x,y1,'-b')x0=[-5:1:5];y0=1./(1+x0.^2);x=[-5:0.1:5];y=lagrange(x0,y0,x);y1=1./(1+x.^2);plot(x,y,'-g')holdonplot(x,y1,'-b')n=5n=10-5-4-3-2-1012345-0.500.511.52n=5和n=10时的插值效果对比图n=10data2n=5data415由图看出，在0附近插值效果较好，即余项较小。另一种现象是插值多项式随节点增多而振动增大;虽然P10(x)在11个插值节点上与所逼近函数f(x)有相同值，但整体逼近效果较差，越靠近端点，逼近的效果就越差。这种现象称为龙格（Runge）现象。这也是高次插值多项式的病态现象。16x0=[-1:0.5:1];y0=1./(1+25*x0.^2);x1=[-1:0.2:1];x2=[-1:0.1:1];y1=lagrange(x0,y0,x1);y2=lagrange(x0,y0,x2);y3=1./(1+25*x2.^2)plot(x0,y0,'-r')holdonplot(x1,y1,'-b')holdonplot(x2,y2,'-g')holdonplot(x2,y3,'-y')-1-0.8-0.6-0.4-0.200.20.40.60.81-0.4-0.200.20.40.60.81步长0.5步长0.2步长0.1原函数从图中看出，在靠近-1或1时，余项会随n的增大而增大18上述现象和定理告诉我们，并不是插值多项式的次数越高（即插值节点越多）精度越高。从数值计算上可解释为高次插值多项式的计算会带来舍入误差的增大。从而引起计算失真。因此，实际应用做插值时一般只用一次、二次最多用三次插值多项式。那么如何提高插值精度呢？采用分段插值是一种办法19练习实验2021分段线性插值xjxj-1xj+1x0xnxoy可以看出，如果以每个小区间的端点为插值点的话，那么相邻区间的两个插值函数在节点处将保持连续。即P(x）是以连续函数，但在节点处一阶二阶导数并不一定连续。22其它,0,,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0xnxoy23求下列数据的多项式插值函数x=[012345678910];y=[12.32.124.64.74.38.19.29.810.3];plot(x,y,'*')012345678910123456789101124求下列数据的多项式插值函数x=[012345678910];y=[12.32.124.64.74.38.19.29.810.3];plot(x,y,'*')0123456789101234567891011plot(x,y,'*',x,y)思考：如何作分段的多项式插值函数？25x0=[-5:2.5:5];%n=4y0=5./(1+x0.^2);x=[-5:0.1:5];y1=5./(1+x.^2);plot(x0,y0,'*',x0,y0,x,y1,'r')-5-4-3-2-101234500.511.522.533.544.55散点图分段线性插值f(x）26x0=[-5:1:5];%n=10y0=5./(1+x0.^2);x=[-5:0.1:5];y1=5./(1+x.^2);plot(x0,y0,'*',x0,y0,x,y1,'r')-5-4-3-2-101234500.511.522.533.544.55散点图分段线性插值f(x）x1=[0:3];y1=[12.32.12];x0=[0:0.1:3];y0=lagrange(x1,y1,x0);x2=[3:6];y2=[24.64.74.3];x20=[3:0.1:6];y20=lagrange(x2,y2,x20);x3=[6:8];y3=[4.38.19.2];x30=[6:0.1:8];y30=lagrange(x3,y3,x30);x4=[8:10];y4=[9.29.810.3];x40=[8:0.1:10];y40=lagrange(x4,y4,x40);plot(x0,y0,'--r',x1,y1,'*',x20,y20,'--r',x2,y2,'*')holdonplot(x30,y30,‘--r’,x3,y3,‘*’,x40,y40,‘--r’,x4,y4,‘*‘)29比分段线性插值更光滑。xyxi-1xiab在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值30三次样条插值},1],,[),({)(1nixxxxsxSiii],[)()3),1,0()()2),1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs)1,,1()()(),()(),()(111nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40nxSxS)(,,,)4)3)2xSdcbaiiiig(x)为被插值函数。lim()()nSxgx31例66,11)(2xxxg用三次样条插值选取11个基点计算插值3266,11)(2xxxgx=-6:6;y=1./(1+x.^2);%插值节点共11个h=-6:0.1:6;%设置新的被插值点y1=interp1(x,y,h,'spline');%用三次样条方法求得相应插值plot(x,y,'+',x,y)%作原函数图形xlabel('x'),ylabel('y')%x轴y轴标签holdon%继续绘图plot(h,y1,'*',h,y1)%作插值点图形-6-4-2024600.10.20.30.40.50.60.70.80.91-6-4-2024600.10.20.30.40.50.60.70.80.91xyh=-6:0.1:6;%设置新的被插值点-6-4-2024600.10.20.30.40.50.60.70.80.91xyh=-6:0.5:6;%设置新的被插值点36用MATLAB作插值计算一维插值函数：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’：最邻近插值‘linear’：线性插值；‘spline’：三次样条插值；‘cubic’：立方插值。缺省时：分段线性插值。注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。3738程序：x=0:.1:1;y=[.3.511.41.61.6.4.81.52];yi0=interp1(x,y,0.025,'linear')xi=0:.02:1;yi=interp1(x,y,xi,'linear');zi=interp1(x,y,xi,'spline');wi=interp1(x,y,xi,'cubic');plot(x,y,'o',xi,yi,'r+',xi,zi,'g*',xi,wi,'k.-')legend('原始点','线性点','三次样条','三次多项式')结果：yi0=0.3500例1：已知数据如表所示，求当xi=0.25时的yi的值。x0.1.2.3.4.5.6.7.8.91y.3.511.41.61.9.6.4.81.523900.10.20.30.40.50.60.70.80.910.20.40.60.811.21.41.61.82原始点线性点三次样条三次多项式40要得到给定的几个点的对应函数值，可用：xi