第四章_多项式与插值--席.ppt

,单击此处编辑母版标题样式,*,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,第四章 多项式与插值,4.1,MATLAB,与多,项式,一、多项式的建立,1,、,直接输入系数向量：,把多项式的系数按降幂从左至右排列,，,存成,行向量,例：建立多项式,p=2,5,0,4,1,4,;,注：以上命令得到的是多项式的数值形式,补充：,将多项式的数值形式转化为字符串或符号表达式的形式,poly2str(p,s),将多项式,p,转化为以,s,为变量的字符串形式,poly2sym(p),将多项式,p,转化为符号表达式,poly2sym(p,s),将,多项式,p,转化为以,s,为变量的符号表达式,p=2,5,0,4,1,4,;,A=poly2sym(p),B=poly2str(p,s,),A=,2*x5+5*x4+4*x2+x+4,B=,2 s5+5 s4+4 s2+s+4,2,、,用多项式的全部根,X,来生成,poly(X),该函数返回以,X,为全部根的一个多项式,P,（,首项系数为,1,），当,X,是一个长度为,m,的向量时，,P,是一个长度为,m+1,的向量。,3,、,用,polyfit,生成,多项式的系数,-,用于多项式插值,给定,n,+1,个点可以唯一确定一个,n,阶多项式。,调用格式为：,p,=polyfit(,x,y,n,),其中,x,y,是同维向量，代表,n+1,个数据点的横、纵坐标，,n,是多项式的阶数。,二、多项式的运算,1,、求根运算,roots(,P,),P,是多项式,p,(,x,),的系数向量，该函数返回,p,(,x,)=0,的全部根（含重根，复根）,2,、求值运算,polyval(,P,x,),求多项式,P,在某点或某些点的函数值；,若,x,为一数值，则求多项式,P,在该点处的值；,若,x,为向量或矩阵，则求多项式,P,在向量或矩阵中的每个元素处的值,polyvalm,(,P,x,),此处,x,必须,为方阵,，它以方阵为自变量求多项式,P,的值,例,已知一个多项式,(1),计算,f,(,x,)=0,的全部根,(2),由方程,f,(,x,)=0,的根构造一个多项式,g,(,x,),，并与,f,(,x,),进行对比,(3),计算,f,(5),、,f,(7.8),、,f,(9.6),、,f,(12.3),的值,P=3,0,4,-5,-7.2,5;,X=roots(P)%,求方程,f(x)=0,的根,G=poly(X)%,求多项式,g(x),与,f(x),的区别是首项系数为,1,X0=5,7.8,9.6,12.3;,f=polyval(P,X0)%,求多项式,f(x),在给定点的值,X1=magic(3);,ff=polyvalm(P,X1),3.,多项式的四则运算,(1),多项式的加减法,当两多项式向量的大小相等，直接两向量相加（减）,当两多项式次数不同时，低阶多项式必须用首,0,填补，使其与高阶多项式次数相同,P1=2,5,0,4,1,4;,P2=0,0,5,1,3,2;,P3=P1+P2,;,P=poly2str(P3,x),P=,2 x5+5 x4+5 x3+5 x2+4 x+6,(2),多项式的乘法,conv(P1,P2),求多项式,P1,和,P2,的乘积,(3),多项式的除法,Q,r=deconv(P1,P2),求,P1/P2,；,其中,Q,为,多项式,P1,除以,P2,的,商式,，,r,为,P1,除以,P2,的,余式,。这里，,Q,和,r,仍是多项式系数向量。,deconv,是,conv,的逆函数，即有,P1=conv(P2,Q)+r,例：,已知三个一次多项式的零点分别为,0.4,、,0.8,、,1.2,，求三个一次多项式的乘积。,法,1,：,X=0.4,0.8,1.2;c=poly(X),P=poly2sym(c),运行后的结果：,c=,1.0000 -2.4000 1.7600 -0.3840,P=,x3-12/5*x2+44/25*x-48/125,法,2,：,P1=poly(0.4);P2=poly(0.8);P3=poly(1.2);,c=conv(conv(P1,P2),P3),P=poly2sym(c),运行后的结果：,c=,1.0000 -2.4000 1.7600 -0.3840,P=,x3-12/5*x2+44/25*x-48/125,例,设有两个多项式，计算：,(1),求,f(x)+g(x),、,f(x)-g(x),。,(2),求,f(x)g(x),、,f(x)/g(x),。,f=3,-5,2,-7,5,6;g=3,5,-3;g1=0,0,0,3,5,-3;,f+g1%,求,f(x)+g(x),f-g1,%,求,f(x)-g(x),conv(f,g)%,求,f(x)*g(x),Q,r=deconv(f,g)%,求,f(x)/g(x),，商式给,Q,，余式给,r,。,4,、多项式的微分与积分,（,1,）对多项式,求导数,：,p=polyder(P),求多项式,P,的导函数,p=polyder(P,Q),求,P*Q,的导函数,p,q=polyder(P,Q),求,P/Q,的导函数，导函数的,分子存入,p,，分母存入,q,（,p,q,为多项式的系数）。,（,2,）对多项式,求积分,：,d=polyint(P),d,是多项式,P,积分后的系数，但不包括 积分常数,例,求有理分式的导数,P=3,5,0,-8,1,-5;,Q=10,5,0,0,6,0,0,7,-1,0,-100;,p,q=polyder(P,Q),若求多项式,P,的积分：,c=polyint(P),4.2,MATLAB,插,值,我们经常会遇到，对于,y=f(x),它是存在的,，,可并不知道它的具体函数表达式，只知道离散数据,希望找到一个简单函数,且满足,被插函数,插值函数（常取多项式）,插值节点,插值条件：,插值区间：,的余项或误差：,若找到,P,(,x,)=,a,0,+,a,1,x,+,a,n,x,n,，则称相应的插值为,代数插值（或多项式插值）,定理,：,若插值节点,x,1,x,2,x,n+1,是,(,n,+1),个,互异,点，则满足插值条件,P,(,x,k,)=,y,k,(,k,=0,1,n,),的,n,次插值多项式,P,(,x,)=,a,0,+,a,1,x,+,a,n,x,n,存在且唯一。,此定理表明：在,n+1,个互异节点构造出来的,n,次多项式是唯一的，也就是说，后面将要学到的幂级数插值多项式、,lagrange,插值多项式、牛顿插值多项式是恒等的。（,n,次插值多项式是唯一的,）,证明,:,由插值条件,P,(,x,1,)=,y,1,P,(,x,2,)=,y,2,P,(,x,n+1,)=,y,n+1,系数行列式,故方程组有唯一解,.,从而插值多项式,P,(,x,),存在而且是唯一的,.,P,(,x,)=,a,0,+,a,1,x,+,a,n,x,n,例 已知误差函数在四个点处函数值,x,00.60001.20001.8000,Erf,(,x,),00.60390.91030.9891,构造,3,次多项式,P,(,x,),逼近,Erf,(,x,),设,P,(,x,)=,a,0,+,a,1,x,+,a,2,x,2,+,a,3,x,3,令,P,(,x,i,)=,Erf,(,x,i,),得,求解，得,a,0,=0,a,1,=1.293,a,2,=-0.5099,a,3,=0.0538,所以,P,(,x,)=1.293,x,0.5099,x,2,+0.0538,x,3,程序一：,V=1.8 1.2 0.6 0;,A=rot90(vander(V),2),b=0 0.6039 0.9103 0.9891;,a=(Ab);%,四元线性方程组的解,bb=rot90(a,2)%,将解向量的元素降幂排列,px=poly2str(bb,x),程序二,:,x=0:0.6:1.8;y=0 0.6039 0.9103 0.9891;,x=x;y=y;,A=ones(4,1)x x.2 x.3;,p=Ay;,a0=p(1);a1=p(2);,a2=p(3);a3=p(4);,%,作插值函数和离散点的图形,t=0:0.2:2;,u=a0+a1*t+a2*t.2+a3*t.3;,plot(x,y,rp,t,u),一、线性插值,(,或一次多项式插值）,线性插值是两个数据点的直线拟合,或,误差,(,余项,),估计：,“点斜式”,“两点式”,在,MATLAB,中，命令,interp1,可做线性插值，,yi,interp1,(,x,y,xi,),其中：,x,表示已知数据点,横坐标的,向量；,y,表示已知数据点,纵坐标,的向量；,xi,为插值点的横坐标，输出的,yi,为插值点的函数值,例,已知数据表如下，分别求,x,0.9,，,0.7,，,0.6,，,0.5,处,y,的值。,x,y,0.9126,0,0.8109,0.25,0.6931,0.50,0.5596,0.75,0.4055,1.00,x=0.9126,0.8109,0.6931,0.5596,0.4055;,y=0.0,0.25,0.5,0.75,1.0;,xi=0.9,0.7,0.6,0.5;,yi=interp1(x,y,xi);,xi,y,i,ans,0.9000 0.0310,0.7000 0.4854,0.6000 0.6743,0.5000 0.8467,另外，,interp1,命令有,5,种可选参数,yi,=,interp1,(,x,y,xi,linear,),线性插值（缺省）,yi,=,interp1,(,x,y,xi,nearst,),最近邻点插值,yi,=,interp1,(,x,y,xi,pchip,),分段三次,Hermite,插值,yi,=,interp1,(,x,y,xi,cubic,),同上,yi,=,interp1(x,y,xi,spline),分段三次样条插值（建议使用三次样条插值）,若已知的数据点,x,是等距节点，插值速度可提高，此时，填写最后一个插值方法时，以星号引导，即为,*,linear,*spline,*cubic,*pchip,*nearst,例：,假设已知的数据点来自函数,试根据生成的数据进行插值,得到较光滑的曲线,解：,x=0:0.12:1;,y=(x.2-3*x+5).*exp(-5*x).*sin(x);,plot(x,y,-o),可看出，由这样的数据直接连线绘制出来的曲线十分粗糙，可再选择一组插值点，然后用,interp1,插值,x=0:,0.12,:1;,y=(x.2-3*x+5).*exp(-5*x).*sin(x);,x0=0:,0.02,:1;,y0=(x0.2-3*x0+5).*exp(-5*x0).*sin(x0);,y01=interp1(x,y,x0);,y02=interp1(x,y,x0,cubic);,y03=interp1(x,y,x0,spline);,y04=interp1(x,y,x0,nearest);,subplot(1,3,1);,plot(x,y,mo,x0,y0,k);,legend(0.12,0.02);,subplot(1,3,2);,plot(x0,y01,r:,x0,y02,g-);,legend(,line,cubic,);,subplot(1,3,3);,plot(x0,y03,r-,x0,y04,b:);,legend(,spline,nearest,),interp1,默认的,linear,插值得到的曲线和,plot,画出来的曲线一样粗糙,而,nearest,选项得到的插值效果就更差了。,而采用,cubic,和,spline,选项得到的插值更接近于理论值。,事实上，采用样条插值算法得到的插值十分逼近理论值，甚至用肉眼难以分辨。,二、用幂级数做多项式插值,给定,n+,1,个数据点：,过,n,+1,个点的,n,阶多项式可写为幂级数形式：,注：过,n,1,个数据点的,n,阶插值多项式是唯一的。,对,n,1,个数据点，设 ，则得到,n,+1,个线性方程，可以表示为矩阵形式,求解该方程组可确定系数（或用,polyfit(,x,y,n,),确定）,例 求用下列数据点拟合出的多项式的系数，并求当,x,取,2.101,、,4.234,时,y,的值，并画出已知数据点和拟合曲线。,clear,x=1.1,2.3,3.9,5.1;,y=3.887,4.276,4.651,2.117;,n=length(x)-1;,c=polyfit(x,y,n);,xi=2.101,4.234;,yi=polyval(c,xi),plot(x,y,ro),hold on,xp=1.1:0.05:5.1;,yp=polyval(c,xp);,pause,plot(xp,yp),yi=,4.1457 4.3007,三、,Lagrange,插值多项式,1,、,给定,n+,1,个数据点：,插值基函数：,满足插值条件：,的次数不超过,n,的,lagrange,多项式,g(x),为：,function y=lagrange1(x0,y0,x),np=length(x0);,y=zeros(size(x);,for i=1:np,z=ones(size(x);,for j=1:np,if j=i,z=z.*(x-x0(j)/(x0(i)-x0(j);,end,end,y=y+y0(i)*z,end,2,、,MATLAB,程,序（法,1,）,:,调用格式：,y=lagrange1(x0,y0,x),clear,x0=1.1,2.3,3.9,5.1;,y0=3.887,4.276,4.651,2.117;,x=2.101,4.234;,y=lagrange1(x0,y0,x),例：写出拟合下面四个数据点的,Lagrange,插值公式，并计算,x,2.101,、,4.234,时,y,的值。,x0,y0,1.1,3.887,2.3,4.276,3.9,4.651,5.1,2.117,结果为,y=,4.1457 4.3007,由于,lagrange,插值基函数,l,i,(x),在节点,x,k,(k=1,2,n),处满足,又由于过,n+1,个数据点的,n,次插值多项式是唯一的,故,l,i,(x),也可看做是拟合下列数据点的,n,次多项式,可编写程序求出,2,、,MATLAB,程,序实现（法,2,）,进而求出,lagrange,插值公式,function p=lag_base(x),np=length(x);,for j=1:np,y=zeros(1,np);y(j)=1;,p(j,:)=polyfit(x,y,np-1);,end,其调用格式为：,p=lag_base(x),其中,:,x,是数据点的横坐标数组，,p,是一个矩阵，它的第,i,行,为,l,i,(x),的幂系数。,先编写程序求出,如前例也可求解如下,：,x=1.1,2.3,3.9,5.1;,y=3.887,4.276,4.651,2.117;,xi=2.101,4.234;,p=lag_base(x);,np=length(x);,s=0;,for i=1:np,s=s+p(i,:).*y(i);,end,s,yi=polyval(s,xi),结果为：,yi,4.1457 4.3007,3,、,Lagrange,插值公式的微分,插值公式,微分,为计算,Lagrange,插值多项式的一阶导数，可用,polyder,函数将,p,的每一行转换为一阶导数的系数数组。,例：求出拟合下面四个数据点的,Lagrange,插值多项式在,x,2.101,、,4.234,处的一阶导数值。,x0,y0,1.1,3.887,2.3,4.276,3.9,4.651,5.1,2.117,x=1.1,2.3,3.9,5.1;,y=3.887,4.276,4.651,2.117;,xi=2.101,4.234;,np=length(x);,p=lag_base(x);,s=0;,for i=1:np,s=,s+,polyder(p(i,:),.*y(i);,end,yi=polyval(s,xi),结果为：,yi,0.6292 -1.4004,例,考虑一个著名的例子,假设已知其中一些点的坐标，则可以采用下面的命令进行,Lagrange,插值，得出相应的插值曲线,x0=-1:0.2:1;,y0=1./(1+25*x0.2);,x=-1:0.01:1;,y=lagrange1(x0,y0,x);,ya=1./(1+25*x.2);,plot(x,y,b:p,x,ya,r),4,、,Runge,（龙格）现象,可看到：由,lagrange,插值得出的效果和精确值相差甚远，这种多项式阶次越高越发散的现象称为,Runge,现象,所以对此例来说，传统的,lagrange,算法失效,下面考虑用,interp1,来解决此问题,x0=-1:0.2:1;,y0=1./(1+25*x0.2);,x=-1:0.01:1;,ya=1./(1+25*x.2);,y1=interp1(x0,y0,x,cubic);,y2=interp1(x0,y0,x,spline);,plot(x,ya,x,y1,:,x,y2,-),interp1,函数效果,可见：,Matlab,中提供的算法不存在,Runge,现象，可大胆放心使用,总结：,（,1,）尽可能,在小区间上,使用多项式插值；,（,2,）只能,在一定范围内,依靠增加插值点个数提高插,值精度，如果插值点个数过多往往会适得其反。,取,x,0,x,1,x,2,求二次插值函数,P,(,x,)=,a,0,+,a,1,(,x x,0,)+,a,2,(,x,x,0,)(,x,x,1,),满足插值条件,P,(,x,0,)=,f,(,x,0,),P,(,x,1,)=,f,(,x,1,),P,(,x,2,)=,f,(,x,2,),引出关于,a,0,a,1,a,2,方程：,四、牛顿插值问题,解上述下三角方程组过程中引入符号,a,0,=,f,(,x,0,),a,1,=,f,x,0,x,1,a,2,=,f,x,0,x,1,x,2,P,(,x,)=,a,0,+,a,1,(,x x,0,)+,a,2,(,x,x,0,)(,x,x,1,),定义,:,若已知函数,f,(,x,),在点,x,1,x,2,x,n+1,处的函数值为,f,(,x,1,),f,(,x,2,),f,(,x,n+1,),如果,i,j,则,n,阶均差,(,j=,1,2,n,),一阶均差,二阶均差,(,j=1,2,n,-1,),x,-2-1013,y,-56-16-2-24,例,由函数表求,各阶均差,解,:,按公式计算各阶差商如下,xf,(,x,),一阶 二阶 三阶 四阶,差商 差商 差商 差商,-2,-56,-1,-16,40,0-,2 14-13,1,-2 0-7 2,3,4 3 1 2 0,计算均差矩阵的,程序：,x=-2-1 0 1 3;,y=-56-16-2-2 4;,n=length(x);,A=zeros(n,n);,A(:,1)=y;,for j=2:n,for i=j:n,A(i,j)=(A(i,j-1)-A(i-1,j-1)/(x(i)-x(i-j+1);,end,end,A,-56 0 0 0 0,-16 40 0 0 0,-2 14 -13 0 0,-2 0 -7 2 0,4 3 1 2 0,xf,(,x,),一阶 二阶 三阶 四阶,差商 差商 差商 差商,-2,-56,-1,-16,40,0-,2 14-13,1,-2 0-7 2,3,4 3 1 2 0,有了均差，就能得到牛顿插值公式（均差插值公式）：,以,3,阶牛顿插值公式为例，改写形式,在编写牛顿插值多项式程序之前，复习几个命令：,例,求三个一次多项式积，它们的零点分别依次为,0.4,、,0.8,、,1.2,。,解,:X1=0.4,0.8,1.2;,l1=poly(X1),L1=poly2sym(l1),运行后输出结果为,l1=,1.0000 -2.4000 1.7600 -0.3840,L1=,x3-12/5*x2+44/25*x-48/125,P1=poly(0.4);P2=poly(0.8);P3=poly(1.2);,C=conv(conv(P1,P2),P3),L1=poly2sym(C),运行后输出的结果与方法,1,相同。,function A,C,N=newton_poly(x,y),%x,为已知点的横坐标，,y,为已知点的纵坐标,%A,为均差矩阵，,C,为牛顿插值多项式的系数向量,%N,为牛顿插值多项式,%,先产生均差矩阵,A,n=length(x);A=zeros(n,n);A(:,1)=y;,for j=2:n,for i=j:n,A(i,j)=(A(i,j-1)-A(i-1,j-1)/(x(i)-x(i-j+1);,end,end,牛顿插值的,Matlab,算法：,%,下面求牛顿插值多项式系数向量,C,C=A(n,n);%3,阶牛顿插值公式，,n=4,for k=(n-1):-1:1,C=conv(C,poly(x(k);,d=length(C);,C(d)=C(d)+A(k,k);,end,%,最后求牛顿插值多项式,N,N=poly2sym(C);,例：,给出节点,x=-2.15 -1.00 0.01 1.02 2.03 3.25,，,y=17.03 7.24 1.05 2.03 17.06 23.05,，作五阶牛顿插值多项式和差商。,x=-2.15 -1.00 0.01 1.02 2.03 3.25;,y=17.03 7.24 1.05 2.03 17.06 23.05;,A,C,N=newton_poly(x,y),结果：,A=,17.0300 0 0 0 0 0,7.2400 -8.5130 0 0 0 0,1.0500 -6.1287 1.1039 0 0 0,2.0300 0.9703 3.5144 0.7604 0 0,17.0600 14.8812 6.8866 1.1129 0.0843 0,23.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.2169,C=,-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954,N=,-0.21686 x5+0.06482 x4+2.1076 x3+3.396 x2-4.5745 x +1.0954,五、三次,Hermite,插值,已知,构造,满足,即拟合两点函数值和导数值。,求解,即可,例,:,已知插值条件,:,求,3,次插值函数,.,解,:,设,得,a,0,=0,a,1,=0,列出方程组,求解,得,a,2,=3,a,3,=2,所以,有,H,(x)=3,x,2,2,x,3,=(3 2,x,),x,2,x,0,1,H,(,x,)0 1,H,(,x,)0 0,六、分段插值,问题：,节点增多，多项式次数增高，逼近精度越,好？未必！多节点高次插值往往在局部误,差更大,Runge,现象。,实用：,采用分段低次插值,有分段线性，分段二次插值，分段三次插值等,分段插值法可以克服在大范围内使用高次插值带来的误差放大的问题。,所谓的,分段插值法,就是首先在插值区间,a,b,上插入节点 ，然后在每个小的子区间 上构造低次插值多项式。再将每个子区间上的多项式连接，作为插值区间的插值函数。,如：分段线性插值：区间,a,b,上的分段线性插值函数为,分段线性插值在,MATLAB,中的函数调用：,1,）横纵坐标都是向量；,2,）,x,是,n,维向量，,y,为矩阵时，行数必须是,x,的维数。,例,：,x0=-6:6;y0=sin(x0);,xi=-6:.25:6;%,插值点,yi=,interp1,(x0,y0,xi);,x=-6:0.001:6;y=sin(x);plot(x0,y0,o,xi,yi,x,y,:),legend(,节点,(x_i,y_i),分段线性插值函数,被插值函数,y=sinx),title(y=sinx,及其分段线性插值函数和节点的图形,Color,w),分段,Hermite,插值,：即在每个小区间上用分段的,三次,Hermite,插值。,在,matlab,里有专门的调用函数：,y1=interp1(x,y,x1,pchip),例,:x0=-6:6;y0=sin(x0);xi=-6:.25:6;,yi=interp1(x0,y0,xi,pchip);,x=-6:0.001:6;y=sin(x);,plot(x0,y0,o,xi,yi,x,y,:),legend(,节点,(xi,yi),分段埃尔米特插值函,数,被插值函数,y=sinx),title(y=sinx,及其分段埃尔米特插值函数和,节点的图形,Color,w),x=-5:5;,y=1./(1+x.2);,plot(x,y,x,y,o),x=-5:5;,y=1./(1+x.2);,xi=-5:.05:5;,yi=spline(x,y,xi);,plot(xi,yi,b,x,y,ro),被插值函数,:,-5,x,5,七、样条插值,分段线性插值只能保证曲线的连续性，不能保证曲线在连接点处的光滑性。,Hermite,插值虽然在连接点处一阶光滑，但整体插值由于节点多次数高而有可能发生,Runge,现象。,若用分段三次,Hermite,插值，节点处一阶导数连续光滑，但在曲线的凹凸性变化比较大的地方，误差较大。,希望：,既想分段插值，又想在节点处保持光滑，甚至二阶光滑,-,三次样条,。,定义,：设有,n,+1,个点的点列,若函数 满足：,(3),在整个区间 内具有二阶连续导数。,(2),在每个小区间,上是三次多项式；,-,此时 叫插值函数；,则称 为点列的,三次样条插值函数,或,三次样条多,项式,，简称,三次样条,。,定义中的一阶导数连续意味着曲线没有急弯，二,阶导数连续意味着曲线每一点的曲率半径有意义。,在,MATLAB,中，三次样条插值运算实现如下：,yi=interp1(,x,y,xi,spline,),或,yi=spline(,x,y,xi,),其中,x,，,y,都是向量形式的点，,xi,是进行插值的点的,横坐标向量，,yi,为插值函数值。,例,求三次样条插值并作图比较。,%,以,s,为参数，分别做,x,，,y,的样条函数,xx=-1-0.866-0.5 0 0.5 0.866 1 1.0402 1.1500 1.3.,1.54 1.8280 2.1736 2.5883 3.0860;,yy=0-0.25-0.433-0.5-0.433-0.25 0 0.15 0.2598 0.3.,0.3 0.3 0.3 0.3 0.3;,s=1:length(xx);sp=1:length(xx)/100:length(xx);,xp=spline(s,xx,sp);,yp=spline(s,yy,sp);%,以,s,为参数，分别对已知的两组数据,%,点进行三次样条插值,plot(xp,yp);hold on,plot(xx,yy,ro);xlabel(x);ylabel(y);,axis(-1 3.5-1 1),%,若不取,s,为参数，直接以,x,，,y,做样条函数，曲线不光滑。,xi=-1,:1/100:3.5;,yi=spline(xx,yy,xi);,plot(xi,yi,xx,yy,ro),x=0,0.0155,0.1485,0.3493,0.6480,1.0547,2.0;,y=0,0.1242,0.3654,0.4975,0.5472,0.4781,0;,求三次样条插值并作图比较。,上机练习：,n=length(x);,t=0:n-1;tt=0:.25:n-1;,xx=spline(t,x,tt);yy=spline(t,y,tt);,plot(xx,yy,x,y,o),总结：前面全都是一维插值，用,interp1,就可完成,若做二维插值，用,interp2,就可完成,若做三维插值，用,interp3,就可完成,二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函数,z=f(x,y),进行插值。为了说明这个附加的维数，考虑一个问题。设人们对海水的温度分布估计感兴趣，给定的温度值取自海水表面均匀分布的格栅。,采集了下列的数据：,width=1:5;%index for width of plate(i.e.,the x-dimension),depth=1:3;%index for depth of plate(i,e,the y-dimension),八、二维插值,temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86%temperature data,temps=,82 81 80 82 84,79 63 61 65 81,84 84 82 85 86,矩阵,temps,表示整个海水的温度分布。,temps,的列与下标,depth,或,y-,维相联系，行与下标,width,或,x-,维相联系。为了估计在中间点的温度，我们就需要用到二维插值。,MESHGRID,命令：,调用格式一：,X,Y=meshgrid(x,y),X,Y=meshgrid(-3:.2:3,-3:.2:3);,Z=7-3*X.4.*exp(-X.2-Y.2),mesh(Z),title(Z=7-3 X4 exp(-X2-Y2),的图形,),调用格式二：,X,Y=meshgrid(x),调用格式三：,X,Y,Z=meshgrid(x,y,z),x=0,1,2,-4;y=0,-1,3,5;z=y;,X,Y,Z=meshgrid(x,y,z),U=2+X.*exp(-X.2-Y.2-Z.2),插值方法：,双线性插值,用两次线性插值对二维函数的数据表格做插值，二维表格是矩形网格点 上的函数值阵列,设需估计如图所示矩形区域,上一些点的函数值。,在,y,方向上做线性插值，则,E,、,F,处的值分别为：,然后做 的线性插值,将这两步结合为一个方程，即为双线性插值：,在,MATLAB,中，双线性插值命令为,table2,g,=table2(tab,x,y,),其中,tab,为一二维数据表格，第一列是 值的数组，第一行是 值的数组，按升序排列，剩余行列为 ，,x,，,y,是待插值点的横、纵坐标值，可以是标量、向量或矩阵。,2.,双,Lagrange,插值,在两个维度上利用两次,Lagrange,插值的方法。,其中形状函数,注：,当,M,N,2,时，方程化简为双线性插值表达式。,三次样条插值：给定区间,a,b,上的一个分划,a=x,0,x,1,x,n,=b,已知,f,(,x,j,)=,y,j,(,j=,0,1,n,),如果,满足,:(1),S,(,x,),在,x,j,，,x,j+,1,上为三次多项式,;,(2),S”,(,x,),在区间,a,，,b,上连续,;,(3),S,(,x,j,)=,y,j,(,j=,0,1,n,).,则称,S,(,x,),为三次样条插值函数,.,当,x,x,j,x,j+,1,(,j=,0,1,n-,1),时,S,j,(,x,)=,a,j,+b,j,x+c,j,x,2,+,d,j,x,3,插值条件,:,S,(,x,j,)=,y,j,(,j=,0,1,n,),连续性条件,:,S,(,x,j,+,0)=S(,x,j,-,0),(,j=,1,n-,1),S,(,x,j,+,0)=S(,x,j,-,0),(,j=,1,n-,1),S”,(,x,j,+,0)=S”(,x,j,-,0),(,j=,1,n-,1),由样条定义,可建立方程,(4,n,-2),个！即：,方程数少于未知数个数？,n,个三次多项式,待定系数共,4,n,个！,(1),自然边界条件,:,S”,(,x,0,)=0,S”,(,x,n,)=0,例,5.7,已知,f,(1)=1,f,(0)=0,f,(1)=1.,求,1,，,1,上的三次自然样条,(,满足自然边界条件,).,解,:,设,则有,:,a,1,+,b,1,c,1,+,d,1,=1,a,2,+,b,2,+,c,2,+,d,2,=1,d,1,=d,2,=0,c,1,=c,2,b,1,=b,2,(2),周期边界条件,:,S,(,x,0,)=,S,(,x,n,),S”,(,x,0,)=,S”,(,x,n,),(3),固定边界条件,:,S,(,x,0,)=,f,(,x,0,),S,(,x,n,)=,f,(,x,n,),由自然边界条件,:6,a,1,+2,b,1,=0,6,a,2,+2,b,2,=0,解方程组,得,a,1,=-,a,2,=,1/2,b,1,=b,2,=3/2,c,1,=c,2,=d,1,=d,2,=,0,问题的解,x=-1,0,1;y=1,0,1;,f1=inline(0.5*x.3+1.5*x.2);,f2=inline(-0.5*x.3+1.5*x.2);,t1=-1:.1:0;t2=0:.1:1;,p1=f1(t1);p2=f2(t2);,plot(x,y,o,t1,t2,p1,p2,r),y=x,2,