数值计算（一）主讲：张森2011-7-9一、矩阵的数值计算相关MATLAB函数提示：二、插值法1、插值有关的MATLAB函数：2、拉格朗日和牛顿插值法（1） 拉格朗日多项式和基函数的MATLAB 程序 求拉格朗日插值多项式和基函数的MATLAB 主程序function [C, L,L1,l]=lagran1(X,Y) m=length(X); L=ones(m,m); for k=1: m V=1;for i=1:m if k~=iV=conv(V,poly(X(i)))/(X(k)-X(i)); endendL1(k,:)=V; l(k,:)=poly2sym (V) endC=Y*L1;L=Y*l例1 给出节点数据03.17)15.2(=-f ，24.7)00.1(=-f ，05.1)01.0(=f ，03.2)02.1(=f ， 06.17)03.2(=f ，05.23)25.3(=f ，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.解 在MATLAB 工作窗口输入程序>> 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]; [C, L ,L1,l]= lagran1(X,Y)运行后输出五次拉格朗日插值多项式L 及其系数向量C ，基函数l 及其系数矩阵L 1如下C =-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x^2+2.1076*x^3+0.0648*x^4-0.2169*x^5 L1 =-0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048 -0.0693 0.2184 0.3961 -1.2116 -0.3166 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097 -0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002 l =[ -0.0056*x^5+0.0299*x^4-0.0323*x^3-0.0292*x^2+0.0382*x-0.0004] [ 0.0331*x^5-0.1377*x^4-0.0503*x^3+0.6305*x^2-0.4852*x+0.0048] [ -0.0693*x^5+0.2184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033] [ 0.0687*x^5-0.1469*x^4-0.5398*x^3+0.6528*x^2+0.9673*x-0.0097] [ -0.0317*x^5+0.0358*x^4+0.2530*x^3-0.0426*x^2-0.2257*x+0.0023] [ 0.0049*x^5+0.0004 *x^4-0.0266*x^3+0.0001*x^2+0.0220*x-0.0002] 估计其误差的公式为)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!6)()6(----++=x x x x x x fξ，)3.25,-2.15(∈ξ.拉格朗日插值及其误差估计的MATLAB 程序 拉格朗日插值及其误差估计的MATLAB 主程序function [y,R]=lagranzi(X,Y,x,M) n=length(X); m=length(x); for i=1:mz=x(i);s=0.0; for k=1:np=1.0; q1=1.0; c1=1.0;for j=1:nif j~=kp=p*(z-X(j))/(X(k)-X(j)); endq1=abs(q1*(z-X(j)));c1=c1*j; ends=p*Y(k)+s; endy(i)=s; endR=M*q1/c1;例 2 已知5.030sin = ，1707.045sin = ，0866.060sin = ，用拉格朗日插值及其误差估计的MATLAB 主程序求 40sin 的近似值，并估计其误差.解 在MATLAB 工作窗口输入程序>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; [y,R]=lagranzi(X,Y,x,M)运行后输出插值y 及其误差限R 为y = R =0.6434 8.8610e-004.（2） 牛顿插值法的MATLAB 综合程序求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X); m=length(x); for t=1:mz=x(t); A=zeros(n,n);A(:,1)=Y';s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j; endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); endy(k)= polyval(C, z); endR=M*q1/c1;L(k,:)=poly2sym(C);例3 将区间 [0,π/2] 分成n 等份)3,2(=n ，用x x f y s i n )(==产生1+n 个节点，求二阶和三阶牛顿插值多项式)(2x P 和)(3x P .用它们分别计算sin (π/7) (取四位有效数字)，并估计其误差.解 首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7; X2=0:pi/6:pi/2; Y2 =sin(X2);[y1,R1,A1,C1,P2]=newdscg(X1,Y1,x,M),[y2,R2,A2,C2,P3]=newdscg(X2,Y2,x,M)运行后输出插值y 1=)7/sin()7/(2π≈πP 和y 2=)7/sin()7/(3π≈πP 及其误差限R 1和R 2，二阶和三阶牛顿插值多项式P 2和P 3及其系数向量C 1和C 2，差商的矩阵A 1和A 2如下y1 =0.4548 R1 =0.0282 A1 =0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357 C1 =-0.3357 1.1640 0 P2 =-3024156947890437/9007199254740992*x^2+163820246322191/140737488355328*xy2 =0.4345 R2 =9.3913e-004 A2 =0 0 0 0 0.5000 0.9549 0 0 0.8660 0.6991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139 C2 =-0.1139 -0.0655 1.0204 0 P3 =-1025666884451963/9007199254740992*x^3-4717668559161127/72057594037927936*x^2+4595602396951547/4503599627370496*x3、高元插值及其MATLAB 程序（1） MESHGRID 命令的功能和调用格式调用格式一 [X,Y] =meshgrid (x,y)例1 已知x =-3:0.2:3；y=x ，计算函数437x z -=e 22y x --的值，并作出函数的图形.解 输入程序>> [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 X^4 exp(-X^2-Y^2)的图形')运行后输出函数值和图形（略）.例2 作出函数xz +=2e22yx --在区域22≤≤-x ，22≤≤-y 上的图形.解 输入程序>> [X,Y] = meshgrid(-2:.2:2, -2:.2:2);Z = 2+X .* exp(-X.^2 - Y.^2); mesh(Z)title('Z = 2+X exp(-X^2 - Y^2)的图形')运行后输出函数值和图形（略）例3 设节点),,(z y x 中的5:0.5:5- =x ,x y =和函数437x Z -=e 22y x --，作Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的双三次插值和二元最近邻插值及其图形.解 （1）双三次插值.输入程序>> [x,y] = meshgrid(-5:0.5:5);z =7-3* x.^4 .* exp(-x.^2 - y.^2); xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; [xi,yi] = meshgrid(xi,yi);zi = interp2(x,y,z,xi,yi, 'cubic'), mesh(xi,yi,zi) hold onplot3(x,y,z,'r.', 'markersize',3*5) hold offxlabel('x'), ylabel('y'), zlabel('z'),title('z =7-3 x^4 exp(-x^2 - y^2) 的双三次插值和数据点的图形')运行后屏幕显示Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的双三次插值及其图形(略).（2）二元最近邻插值.输入程序>> [x,y] = meshgrid(-5:0.5:5); z =7-3* x.^4 .*exp(-x.^2 - y.^2);xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; [xi,yi] =meshgrid(xi,yi);zi = interp2(x,y,z,xi,yi, 'nearest'), mesh(xi,yi,zi) hold on,plot3(x,y,z,'r.', 'markersize',3*5), holdoffxlabel('x'), ylabel('y'), zlabel('z')title('z =7-3 x^3 exp(-x^2-y^2) 的二元最近邻插值和数据点的图形')运行后屏幕显示Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的二元最近邻插值及其图形(略).（2) 三元插值及其MATLAB 程序例3 设节点),,(z y x 的坐标为y z y x =-=-=,15,3,0,1,12,1,0,4，计算函数x V +=2e 222zy x---在插值点13:25.0:3,3:25.0:3,10:25.0:3-=-=-=i i i z y x 处的三元线性插值，并作其图形.解 输入程序>> x=[-4,0,1,12];y=[-1,0,3,15];z=y; [X,Y,Z]= meshgrid(x,y,z);V= 2+X .* exp(-X.^2 - Y.^2- Z.^2); [xi,yi,zi] =meshgrid(-3:.25:10,-3:.25:3,-3:.25:13);vi = interp3(X,Y,Z,V,xi,yi,zi),slice(xi,yi,zi,vi,[-1 6 9.5],9,[-2 .2 9]), shading flat,lighting flatxlabel('x'), ylabel('y'), zlabel('z'),title('V=2+ x exp(-x^2 - y^2-z^2) 的三元线性插值图形') hold on,colorbar('horiz'), view([-30 45])运行后屏幕显示三元线性插值及其图形（略）.实例：1994年全国数学建模竞赛A 题练习题：（二维插值） 在一丘陵地带测量高程，x 和y 方向每隔100米测一个点，得高程数据如下。