数值计算方法(2)

1、数值计算（一）主讲：张森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: mV=1;for i=1:mif 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 给出节点数据f(

2、-2.15)=17.03，f(-1.00)=7.24，f(0.01)=1.05，f(1.02)=2.03， f(2.03)=17.06，f(3.25)=23.05，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.解 在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及其系数矩阵L1如下C =-0.2169 0.0648 2.1076 3.3960

3、 -4.5745 1.0954 L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5 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.

4、0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002 l = -0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.1377*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0.2184*x4+0.3961*x3-1.2116*x2-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+

5、0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004 *x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002 估计其误差的公式为R5(x)=f(6)()6!(x+2.15)(x+1.00)(x-0.01)(x-1.02)(x-2.03)(x-3.25)，(-2.15,3.25). 2拉格朗日插值及其误差估计的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;

6、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 已知sin30 =0.5，sin45 =0.7071，sin60 =0.8660，用拉格朗日插值及其误差估计的MATLAB主程序求sin40 的近似值，并估计其误差.解 在MATLAB工作窗口输入程序>> x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.866

7、0; 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:nfor i=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1)/(X(i)-X

8、(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等份(n=2,3)，用y=f(x)=nisx产生n+1个节点，求二阶和三阶牛顿插值多项式P2(x)和P3(x).用它们分别计算sin (/7) (取四位有效数字)，并估计其误差.解 首先将名为

9、newdscg.m的程序保存为M文件，然后在MATLAB工作窗口输入程序>> 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)运行后输出插值y1=P2(/7)sin(/7)和y2=P3(/7)sin(/7)及其误差限R1和R2，二阶和三阶牛顿插值多项式P2和P3及其系数向量C1和C2，差商的矩阵A1和A2如下y1 =0.4548R1 =0.0282A1 =0 0

10、 00.7071 0.9003 01.0000 0.3729 -0.3357C1 =-0.3357 1.1640 0P2 =-3024156947890437/9007199254740992*x2+163820246322191/140737488355328*xy2 =0.4345R2 =9.3913e-004A2 =0 0 0 00.5000 0.9549 0 00.8660 0.6991 -0.2443 01.0000 0.2559 -0.4232 -0.1139C2 =-0.1139 -0.0655 1.0204 0P3 =-1025666884451963/900719925474

11、0992*x3-4717668559161127/72057594037927936*x2+4595602396951547/4503599627370496*x3、高元插值及其MATLAB程序（1） MESHGRID命令的功能和调用格式调用格式一 X,Y =meshgrid (x,y)例1 已知x=-3:0.2:3；y=x，计算函数z=7-3x4e-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 X

12、4 exp(-X2-Y2)的图形')运行后输出函数值和图形（略）.例2 作出函数z=2+xe-x-y2222在区域-2x2，-2y2上的图形. 解 输入程序>> X,Y = meshgrid(-2:.2:2, -2:.2:2);Z = 2+X .* exp(-X.2 - Y.2); mesh(Z)title('Z = 2+X exp(-X2 - Y2)的图形')运行后输出函数值和图形（略）例3 设节点(x,y,z)中的x= -5:0.5:5,y=x和函数Z=7-3x4e-x-y，作Z在插值点X=-3.9:0.5:5,Y=-4.9:0.5:4.5处的双三次插值

13、和二元最近邻插值及其图形.解 （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(&

14、#39;y'), zlabel('z'),title('z =7-3 x4 exp(-x2 - y2) 的双三次插值和数据点的图形')运行后屏幕显示Z在插值点X=-3.9:0.5:5,Y =-4.9:0.5:4.5处的双三次插值及其图形(略).（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);5 22zi = interp2(x,y,

15、z,xi,yi, 'nearest'), mesh(xi,yi,zi) hold on,plot3(x,y,z,'r.', 'markersize',3*5), hold offxlabel('x'), ylabel('y'), zlabel('z')title('z =7-3 x3 exp(-x2-y2) 的二元最近邻插值和数据点的图形')运行后屏幕显示Z在插值点X=-3.9:0.5:5,Y =-4.9:0.5:4.5处的二元最近邻插值及其图形(略).（2) 三元插值及其MATLA

16、B程序例3 设节点(x,y,z)的坐标为x=-4,0,1,12,y=-1,0,3,15,z=y，计算函数V=2+xe-x2-y-z22在插值点xi=-3:0.25:10,yi=-3:0.25:3,zi=-3:0.25:13处的三元线性插值，并作其图形.解 输入程序>> 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,z

17、i),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(-x2 - y2-z2) 的三元线性插值图形') hold on,colorbar('horiz'), view(-30 45)运行后屏幕显示三元线性插值及其图形（略）.实例：1994年全国数学建模竞赛A题练习题：（二维插值）在一丘陵地带测量高程，x和y方向每隔100米测一

18、个点，得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值，并由此找出最高点和该点三、函数逼近和曲线拟合1、多项式拟合下面的MATLAB程序给出了该函数的二次和三次拟合多项式。 x=-1:0.2:1;y=1/(1+25*x.*x);xx=-1:0.02:1;p2=polyfit(x,y,2);yy=polyval(p2,xx);plot(x,y,o,xx,yy);xlabel(x);ylabel(y);hold on;p3=polyfit(x,y,3);yy=polyval(p3,xx);plot(x,y,o,xx,yy);hold off;2、最佳平方逼近函数逼近及其MA

19、TLAB程序最佳均方逼近的MATLAB主程序function yy1,a,WE=zjjfbj(f,X,Y,xx)m=size(f);n=length(X);m=m(1);b=zeros(m,m); c=zeros(m,1);if n=length(Y)error('X和Y的维数应该相同')endfor j=1:mfor k=1:mb(j,k)=0;for i=1:nb(j,k)=b(j,k)+feval(f(j,:),X(i)*feval(f(k,:),X(i);endendc(j)=0;for i=1:nc(j)=c(j)+feval(f(j,:),X(i)*Y(i);end

20、enda=bc;WE=0;for i=1:nff=0;for j=1:mff=ff+a(j)*feval(f(j,:),X(i);endWE=WE+(Y(i)-ff)*(Y(i)-ff);endif nargin=3return;endyy=;for i=1:ml=;for j=1:length(xx)l=l,feval(f(i,:),xx(j);endyy=yy l'endyy=yy*a; yy1=yy' a=a'WE;例1 对数据X和Y，用函数y=1,y=x,y=x2，y=cosx，y=ex，y=sinx进行逼近，其中X=（0 0.50 1.00 1.50 2.00

21、 2.50 3.00），Y=（0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645）.解 在MATLAB工作窗口输入程序>> X= 0 0.50 1.00 1.50 2.00 2.50 3.00;Y=0 0.4794 0.8415 0.9815 0.9126 0.5985 0.1645;f='fun0''fun1''fun2''fun3''fun4''fun5'xx=0:0.2:3;yy,a,WE=zjjfbj(f,X,Y, xx), plot(X,Y,

22、'ro',xx,yy,'b-')运行后屏幕显示如下（图略）yy = Columns 1 through 7-0.0005 0.2037 0.3939 0.5656 0.7141 0.8348 0.9236Columns 8 through 140.9771 0.9926 0.9691 0.9069 0.8080 0.6766 0.5191Columns 15 through 160.3444 0.1642a = 0.3828 0.4070 -0.3901 0.0765 -0.4598 0.5653WE = 1.5769e-004即，最佳逼近函数为y=0.3828

23、+0.4070*x-0.3901*x2+0.0765*exp(x) -0.4598*cos(x) +0.5653*sin(x).3、最小二乘拟合四、数值积分与数值微分相关MATLAB函数提示：1 数值积分的思想及其MATLAB程序矩形公式的MATLAB程序（一） 函数sum的调用格式调用格式一：sum(X)调用格式二：sum (X,DIM)例1 用MATLAB和矩形公式（9.3）、（9.4）计算2esinxdx，并与精确值0比较.解 将0,/2分成20等份，步长为/40，输入程序>> h=pi/40; x=0:h:pi/2; y=exp(sin(x);z1=sum(y(1:20)*

24、h, z2=sum(y(2:21)*h,运行后屏幕显示矩形公式计算结果分别如下z1 = z2 =3.0364 3.1713求定积分的精确值，输入程序>> syms xF=int(exp(sin(x),x,0, pi/2), Fs= double (F),wz1=abs( Fs-z1), wz2= abs( Fs-z2)运行后屏幕显示定积分的精确值Fs和与用矩形公式（9.3），（9.4）计算结果的绝对误差wz1、wz2.（二） 函数cumsum的调用格式调用格式一：cumsum(X)调用格式二：cumsum (X,DIM)例2 用MATLAB的函数sum 和 cumsum及矩形公式计

25、算20e-xsinxdx，并与精确值比较.解 将0,/2分成20等份，步长为/40，输入程序如下（注意sum 和 cumsum的用法）>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x);z1=sum(y(1:20)*h,z2=sum(y(2:21)*h,z=cumsum(y); z11=z(20)*h, z12=(z(21)-z(1)*h,运行后屏幕显示计算结果分别如下z1 = z2 = z11 = z12 =0.3873 0.4036 0.3873 0.4036求定积分的精确值，输入程序>> syms xF=int(exp(-x)*si

26、n(x),x,0, pi/2)Fs= double (F) ,wz1=abs( Fs-z1), wz2= abs( Fs-z2)运行后屏幕显示定积分的精确值Fs和用矩形公式（9.3），（9.4）计算结果的绝对误差wz1、wz2分别如下F = Fs =1/2*(-1+exp(pi)(1/2)/exp(pi)(1/2) 0.3961wz1 = wz2 =0.0088 0.00752 插值型数值积分及其MATLAB 程序梯形公式的MATLAB程序(一) 根据梯形公式和估计误差公式自己编写MATLAB程序计算定积分例3 分别取h=/8000,/800,/80，用梯形公式计算定积分I=/20esinxd

27、x，并与精确值比较.然后观察h对计算结果的有效数字和绝对误差的影响.解 编写并输入如下程序>>h=pi/8000;a=0;b=pi/2;x=a:h:b;n=length(x),y=exp(sin(x);z1=(y(1)+y(n)*h/2; z2=sum(y(2:n-1)*h; z8000=z1+z2, syms tf=exp(sin(t); intf=int(f,t,a,b), Fs=double(intf),Juewucha8000=abs(z8000-Fs)运行后屏幕显示取h=/8000时，积分区间0,/2上等距节点的个数n，用梯形公式计算定积分I的值z8000和精确值intf

28、的近似值Fs及其绝对误差Juewucha8000.(二) 用函数trapz计算定积分调用格式一：Z =trapz(Y)调用格式二：Z =trapz(X,Y)调用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM)(三) 用函数cumtrapz计算定积分调用格式一：Z =cumtrapz (Y)调用格式二：Z =cumtrapz (X,Y)调用格式三：Z = cumtrapz (X,Y,DIM) 或 cumtrapz (Y,DIM) 例.4 用MATLAB的函数trapz 和cumtrapz分别计算/20e-xsinxdx，精确到10-4，并与矩形公式（9.3），（9.

29、4）比较.解 将0,/2分成20等份，步长为/40，输入程序如下（注意trapz(y)是单位步长, trapz(y)*h=trapz(x,y)）：>> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x);z1=sum(y(1:20)*h, z2=sum(y(2:21)*h, z=(z1+z2)/2 z3=trapz(y)*h, z3h=trapz(x,y), z3c=cumtrapz(y)*h, 运行后屏幕显示用矩形公式（9.3），（9.4）计算结果z1、z2和二者的平均数z、函数trapz 和cumtrapz分别计算结果z3、z3c.(四)梯形数值积分的

30、MATLAB主程序梯形数值积分的MATLAB主程序function T=rctrap(fun,a,b,m)n=1;h=b-a; T=zeros(1,m+1); x=a;T(1)=h*(feval(fun,a)+feval(fun,b)/2;for i=1:mh=h/2; n=2*n; s=0;for k=1:n/2x=a+h*(2*k-1); s=s+feval(fun,x);endT(i+1)=T(i)/2+h*s;endT=T(1:m);例5 用rctrap计算I=1220e-x22dx，递归14次，并将计算结果与精确值比较.解 输入程序>>T=rctrap(fun,0,pi/

31、2,14), syms tfi=int(exp(-t2)/2)/(sqrt(2*pi),t,0, pi/2); Fs= double(fi), wT= double(abs(fi-T)运行后屏幕显示I精确值Fs，用rctrap计算I的递归值T和T与精确值Fs的绝对误差wT3 辛普森公式及其误差分析例6 用辛普森公式计算I=1210e-x22dx，取n=20001个等距节点，并将计算结果与精确值比较，然后再取n=13计算，观察n对误差的影响.解 由n=2m+1=20001，得m=100（Simpson）公式编写并00.根据辛普森输入下面的程序>> a=0;b=1;m=10000; h

32、=(b-a)/(2*m); x=a:h:b;y=exp(-x.2)./2)./(sqrt(2*pi);z1=y(1)+y(2*m+1); z2=2*sum(y(2:2:2*m);z3=4*sum(y(3:2:2*m);z=(z1+z2+z3)*h/3, syms t,f=exp(-t2)/2)/(sqrt(2*pi);intf=int(f,t,a,b), Fs=double(intf); Juewucha=abs(z-Fs)运行后屏幕显示用辛普森公式(9.11)计算定积分I的近似值z和精确值intf及其绝对误差Juewucha（取n=20001个等距节点）.4 辛普森（Simpson）数值积分

33、的MATLAB程序调用格式一：quad(fun,a,b)调用格式二：quad(fun,a,b,tol)调用格式三：Q,FCNT = quad (.)调用格式四：quad(fun,a,b, tol,TRACE)调用格式五：quad(fun,a,b, tol,TRACE,P1,P2, )复合辛普森（Simpson）数值积分的MATLAB主程序function y=comsimpson(fun,a,b,n)z1=feval (fun,a)+ feval (fun,b);m=n/2;h=(b-a)/(2*m); x=a;z2=0; z3=0; x2=0; x3=0;for k=2:2:2*mx2=x+

34、k*h; z2= z2+2*feval (fun,x2);endfor k=3:2:2*mx3=x+k*h; z3= z3+4*feval (fun,x3);endy=(z1+z2+z3)*h/3;例8 用comsimpson.m和quad.m分别计算定积分I=1210e-x22dx，取精度为10-4，并与精确值比较.解 输入程序>> Q1,FCNT14 = quad(fun,0,1,1.e-4,3),Q2 =comsimpson (fun,0,1,10000)syms xfi=int(exp( (-x.2)./2)./(sqrt(2*pi),x,0, 1); Fs= double

35、 (fi)wQ1= double (abs(fi-Q1) ), wQ2= double (abs(fi-Q2) ) 运行后屏幕显示I的精确值Fs，用comsimpson.m和quad.m分别计算I的近似值Q2、Q1和迭代次数FCNT14，取精度分别为10-4，Q2、Q1分别与精确值Fs的绝对误差wQ2, wQ1如下9 0.0000000000 2.71580000e-0010.107027510011 0.2715800000 4.56840000e-001 120.159794224213 0.7284200000 2.71580000e-001 0.0745230082Q1 = FCNT1

36、4 = Q2 =0.3413 13 0.3413Fs = wQ1 = wQ2 =0.3413 3.6619e-009 3.7061e-005五、微分方程数值解问题提出：考虑著名的Lorenz方程dxdt=s(y-x)dy =rx-y-xzdtdz=xy-bzdt(5.1)其中s,r,b为变化区域有一定限制的实参数。该方程形式简单，表面上看并无惊人之处，但由该方程提示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。实验内容：先取定初值y0=(0,0,0),参数s=10,r=28,b=8/3,用MATLAB的数值求常微分方程函数ods45编程对（5.1）进行求解六

37、非线性方程（组）的求解1 求解方程（组）的solve命令求方程f(x)=q(x)的根可以用MATLAB命令:>> x=solve('方程f(x)=q(x)','待求符号变量x')求方程组fi(x1,xn)=qi(x1,xn) (i=1,2,n)的根可以用MATLAB命令:>>E1=sym('方程f1(x1,xn)=q1(x1,xn)');.En=sym('方程fn(x1,xn)=qn(x1,xn)');x1,x2,xn=solve(E1,E2,En, x1,xn)求解方程（组）的fsolve命令fsolve

38、的调用格式： X=fsolve(F,X0)2 逐步搜索法求解方程的根及其MATLAB程序逐步搜索法的MATLAB主程序function k,r=zhubuss(a,b,h,tol)% 输入的量- a和b是闭区间a,b的左、右端点；%-h是步长；%-tol是预先给定的精度.% 运行后输出的量-k是搜索点的个数；% - r是方程 在a,b上的实根的近似值，其精度是tol；X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n);for k=2:nX(k)=a+k*h;Y(k)=funs(X(k); %程序中调用的funs.m为函数sk=Y(

39、k)*Y(k-1);if sk<=0,m=m+1;r(m)=X(k);endxielv=(Y(k+1)-Y(k)*(Y(k)-Y(k-1);if (abs(Y(k)<tol)&( xielv<=0)m=m+1;r(m)=X(k);endend例1 用逐步搜索法的MATLAB程序分别求方程2x3+2x2-3x-3=0和3sin(cos2x)=0在区间-2,2上的根的近似值，要求精度是0.000 1.解 将逐步搜索法的MATLAB程序保存名为zhubuss.m的M文件.建立M文件funs.mfunction y=funs(x)y=2.*x.3+2.*x.2-3.*x-3在

40、MATLAB工作窗口输入如下程序>> k,r=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k =4001r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250 即搜索点的个数为k =4 001，其中有5个是方程2x3+2x2-3x-3=0的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度为0.000 1.在程序中将y=2.*x.3+2.*x.2-3.*x-3用y=sin(cos(2.*x.3) 代替，可得到方程sin(cos2x3)=0在区间-2,2上的根的近似值如

41、下r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.92301.3310 1.5780 1.7650 1.92003、二分法二分法的MATLAB主程序function k,x,wuca,yx=erfen(a,b,abtol)a(1)=a; b(1)=b;ya=fun(a(1); yb=fun(b(1); %程序中调用的fun.m 为函数 if ya* yb>0,disp('注意：ya*yb>0,请重新调整区间端点a和b.'), return endmax1=-1+ceil(log(b-a)- log(abtol)/ log(2); %ceil是向+ 方向取整for k=1: max1+1a;ya=fun(a); b;yb=fun(b); x=(a+b)/2;yx=fun(x); wuca=abs(b-a)/2; k=k-1;k,a,b,x,wuca,ya,yb