《多项式与插值》PPT课件.ppt

《《多项式与插值》PPT课件.ppt》由会员分享，可在线阅读，更多相关《《多项式与插值》PPT课件.ppt（109页珍藏版）》请在三一办公上搜索。

1、,第四章 多项式与插值,4.1 MATLAB与多项式,一、多项式的建立1.MATLAB中多项式用行向量表示，其元素为 多项式的系数，且从左至右按降幂排列；2.已知一个多项式的全部根 X，求多项式系数的函数是poly(X)，该函数返回以 X为全部根的一个多项式 P（首项系数为1），当 X是一个长度为m的向量时，P是一个长度为 m+1的向量。,3.给定 n+1个点可以唯一确定一个 n 阶多项式，利 用polyfit可以确定多项式的系数。调用格式为：p=polyfit(x,y,n)其中 x,y是同维向量，代表数据点的横、纵坐标，n 是多项式的阶数。,二、多项式计算1.多项式求根 求多项式 p(x)的

2、根的函数是roots(P)，这里，P是 p(x)的系数向量，该函数返回方程 p(x)=0 的全部根(含重根，复根)。2.多项式求值 求多项式 p(x)在某点或某些点的函数值的函数是polyval(P,x)。若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。,例1 已知一个多项式(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

3、(X)%求多项式g(x)X0=5,7.8,9.6,12.3;f=polyval(P,X0)%求多项式f(x)在给定点的值,3.多项式的四则运算(1)多项式的加减法,注：多项式求值还有一个函数是polyvalm，其调用格式与polyval相同，但含义不同。polyvalm函数要求x为方阵，它以方阵为自变量求多项式的值。,function p3=poly_add(p1,p2)n1=length(p1);n2=length(p2);if n1=n2 p3=p1+p2;endif n1n2，p3=p1+zeros(1,n1-n2),p2;endif n1n2，p3=zeros(1,n2-n1),p1+

4、p2;end,加法：c=poly_add(a,b)减法:c=poly_add(a,-b),(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。,例2 设有两个多项式，计算：(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;po

5、ly_add(f,g)%求f(x)+g(x)poly_add(f,-g)%求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。（2）对多项式的积分函数：d=poly_itg(c)d是多项式c积分后的系数，但 不包括积分常数,function py=poly_itg(p)n=length(p

6、);py=p.*n:-1:1.(-1),0;,例3 求有理分式的导数。,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=poly_itg(P),4.2 MATLAB插值,通常取 为多项式函数代数插值（多项式插值）,已知f(x)在点xi上的函数值 yi=f(xi),(i=0,1,2,n),则称 P(x)为 f(x)的 n 次代数插值多项式.称x0,x1,xn为插值结点;P(x)为插值函数;条件P(xk)=yk(k=0,1,n)为插值条件;f(x)为被插值函数.,如果 P(x)=a0+a1x+an

7、xn满足:P(xk)=yk(k=0,1,n),设 f(x)C a,b,取点 a x0 x1xnb,代数插值问题,定理：若插值结点x0,x1,xn 是(n+1)个互异点，则满足插值条件 P(xk)=yk(k=0,1,n)的n次插值多项式 P(x)=a0+a1x+anxn存在而且是唯一的。,方程组系数矩阵取行列式,故方程组有唯一解.从而插值多项式P(x)存在而且是唯一的.,构造3次多项式P(x)逼近 Erf(x),设P(x)=a0+a1x+a2x2+a3x3,令 P(xk)=Erf(xk),得,求解，得a0=0,a1=1.293,a2=-0.5099,a3=0.0538所以,P(x)=1.293

8、x 0.5099 x2+0.0538 x3,MATLAB计算程序x=0:.6:1.8;y=erf(x);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:.2:2;u=a0+a1*t+a2*t.2+a3*t.3;plot(x,y,o,t,u),注:一次多项式插值-过两点直线;二次多项式插值-过三点抛物线;不用待定系数法-(1)计算量大;(2)不易讨论误差;,几何意义：两条曲线有交点（公共点）,一、线性插值,线性插值是两个数据点的直线拟合,或,误差估计：,在MATLAB中，命令 interp1可做线性插值，

9、调用格式为：yiinterp1(x,y,xi),其中 x 表示数据点横坐标的列数组，y 表示数据纵坐标的列数组（可以有多列）。,另外，interp1命令有三种可选参数,yi=interp1(x,y,xi,linear)线性插值（缺省）yi=interp1(x,y,xi,spline)三次样条yi=interp1(x,y,xi,cubic)三次插值,例3 已知数据表如下，分别求 y0.9，0.7，0.6，0.5 处 x 的值。,x=0.0,0.25,0.5,0.75,1.0;y=0.9126,0.8109,0.6931,0.5596,0.4055;yi=0.9,0.7,0.6,0.5xi=int

10、erp1(y,x,yi,linear);yi,xi,ans 0.9000 0.0310 0.7000 0.4854 0.6000 0.6743 0.5000 0.8467,二、用幂级数做多项式插值,给定 n+1 个数据点：,过 n+1个点的 n 阶多项式可写为幂级数形式：,注：过 n1个数据点的 n 阶插值多项式是唯一的。,对 n1个数据点，设，则得到 n+1个线性方程，可以表示为矩阵形式,求解该方程组可确定系数（或用 polyfit(x,y,n)确定）,例3 求下列数据点拟合多项式的系数，并求当x2.101和4.234处 y 的值，并画出数据点和曲线。,x=1.1,2.3,3.9,5.1;y

11、=3.887,4.276,4.651,2.117;n=length(x)-1;a(:,n+1)=ones(size(x);a(:,n)=x;for j=n-1:-1:1 a(:,j)=a(:,j+1).*x;endcoeff=ay;,xi=2.101,4.234;yi=zeros(size(xi);for k=1:n+1 yi=yi+coeff(k)*xi.(n+1-k)endxp=1.1:0.05:5.1;yp=zeros(size(xp);for k=1:n+1 yp=yp+coeff(k)*xp.(n+1-k);endplot(xp,yp,x,y,ro),三、Lagrange插值多项式,

12、1.插值基函数,形状函数的图形如下（n8）,2.Lagrange插值多项式,function fi=Lagran_(x,f,xi)fi=zeros(size(xi);np1=length(f);for i=1:np1 z=ones(size(xi);for j=1:np1 if i=j z=z.*(xi-x(j)/(x(i)-x(j);end end fi=fi+z*f(i);endreturn,3.MATLAB程序实现:,调用格式：yi=Lagran_(x,y,xi),clearx=1.1,2.3,3.9,5.1;y=3.887,4.276,4.651,2.117;xi=2.101,4.23

13、4;yi=Lagran_(x,y,xi),例4：写出拟合下面三个数据点的Lagrange插值公式，并计算 x2.101、4.234时 y 的值。,4.截断误差：,例5 用 的5个等距点对函数进行插值估计。,插值结果及误差分布如下图：,可见，误差峰值出现在端点附近的区间里，这是由于 的局部峰值在端点附近。,Runge反例:,(-5x5),取xk=5+k 计算:f(xk)(k=0,1,10)构造L10(x).取:tk=5+0.05k(k=0,1,200),计算:L10(tk),x=-5:5;y=1./(1+x.2);t=-5:0.05:5;y1=1./(1+t.2);n=length(t);for

14、 i=1:n z=t(i);s=0;for k=1:11 Lk=1;u=x(k);for j=1:11 if j=k,Lk=Lk*(z-x(j)/(u-x(j);end end s=s+Lk*y(k);end y2(i)=s;endplot(x,y,ko,t,y1,t,y2,r),减小误差的方案：（1）减小插值区域，即 ba；（2）增加数据点个数；（3）使用可变间距的数据点(Chebyshev点)。,总结：（1）尽可能在小区间上使用多项式插值；（2）只能在一定范围内依靠增加插值点个数提高插值精度，如果插值点个数过多往往会适得其反。,5.Lagrange 插值公式的微分与积分,插值公式,微分,实

15、际上，是一个拟合如下数据点的 n 阶多项式,一般地，是拟合数据点的多项式,所以，多项式 可通过拟合待定数据点的 n 阶多项式表示为幂级数形式。,对所有 i，的幂级数形式可用函数 shape_pw 计算,function p=shape_pw(x)np=length(x);for j=1:np y=zeros(1,np);y(j)=1;p(j,:)=polyfit(x,y,np-1);end,如例4也可求解如下：,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=shape_pw(x);s=0;fo

16、r i=1:np s=s+p(i,:).*y(i);endsyi=polyval(s,xi),结果为：yi 4.1457 4.3007,为计算Lagrange插值多项式的一阶导数，可用polyder函数将 p 的每一行转换为一阶导数的系数数组。,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=shape_pw(x);s=0;for i=1:np s=s+polyder(p(i,:).*y(i);endyi=polyval(s,xi),结果为：yi 0.6292-1.4004,取x0,x1,x2,

17、求二次函数 P(x)=a0+a1(x x0)+a2(x x0)(x x1)满足条件 P(x0)=f(x0),P(x1)=f(x1),P(x2)=f(x2),插值条件引出关于a0,a1,a2方程,四、牛顿插值问题,解下三角方程组过程中引入符号,a0=f(x0),a1=fx1,x2,a2=fx0,x1,x2,P(x)=a0+a1(x x0)+a2(x x0)(x x1),定义:若已知函数 f(x)在点 x0,x1,xn 处的值 f(x0),f(x1),f(xn),如果 i j,则,n阶均差,例 由函数表求各阶均差,解:按公式计算一阶差商、二阶差商、三阶差商如下,MATLAB程序计算,-56-16-

18、2-2 4-56 40 14 0 3-56 40-13-7 1-56 40-13 2 2-56 40-13 2 0,x=-2-1 0 1 3;y=-56-16-2-2 4;f=yn=length(x);for k=2:n for j=n:-1:k f(j)=(f(j)-f(j-1)/(x(j)-x(j+1-k);end fend,另外一个程序：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);en

19、dend 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,在编写牛顿插值多项式程序之前，复习几个命令：例 6.1.2 求三个一次多项式积。它们的零点分别依次为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.3840L1=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

20、=poly2sym(C)运行后输出的结果与方法1相同.,牛顿插值多项式的函数文件：function A,C,L=newploy(X,Y)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);endend C=A(n,n);for k=(n-1):-1:1 C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C);,例：给出节点X=-2.15-1.00 0.01 1.02 2.

21、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,L=newploy(X,Y),五、Chebyshev多项式,等距Lagrange插值的误差在中间部分最小，越靠近端点处误差越大，降低误差的一种方法是改变插值点的分布，即增大定义域中间部分的插值步长，同时减少定义域两端的插值步长，但是最优的插值点分布取决于插值多项式的目的。,若目的是近似一个函数，则最优插值点Chebyshev多项式的零点

22、，因为（1）此时L(x)分布是最平坦的分布；（2）误差值不会像等距插值那样随不同插值区间 而变化。,Lagrange插值余项定理:取插值结点ax0 x1xnb则满足Ln(xk)=f(xk)的 n 次插值多项式Ln(x)误差,其中,选择:x0,x1,xn,使,结论:Chebyshev多项式Tn+1(x)的全部零点,Chebyshev多项式：,当 时，Chebyshev多项式的图形如下,幂级数形式Chebyshev多项式的系数可由函数Cheby_ pw计算,function pn=Cheby_pw(n)pbb=1;if n=0,pn=pbb;break;endpb=1 0;if n=1,pn=pb

23、;break;endfor i=2:n;pn=2*pb,0-0,0,pbb;pbb=pb;pb=pn;end,调用格式：p=Cheby_pw(n)n为多项式阶数，p是系数的行数组,Chebyshev多项式的根系数可由 roots命令计算 roots(Cheby_ pw(n),也可如下求解，令,得,有 k 个根，都在-1,1，也可用下面公式计算：,若插值区间为a,b，可建立映射,下图比较了区间 上9个Chebyshev点和等距点条件下L(x)的分布,例.函数,取等距插值结点:-5,-4,-3,-2,-1,0,1,2,3,4,5,x-5,5,w11(x)=(x+5)(x+4)(x+3)(x+2)(

24、x+1)x(x-1)(x-2)(x-3)(x-4)(x-5),w11(x)的图形,-4.9491-4.5482-3.7787-2.7032-1.4087 0.0000 1.4087 2.7032 3.7787 4.5482 4.9491,w11(x),w11(x)=(x x0)(x x1)(x x2)(x x10),插值函数L10(x)取Chebyshev结点,插值函数L10(x)取等距结点插值,六、Lobatto点,Lobatto点是Chebyshev多项式的零点加上 x-1和 x1对区间a,b，k+1个Lobatto点可由 k 阶Chebyshev多项式得到,七、Legendre多项式,L

25、egendre多项式：,n 阶Legendre多项式的系数可由函数Legen_pw计算：,function pn=Legen_pw(n)pbb=1;if n=0,pn=pbb;break;endpb=1 0;if n=1,pn=pb;break;endfor i=2:n;pn=(2*i-1)*pb,0-(i-1)*0,0,pbb)/i;pbb=pb;pb=pn;end,调用：p=Legen_pw(n),n 阶Legendre多项式的根可由roots计算。,已知节点x0和x1处的函数值及导数值,求三次插值函数,八、三次Hermite插值问题,例:已知插值条件:求3次插值函数.,解:设,得 a0=

26、0,a1=0,列出方程组,求解,得 a2=3,a3=2,所以,有 H(x)=3x2 2x3=(3 2x)x2,利用基函数表示Hermite插值,定理:两点Hermite插值的误差估计式,反复应用Roll定理,得F(4)(t)有一个零点设为,由,得,所以,显然,F(t)有三个零点x0,x,x1,由Roll定理知,存在F(t)的两个零点t0,t1满足x0t0t1x1,而x0和x1也是F(x)的零点,故F(x)有四个相异零点.,当插值节点含有奇点时的处理方法，以三次Hermite多项式为例。三次Hermite多项式：,可以拟合两点函数值和导数值。,已知两点 和 及它们的函数值和导数值，可建立方程组：

27、,若用有限差分近似表示导数边界条件,可通过Lagrange插值确定该方程组的系数。,当y(x)、x(y)含有奇点时，常需引入一个参数 s，并将 x 和 y 表示为 x(s)和 y(s)，假设点A处 s=0，B处 s=1，则 x，y 表示为关于 s 的三次多项式：,x(s)和 y(s)的边界条件可写为：,其中a、b是任意参数，在一定程度上影响曲线形状。,例6 确定一条经过点A、B的曲线，使满足,解：引入 s，按上述方法可得到两方程组,分别求解：,a=3;b=3;c=0 0 0 1;0 0 1 0;1 1 1 1;3 2 1 01;a;4;0d=0 0 0 1;0 0 1 0;1 1 1 1;3

28、2 1 01;0;2;bs=0:0.01:1;x=polyval(c,s)y=polyval(d,s);plot(x,y),结果：c=-3 3 3 1 d=1 0 0 1,若用差分近似，设置四个数据点：,取，求解,z=0.01;a=3;b=3;s(1)=0;x(1)=1;y(1)=1;s(2)=z;x(2)=1+z*a;y(2)=1;s(3)=1-z;x(3)=4;y(3)=2-z*b;s(4)=1;x(4)=4;y(4)=2;,c=polyfit(s,x,length(s)-1);d=polyfit(s,y,length(s)-1);ss=0:0.1:1;xp=polyval(c,ss);y

29、p=polyval(d,ss);plot(xp,yp)ylabel(y)xlabel(x),结果是：c=-3.9021 3.1231 2.9691 1d=1.0307-0.0309 0.0002 1,九、分段插值问题：结点增多，多项式次数增高，逼近精度越 好？未必！多结点高次插值往往在局部误 差更大Runge（龙格）现象。实用：采用分段低次插值 有分段线性，分段二次插值等缺点：分段插值函数只能保证连续性，不能保证光滑性。,分段插值法可以克服在大范围内使用高次插值带来的误差放大的问题。所谓的分段插值法就是首先在插值区间a,b上插入节点然后在每个小的子区间xi-1,xi上构造低次插值多项式。再将每

30、个子区间上的多项式连接，作为插值区间的插值函数。1.分段线性插值：区间a,b上的分段线性插值函数为,分段线性插值在MATLAB中的函数调用：yi=interp1(x0,y0,xi)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(节点(xi,yi),分段线性插值函数,被 插值函数y=sinx)title(y=sinx及其分段线性插值函数和节点的图形),分段H

31、ermite插值：即在每个小区间上用分段的三次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及其分段埃尔米特插值函数和 节点的图形),x=-5:5;y=1./(1+x.2);plot(x,y,x,y,o),x=-5:5;y

32、=1./(1+x.2);xi=-5:.05:5;yi=spline(x,y,xi);plot(xi,yi,b,x,y,ro),被插值函数:,-5 x 5,十、样条插值,分段线性插值可以得到整体连续函数，但在连接点处一般不光滑，而Hermite 插值虽然在连接点处一阶光滑，但整体插值由于结点多次数高而有可能发生龙格现象，若用分段三次Hermite插值，节点处一阶导数连续光滑，但在曲线的凹凸性变化比较大的地方，误差较大。,既想分段插值，又想在结点处保持光滑，甚至二阶光滑三次样条。,希望：,定义：设有 n+1个点的点列,若函数 满足：,(3)在整个区间 内具有二阶连续导数。,(2)在每个小区间,上是

33、三次多项式；,-此时 叫插值函数；,则称 为点列的三次样条插值函数或三次样条多 项式，简称三次样条。,定义中的一阶导数连续意味着曲线没有急弯，二阶导数连续意味着曲线每一点的曲率半径有意义。,则样条函数可写为,考虑区间,令,列方程组,求解得：,所以,则在区间 上,同理，在区间 上,由在节点 处连续，式、相等，消去，得：,式除了对第一个点和最后一个点不成立外，其它各点均成立，所以共有 n-2 个方程，而待定的 个数为 n，所以还需两个额外条件，可由下述边界条件提供。,常见到的边界条件有三种：,（a）在端点处指定，若已知，则方程组为：,该方程组为含 n-2 个未知量的 n-2 个方程，其系数矩阵为三

34、对角阵，可用追赶法求解。,大多数情形中，可令，等价于在几何上趋于端点的曲线变为直线。,（b）从内部外推，利用 对 进行外推,将其代入方程组的第一个方程中，整理得：,类似，利用 对 进行外推,将其代入方程组的最后一个方程，整理得：,此时，同样得到一三对角方程组，求解即可。,（c）循环边界条件 适用于封闭曲线，即第一个点和最后一个点相 同，它们的导数值也相同，此时相当于有 n-1个点，相应有 n-1个方程，可以求解。,在MATLAB中，三次样条插值运算实现如下：yi=interp1(x,y,xi,spline)或 yi=spline(x,y,xi)其中 x，y 都是向量形式的点，xi是进行插值的点

35、的横坐标向量，yi 为插值函数值。,例7 求三次样条插值并作图比较。,%以s为参数，分别做x，y的样条函数xx=-1-0.866-0.5 0 0.5 0.866 1 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 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);plot(xp,y

36、p);hold onplot(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);

37、plot(xx,yy,x,y,o),二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函数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

38、 84 82 85 86%temperature datatemps=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=

39、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为一二维数据表格，第一列是

40、 值的数组，第一行是 值的数组，按升序排列，剩余行列为，x，y是待插值点的横、纵坐标值，可以是标量、向量或矩阵。,2.双Lagrange插值 在两个维度上利用两次Lagrange插值的方法。,其中形状函数,当MN2时，方程化简为双线性插值表达式。,三次样条插值：给定区间a,b上的一个分划 a=x0 x1 xn=b已知 f(xj)=yj(j=0,1,n),如果,满足:(1)S(x)在 xj，xj+1上为三次多项式;(2)S”(x)在区间a，b上连续;(3)S(xj)=yj(j=0,1,n).则称 S(x)为三次样条插值函数.,当xxj,xj+1(j=0,1,n-1)时 Sj(x)=aj+bj x

41、+cj x2+dj x3,插值条件:S(xj)=yj(j=0,1,n)连续性条件:S(xj+0)=S(xj-0)(j=1,n-1)S(xj+0)=S(xj-0)(j=1,n-1)S”(xj+0)=S”(xj-0)(j=1,n-1),由样条定义,可建立方程(4n-2)个！即：,方程数少于未知数个数？,n个三次多项式,待定系数共4n个！,(1)自然边界条件:S”(x0)=0,S”(xn)=0,例 5.7 已知f(1)=1,f(0)=0,f(1)=1.求1，1 上的三次自然样条(满足自然边界条件).,解:设,则有:a1+b1c1+d1=1,a2+b2+c2+d2=1 d1=d2=0,c1=c2,b1=b2,(2)周期边界条件:S(x0)=S(xn),S”(x0)=S”(xn)(3)固定边界条件:S(x0)=f(x0),S(xn)=f(xn),由自然边界条件:6a1+2b1=0,6a2+2b2=0,解方程组,得 a1=-a2=1/2,b1=b2=3/2,c1=c2=d1=d2=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=x2,