北京科技大学计算方法大作业概要

1、计算方法大作业机械电子工程系老师：廖福成注：本文本只有程序题，证明题全部在手写已交到理化楼204了。2. 证明方程 在上有一实根，并用二分法求这个根。要求。请给出程序和运行结果。证明：设f(x)=x3-x-1则f(1)= -1， f(2)= 5，f(1)*f(2)= -5<0因此，方程在1,2上必有一实根。二分法求解程序：%预先定义homework2.m文件如下：function lc=homework2(x)lc=x3-x-1;在MALAB窗口运行：cleara=1;b=2;tol=10(-3);N=10000;k=0; fa=homework2(a); % f 需事先定义for k=

2、1:Np=(a+b)/2; fp=homework2(p);if( fp=0 | (b-a)/2<tol)breakendif fa*fp<0 b=p; else a=p; endendk,p程序运行结果：k = 103. 用Newton迭代法求方程 的一个正根，计算结果精确到7位有效数字. 要求给出程序和运行结果. 解：取迭代初值 ,并设,则牛顿迭代函数为牛顿迭格式为：Matlab程序如下：%定义zuoye3.m文件function x=zuoye3(fname,dfname,x0,e,N)if nargin<5,N=500;endif nargin<4,e=1e-7

3、;endx=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<N, k=k+1; x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x)endif k=N,warning('已达上限次数');end在Matlab窗口中执行：zuoye3(inline('x3+2*x2+10*x-20'),inline('3*x2+4*x+10'),1,1e-7)结果如下：1.369336470588241.368808188617531.36880810782137ans

4、= 1.368808107821374. 用牛顿迭代法求方程在附近的根. 要求给出程序和运行结果. 解：令：,则牛顿迭代函数为牛顿迭格式为：Matalb程序如下：%定义zuoye4.m文件function x=zuoye4(fname,dfname,x0,e,N)if nargin<5,N=500;endif nargin<4,e=1e-7;endx=x0;x0=x+2*e;k=0;while abs(x0-x)>e&k<N, k=k+1; x0=x;x=x0-feval(fname,x0)/feval(dfname,x0); disp(x)endif k=N,

5、warning('已达上限次数');end在Matlab窗口执行：zuoye4(inline('x3-x-1'),inline('3*x2-1'),1,1e-7)结果如下：1.500000000000001.347826086956521.325200398950911.324718173999051.324717957244791.32471795724475ans = 1.324717957244756. 编写用全主元Gauss消去法解线性方程组的程序，并求解解：Matlab程序如下：A=2 -1 4 -3 1;-1 1 2 1 3;4 2

6、3 3 -1;-3 1 3 2 4;1 3 -1 4 4b=11 14 4 16 18function x=zuoye6(A,b)n,v=size(b);D=A,b;eye(n),zeros(n,v),s1,m=size(D);for k=1:(n-1) s=abs(A(k,k);p=k;q=k; for i=k:n for j=k:n if abs(A(i,j)>s s=abs(A(i,j);p=i;q=j; end end end if p>k t=D(k,:); D(k,:)=D(p,:); D(p,:)=t; end if q>k t1=D(:,k); D(:,k)=

7、D(:,q); D(:,q)=t1; end h=D(k+1:n,k)/D(k,k); D(k+1:n,k+1:m)=D(k+1:n,k+1:m)-h*D(k,k+1:m); D(k+1:n,k)=zeros(n-k,1);endfor k=n:-1:1 D(k,k:m)=D(k,k:m)/D(k,k); for r=1:k-1 D(r,:)=D(r,:)-D(r,k)*D(k,:); endendP=D(n+1:2*n,1:n);Q=D(1:n,n+1:m);x=P*Q在Matlab窗口中执行：A=0.02 -1 4 -3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4

8、;1 3 -1 4 4;b=11 14 4 16 18'zuoye6(A,b)运行结果如下：x = -3.82352941176472 1.000000000000007. 用追赶法解线性方程组 要求给出程序和运行结果. 解：于是有求解Ax=b即为求解，式中b=(1 0 0 0 0)T据 y=据= x=Matlab程序如下：%定义zuoye7.m文件function x=zuoye7(a,b,c,d)a1=0;a;n=length(b);q=zeros(n,1);p=zeros(n,1);%LU分解q(1)=b(1);for k=2:n,p(k)=a1(k)/q(k-1); q(k)=

9、b(k)-p(k)*c(k-1); end%解Ly=dy=zeros(n,1);y(1)=d(1);for k=2:n, y(k)=d(k)-p(k)*y(k-1);end%解Ux=yx=zeros(n,1); x(n)=y(n)/q(n);for k=n-1:-1:1,x(k)=(y(k)-c(k)*x(k+1)/q(k);endx在Matlab窗口中执行： a=-1 -1 -1 -1' b=2 2 2 2 2' c=-1 -1 -1 -1' d=1 0 0 0 0' x=zuoye7(a,b,c,d)运行结果如下：x = 0.83333333333333 0

10、.66666666666667 0.50000000000000 0.333333333333339. 分别用Jacobi迭代法和Gauss-Seidel迭代法求解程组（编写程序）取初值，精确到小数后面四位。解：（1） Jacobi迭代法的Matlab程序：% x0为初始向量，ep为精度，N为最大次数，x是近似解向量A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;n=length(b);N=500;ep=1e-6;x0=zero(n,1);n=length(b);x0=zeros(n,1);x=zeros(n,1);k=0while K<Nfor i=1:nx(i)

11、=(b(i)-A(I,1:i-1,i+1:n)*x0(1:i-1,i+1:n)/A(i,i);endif norm(x-x0,inf)<ep,break;end;x0=x;k=k+1;endif k=N,Warning(已达到迭代次数上限);enddisp(k=,num2str(k),x运行结果如下：k=15x= 1.000000327906423 0.999999801006905 1.000000327906432（2）Gauss-Seidel迭代法的Matlab程序：% x0为初始向量，ep为精度，N为最大次数，x是近似解向量Format long;clear;A=10 3 1;2

12、 -10 3;1 3 10;b=14 -5 14;n=length(b);N=500;ep=1e-6;x0=zero(n,1);P=inf;%以下是Guass-Seidal迭代法程序D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);dD=det(D);if dD=0disp(请注意：因为对角矩阵D奇异，所以此方程组无解。)elsedisp(请注意：因为对角矩阵D非奇异，所以此方程组无解。)iD=inv(D-L);B_GS=ID*U;d_GS=iD*b;endk=0;while k<Nx=B_GS*x0+d_GS;if norm(x-x0,P)<ep

13、,break;endx0=x;k=k+1;endif k=N,warning(已达到迭代次数上限)；enddisp(k=,num2str(k),x,%D,U,L,B_GS,d_GS,%jx=Ab,运行结果如下：请注意：因为对角矩阵D非奇异，所以此方程组有解。k=9x=1.000000043651052 1.000000031224555 0.99999998626752810编写幂法程序求矩阵 按模最大的特征值和对应的特征向量。解：幂法程序如下：% A 为n 阶方阵，u0 为初始向量，%epsilon 为上限，max1 为循环次数。lambda 返回按模最大的特征值，u 返回对应的特征向量。c

14、lear;format long;clc;max1=1000; epsilon=1e-6;A=1 1 0.5;1 1 0.25;0.5 0.25 2; u0=1;1;1; u=u0;v0=u0;lambda=0;k=0;err=1;R=;while (k<max1)&(err>epsilon)v=A*u; m,j=max(abs(v); m=v(j); u=v/m;err=abs(lambda-m);lambda=m; k=k+1;Temp=k;m;v;u;R=R Temp;endk, lambda, u, disp(R), xlswrite('d:1.xls

15、9;, R,'sheet1', 'b1') ;%D,Lmd=eig(A) %D 之列为A 的特征向量，Lmd 之对角线上为A 的特征值运行结果如下：k = 23lambda = 2.536527202038736u = 0.748222175139906 0.649662222855908 1.00000000000000011. 编写用原点位移加速反幂法程序求矩阵最接近于的特征值和相应的特征向量。 解：反幂法程序如下：% A 为n 阶方阵，v0,u0 为初始向量，%epsilon 为上限，max1 为循环次数，q 为近似特征值。lambda 为提高精度的特征值

16、，v 给出对应的特征向量。clear;format long;max1=100; epsilon=1e-10;A=7 3 -2;3 4 -1;-2 -1 3; u0=1;1;1; u=u0;v=u0;q=1.9;lambda=0;k=0;err=1;mu=0.5;n=length(u0); z=zeros(n,1);A=A-q*eye(n);%LU 分解n=length(u0);UU=zeros(n,n);L=eye(n,n);UU(1,:)=A(1,:); L(2:n,1)=A(2:n,1)/UU(1,1);for s=2:nUU(s,s:n)=A(s,s:n)-L(s,1:s-1)*UU(

17、1:s-1,s:n);L(s+1:n,s)=(A(s+1:n,s)-L(s+1:n,1:s-1)*UU(1:s-1,s)/UU(s,s);endwhile (k<max1)&(err>epsilon)m,j=max(abs(v); m=v(j); u=v/m;%解下三角方程组Lz=uz(1)=u(1);for j=2:n, z(j)=u(j)-L(j,1:j-1)*z(1:j-1); end%解上三角方程组Uv=zv(n)=z(n)/UU(n,n);for j=n-1:-1:1, v(j)=(z(j)-UU(j,j+1:n)*v(j+1:n)/UU(j,j); enderr

18、=abs(1/m-1/mu); k=k+1; mu=m;enddisp('迭代次数k= ',num2str(k)disp('要求的矩阵A 的特征值lambda='), lambda= q+1/mdisp('要求的矩阵A 的特征向量u'), u运行结果如下：迭代次数k= 17要求的矩阵A 的特征值lambda=lambda = 2.000000000009091要求的矩阵A 的特征向量uu = 0.999999999973922 -0.999999999956935 1.00000000000000012. 已知插值点(-2.00,17.00),

19、(0.00,1.00), (1.00,2.00), (2.00,17.00),求三次插值多项式，并计算. 解：其Matlab程序如下：function yy=zuoye10(x,y,xx) m=length(x);n=length(y); if m=n, error('向量x与y的长度必须一致');end s=0; for i=1:n t=ones(1,length(xx); for j=1:n if j=i t=t.*(xx-x(j)/(x(i)-x(j); end end s=s+t*y(i); end yy=s;在Matlab窗口中执行：x=-2.00 0.00 1.00

20、 2.00;y=17.00 1.00 2.00 17.00;xx=0.6;zuoye10(x,y,xx)结果如下：ans = 0.2560000000000013编制Newton插值法的通用Matlab程序，并求的近似值. 已知的数值如下表所示解：Newton插值法的通用Matlab程序：x=0 0.2 0.4 0.6 0.8'y=0.1995 0.3965 0.5881 0.7721 0.9461'xx=0.5;m=length(x);n=length(y);if m=n, error('向量x与y的长度必须一致');endfor j=2:n，for i=j:

21、nnewton_chazhi(i,j)=(newton_chazhi(i-1,j-1)-newton_chazhi(i,j-1)/(x(1+i-j)-x(i);endendyy= newton_chazhi(1,1);t=1;for i=2:n，t=t*(xx-x(i-1);yy=yy+ newton_chazhi(i,i)*t;end%计算误差近似值err=newton_chazhi(n,n);for i=1:n,err=err*(xx-x(i);end,erryy运行结果为：err =-2.343750000006213e-0620编写解常微分方程的四阶龙格库塔法通用程序. 解：四阶龙格库

22、塔法通用程序：function x,y=zuoye20(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x'y=y'21. 利用上题的程序求解初值问

23、题:的数值解,并将结果与准确解进行比较. 取.解：Matlab程序：function x,y=zuoye21(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2); k4=feval(dyfun,x(n+1),y(n)+h*k3); y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x&#

24、39;y=y'在Matlab窗口中执行：dyfun=inline('2*x*y(-2)/3');xspan=0,1.2;y0=1;h=0.4;nark18(dyfun,xspan,y0,h)结果如下：x = 0 0.40000000000000 0.80000000000000 1.20000000000000y =1.00000000000000 1.05075062243354 1.34631542500313计算值与准确值比较如下： 22. 编写四阶龙格库塔法程序，并求解洛伦兹型系统（不许调用现成的龙格库塔法程序软件包）：解：为了便于求解，这里首先为各个参数赋予初值，实际中根据需要修改程序中的参数即可。取，.取初始点为.程序如下：%预定义homework22.m文件function f=homework22(t,y)sgm=0.25;tao=0.06;xgm=0.