北京科技大學(xué)計(jì)算方法大作業(yè)

1、計(jì)算方法大作業(yè)機(jī)械電子工程系老師：廖福成注：本文本只有程序題，證明題全部在手寫(xiě)已交到理化樓204了。2. 證明方程 在上有一實(shí)根，并用二分法求這個(gè)根。要求。請(qǐng)給出程序和運(yùn)行結(jié)果。證明：設(shè)f(x)=x3-x-1則f(1)= -1， f(2)= 5，f(1)*f(2)= -5<0因此，方程在1,2上必有一實(shí)根。二分法求解程序：%預(yù)先定義homework2.m文件如下：function lc=homework2(x)lc=x3-x-1;在MALAB窗口運(yùn)行：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程序運(yùn)行結(jié)果：k = 10p = 1.3251953125000003. 用Newton迭代法求方程 的一個(gè)正根，計(jì)算結(jié)果精確到7位有效數(shù)字. 要求給出程序和運(yùn)行結(jié)果. 解：取迭代初值 ,并設(shè),則牛頓迭代函數(shù)為牛頓迭格式為：Matlab程序如下：%定義zuoye3.m文件function x=zuoye3(fname,dfname,x0,e,N)if nargin<5,N=500;end

3、if 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,warning('已達(dá)上限次數(shù)');end在Matlab窗口中執(zhí)行：zuoye3(inline('x3+2*x2+10*x-20'),inline('3*x2+4*x+10'),1,1e-7)結(jié)果如下：1.411764705882351.3693364705882

4、41.368808188617531.36880810782137ans = 1.368808107821374. 用牛頓迭代法求方程在附近的根. 要求給出程序和運(yùn)行結(jié)果. 解：令：,則牛頓迭代函數(shù)為牛頓迭格式為：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

5、)/feval(dfname,x0); disp(x)endif k=N,warning('已達(dá)上限次數(shù)');end在Matlab窗口執(zhí)行：zuoye4(inline('x3-x-1'),inline('3*x2-1'),1,1e-7)結(jié)果如下：1.500000000000001.347826086956521.325200398950911.324718173999051.324717957244791.32471795724475ans = 1.324717957244756. 編寫(xiě)用全主元Gauss消去法解線性方程組的程序，并求解解：Mat

6、lab程序如下：A=2 -1 4 -3 1;-1 1 2 1 3;4 2 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,:

7、)=t; end if q>k t1=D(:,k); D(:,k)=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窗口中執(zhí)行：A=0.02 -1 4

8、-3 1;-1 1 2 1 3;4 2 3 3 -1;-3 1 3 2 4;1 3 -1 4 4;b=11 14 4 16 18'zuoye6(A,b)運(yùn)行結(jié)果如下：x = 2.94117647058824 -3.82352941176472 1.00000000000000 0.94117647058824 5.941176470588247. 用追趕法解線性方程組 要求給出程序和運(yùn)行結(jié)果. 解：于是有求解Ax=b即為求解，式中b=(1 0 0 0 0)T據(jù) y=據(jù)= x=Matlab程序如下：%定義zuoye7.m文件function x=zuoye7(a,b,c,d)a1=0;a

9、;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)=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窗口中執(zhí)行： a=-1 -1 -1 -1' b=2 2 2 2 2

10、' c=-1 -1 -1 -1' d=1 0 0 0 0' x=zuoye7(a,b,c,d)運(yùn)行結(jié)果如下：x = 0.83333333333333 0.66666666666667 0.50000000000000 0.33333333333333 0.166666666666679. 分別用Jacobi迭代法和Gauss-Seidel迭代法求解程組（編寫(xiě)程序）取初值，精確到小數(shù)后面四位。解：（1） Jacobi迭代法的Matlab程序：% x0為初始向量，ep為精度，N為最大次數(shù)，x是近似解向量A=10 3 1;2 -10 3;1 3 10;b=14 -5 14;n

11、=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)=(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(已達(dá)到迭代次數(shù)上限);enddisp(k=,num2str(k),x運(yùn)行結(jié)果如下：k=15x= 1.000000327906423 0.99999980100690

12、5 1.000000327906432（2）Gauss-Seidel迭代法的Matlab程序：% x0為初始向量，ep為精度，N為最大次數(shù)，x是近似解向量Format long;clear;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);P=inf;%以下是Guass-Seidal迭代法程序D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);dD=det(D);if dD=0disp(請(qǐng)注意：因?yàn)閷?duì)角矩陣D奇異，所以此方程組無(wú)解。)elsedisp(請(qǐng)注意：因?yàn)閷?duì)

13、角矩陣D非奇異，所以此方程組無(wú)解。)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,break;endx0=x;k=k+1;endif k=N,warning(已達(dá)到迭代次數(shù)上限)；enddisp(k=,num2str(k),x,%D,U,L,B_GS,d_GS,%jx=Ab,運(yùn)行結(jié)果如下：請(qǐng)注意：因?yàn)閷?duì)角矩陣D非奇異，所以此方程組有解。k=9x=1.000000043651052 1.000000031224555 0.99999998626752810編寫(xiě)冪法程序

14、求矩陣 按模最大的特征值和對(duì)應(yīng)的特征向量。解：冪法程序如下：% A 為n 階方陣，u0 為初始向量，%epsilon 為上限，max1 為循環(huán)次數(shù)。lambda 返回按模最大的特征值，u 返回對(duì)應(yīng)的特征向量。clear;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;e

15、rr=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', R,'sheet1', 'b1') ;%D,Lmd=eig(A) %D 之列為A 的特征向量，Lmd 之對(duì)角線上為A 的特征值運(yùn)行結(jié)果如下：k = 23lambda = 2.536527202038736u = 0.748222175139906 0.649662222855908 1.00000000000000011. 編寫(xiě)用原點(diǎn)位移加速反冪法

16、程序求矩陣最接近于的特征值和相應(yīng)的特征向量。 解：反冪法程序如下：% A 為n 階方陣，v0,u0 為初始向量，%epsilon 為上限，max1 為循環(huán)次數(shù)，q 為近似特征值。lambda 為提高精度的特征值，v 給出對(duì)應(yīng)的特征向量。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=z

17、eros(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(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

18、(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=abs(1/m-1/mu); k=k+1; mu=m;enddisp('迭代次數(shù)k= ',num2str(k)disp('要求的矩陣A 的特征值lambda='), lambda= q+1/mdisp('要求的矩陣A 的特征向量u'), u運(yùn)行結(jié)果如下：迭代次數(shù)k= 17要求的矩陣A 的特征值lambda=lambda = 2.00000

19、0000009091要求的矩陣A 的特征向量uu = 0.999999999973922 -0.999999999956935 1.00000000000000012. 已知插值點(diǎn)(-2.00,17.00), (0.00,1.00), (1.00,2.00), (2.00,17.00),求三次插值多項(xiàng)式，并計(jì)算. 解：其Matlab程序如下：function yy=zuoye10(x,y,xx) m=length(x);n=length(y); if m=n, error('向量x與y的長(zhǎng)度必須一致');end s=0; for i=1:n t=ones(1,length(xx

20、); 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窗口中執(zhí)行：x=-2.00 0.00 1.00 2.00;y=17.00 1.00 2.00 17.00;xx=0.6;zuoye10(x,y,xx)結(jié)果如下：ans = 0.2560000000000013編制Newton插值法的通用Matlab程序，并求的近似值. 已知的數(shù)值如下表所示解：Newton插值法的通用Matlab程序：x=0 0.2 0.4 0.6 0.8'y=0.1995 0.3965 0.5881 0.7

21、721 0.9461'xx=0.5;m=length(x);n=length(y);if m=n, error('向量x與y的長(zhǎng)度必須一致');endfor j=2:n，for i=j: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%計(jì)算誤差近似值err=newton_chazhi(

22、n,n);for i=1:n,err=err*(xx-x(i);end,erryy運(yùn)行結(jié)果為：err =-2.343750000006213e-06yy=0.68118750000000020編寫(xiě)解常微分方程的四階龍格庫(kù)塔法通用程序. 解：四階龍格庫(kù)塔法通用程序：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

23、(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. 利用上題的程序求解初值問(wèn)題:的數(shù)值解,并將結(jié)果與準(zhǔn)確解進(jìn)行比較. 取.解：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

24、)+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'在Matlab窗口中執(zhí)行：dyfun=inline('2*x*y(-2)/3');xspan=0,1.2;y0=1;h=0.4;nark18(dyfun,xspan,y0,h)結(jié)果如下：x = 0 0.40000000000000 0.80000000000000 1.20000000000000y =1.00000000000000 1.05075062243354 1.17933175912031 1.34631542500313計(jì)算值與準(zhǔn)確值比較如下： 22. 編寫(xiě)四階龍格庫(kù)塔法程序，并求解洛倫茲型系統(tǒng)（不許調(diào)用現(xiàn)成的龍格庫(kù)塔法程序軟件包）：解：為了便于求解，這里首先為各個(gè)參數(shù)賦予初值，實(shí)際中根據(jù)需要修改程序中的參數(shù)即可。取，.取初始點(diǎn)為.程序如下：%預(yù)定義homework22.m文件function f=homework22(t,y)sgm=0.25;tao