数值分析上机作业最强版

东南大学《数值分析》上机题数值分析上机题1(1) 编制按从大到小的顺序几=亠+42- -1 3~ — 1计算几的通用程序。

(2 )编制按从小到大由-走+詔E +H 计算“的通用程月(3) 按两种顺序分别计算％, %, %, 有效位数。

(编制程序时用单精度)(4) 通过本上机题，你明白了什么?程序代码(matlab 编程):clc cleara=single(1・/([2:10A 7]・ A 2-l)); Si (1)=single(0); SI (2)=1/(2A 2-1); for N=3:10A 2Sl(N)=a(l)； for i=2:N-lSI (N)=S1 (N)+a(i);endendS2 (l)=single(0); S2 (2)=1/(2A 2-1); for N=3:10A 2S2(N)=a(N-l);for i=linspace(N-2,1,N-2)S2(N)=S2(N)+a(i);endend其精确值为俣怙卜N —l顺序并指出S1表示按从大到小的顺序的S NS2表示按从小到大的顺序的S N 计算结果通过本上机题，看出按两种不同的顺序计算的结果是不相同的，按从大到小的顺序计算的值与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。

从大到小的顺序计算得到的结果的有效位数少。

计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，我们在计算机中进行同号数的加法时，采用绝对值较小者先加的算法，其结果的相对误差较小。

数值分析上机题220・（上机题）Newton迭代法（1）给定初值、及容许误差，，编制Newton法解方程/⑴“根的通用程序。

（2）给定方程弘—,易知其有三个根1.由Newton方法的局部收敛性可知存在5>o,当“（—恥）时,Newton迭代序列收敛于根工;。

试确定尽可能大的恥2•试取若干初始值,观察当x0 e （-00,-1）9（一1,一»）, （一恥），°1）,（is时Niwton序列是否收敛以及收敛于哪一个根。

数值分析上机实验一、解线性方程组直接法（教材49页14题）追赶法程序如下：function x=followup(A,b)n = rank(A);for(i=1:n)if(A(i,i)==0)disp('Error: 对角有元素为0');return;endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n)d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1) = (b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end主程序如下：function zhunganfaA=[2 -2 0 0 0 0 0 0;-2 5 -2 0 0 0 0 0;0 -2 5 -2 0 0 0 0;0 0 -2 5 -2 0 0 0;0 0 0 -2 5 -2 0 0;0 0 0 0 -2 5 -2 0;0 0 0 0 0 -2 5 -2;0 0 0 0 0 0 -2 5];b=[220/27;0;0;0;0;0;0;0];x=followup(A,b)计算结果：x =8.14784.07372.03651.01750.50730.25060.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor i=1:1:nA(i,j)=(1+0.1*i)^(j-1);endendc=cond(A)d=rcond(A)当n=5时c =5.3615e+005d =9.4327e-007当n=10时c =8.6823e+011d =5.0894e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。

数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。

但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。

对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。

求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。

使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。

求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。

实用标准文案文档大全上机作业题报告2015.1.9 USER1.Chapter 11.1题目设S N =∑1j 2−1N j=2，其精确值为)11123(21+--N N 。

（1）编制按从大到小的顺序11131121222-+⋯⋯+-+-=N S N ，计算S N 的通用程序。

（2）编制按从小到大的顺序1211)1(111222-+⋯⋯+--+-=N N S N ，计算S N 的通用程序。

（3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。

（编制程序时用单精度） （4）通过本次上机题，你明白了什么？1.2程序1.3运行结果1.4结果分析按从大到小的顺序，有效位数分别为：6，4，3。

按从小到大的顺序，有效位数分别为：5，6，6。

可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。

当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。

因此，采取从小到大的顺序累加得到的结果更加精确。

2.Chapter 22.1题目（1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。

（2）给定方程03)(3=-=x xx f ,易知其有三个根3,0,3321=*=*-=*x x x○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。

试确定尽可能大的δ。

○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？2.2程序2.3运行结果（1）寻找最大的δ值。

算法为：将初值x0在从0开始不断累加搜索精度eps，带入Newton迭代公式，直到求得的根不再收敛于0为止，此时的x0值即为最大的sigma值。

运行Find.m，得到在不同的搜索精度下的最大sigma值。

东南⼤学数值分析上机题答案数值分析上机题第⼀章17.（上机题）舍⼊误差与有效数设∑=-=Nj N j S 2211，其精确值为）111-23（21+-N N 。

（1）编制按从⼤到⼩的顺序1-1···1-311-21222N S N +++=，计算N S 的通⽤程序；（2）编制按从⼩到⼤的顺序121···1)1(111222-++--+-=N N S N ，计算NS 的通⽤程序；（3）按两种顺序分别计算210S ,410S ,610S ,并指出有效位数（编制程序时⽤单精度）；（4）通过本上机题，你明⽩了什么？解：程序：（1）从⼤到⼩的顺序计算1-1···1-311-21222N S N +++=：function sn1=fromlarge(n) %从⼤到⼩计算sn1format long ; sn1=single(0); for m=2:1:nsn1=sn1+1/(m^2-1); end end（2）从⼩到⼤计算121···1)1(111222-++--+-=N N S N function sn2=fromsmall(n) %从⼩到⼤计算sn2format long ; sn2=single(0); for m=n:-1:2sn2=sn2+1/(m^2-1); end end（3）总的编程程序为： function p203()clear allformat long;n=input('please enter a number as the n:') sn=1/2*(3/2-1/n-1/(n+1));%精确值为sn fprintf('精确值为%f\n',sn);sn1=fromlarge(n);fprintf('从⼤到⼩计算的值为%f\n',sn1);sn2=fromsmall(n);fprintf('从⼩到⼤计算的值为%f\n',sn2);function sn1=fromlarge(n) %从⼤到⼩计算sn1 format long;sn1=single(0);for m=2:1:nsn1=sn1+1/(m^2-1);endendfunction sn2=fromsmall(n) %从⼩到⼤计算sn2 format long;sn2=single(0);for m=n:-1:2sn2=sn2+1/(m^2-1);endendend运⾏结果：从⽽可以得到N值真值顺序值有效位数2 100.740050 从⼤到⼩0.740049 5从⼩到⼤0.740050 64 100.749900 从⼤到⼩0.749852 3从⼩到⼤0.749900 66 100.749999 从⼤到⼩0.749852 3从⼩到⼤0.749999 6（4）感想：通过本上机题，我明⽩了，从⼩到⼤计算数值的精确位数⽐较⾼⽽且与真值较为接近，⽽从⼤到⼩计算数值的精确位数⽐较低。

第一章第二题(1) 截断误差为104-时：k=1;n=0;m=0;x=0;e=1e-4;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141792613595791m =5001（2）截断误差为108-时：k=1;n=0;m=0;x=0;e=1e-8;while k==1x1=x+(-1)^n/(2*n+1);if abs(x-x1)<ey=4*x1;m=n+1;break;endx=x1;k=1;n=n+1;endformat longy,my =3.141592673590250m =50000001由以上计算可知，截断误差小于104-时，应取5001项求和，π=3.141792613595791；截断误差小于108-时，应取50000001项求和，π=3.141592673590250。

第二章第二题a=[0 -2 -2 -2 -2 -2 -2 -2];b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0];v=220;r=27;d=[v/r 0 0 0 0 0 0 0];n=8;for i=2:na(i)=a(i)/b(i-1);b(i)=b(i)-c(n-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(n)=d(n)/b(n);for i=n-1:-1:1d(i)=(d(i)-c(i)*d(i+1));end;I=d'I =1.0e+002 *1.490717294184090.704617906351300.311568212434910.128623612390290.049496991380330.017168822994210.004772412363470.00047741598598第三章第一题(1)Jacobi迭代法：b=[12;-27;14;-17;12]x = [0;0;0;0;0;]k = 0;r = 1;e=0.000001A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;] D = diag(diag(A));B = inv(D)*(D-A);f = inv(D)*b;p = max(abs(eig(B)));if p >= 1'迭代法不收敛'returnendwhile r >ex0 = x;x = B*x0 + f;k = k + 1;r = norm (x-x0,inf);endxk计算结果：x =1.0000-2.00003.0000-2.00001.0000k =65(2) 高斯赛德尔迭代：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;]x=[0;0;0;0;0];b=[12;-27;14;-17;12]c=0.000001L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-L)*U*x+inv(D-L)*b;k=1;while norm(X-x,inf)>= cx=X;X=inv(D-L)*U*x+inv(D-L)*b;k=k+1;endXk计算结果：X =1.0000-2.00003.0000-2.00001.0000k =37(3) SOR：A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15] x=[0;0;0;0;0];b=[12;-27;14;-17;12]e=0.000001w=1.44;L=-tril(A,-1)U=-triu(A,1)D=(diag(diag(A)))X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*bn=1;while norm(X-x,inf)>=ex=X;X=inv(D-w*L)*((1-w)*D+w*U)*x+w*inv(D-w*L)*b;n=n+1;endXn计算结果：X =1.0000-2.00003.0000-2.00001.0000n =22由以上可知，共轭梯度法收敛速度明显快于Jacobi法和G-S法。

数值分析上机题作业电器学院交通信息工程及控制罗宁20050290010110第一章源程序：REAL SS=1DO I=1,71S=S*29/IEND DOWRITE(*,*) "Y=",SEND结果：Y=79.5416第三章源程序：REAL X,Y,N,K,DX,NEWTON,NE,S,XKDIMENSION X(6),Y(6),NE(6,6),NEWTON(21)DATA X/0.2,0.24,0.28,0.32,0.36,0.4/DATA Y/0.1987,0.2377,0.2764,0.3146,0.3523,0.3894/!计算差商表：DX=0.04DO I=1,6NE(I,1)=Y(I)END DODO J=2,6DO I=J,6NE(I,J)=(NE(I,J-1)-NE(I-1,J-1))/(DX*(J-1))END DOEND DODO I=1,21XK=0.2+(I-1)*0.01NEWTON(I)=NE(1,1)S=1DO J=2,6S=S*(XK-X(J-1))NEWTON(I)=NEWTON(I)+S*NE(J,J)END DOEND DOWRITE(*,*) "NEWTON=",NEWTONEND结果：NEWTON=0.198700 0.208451 0.218209 0.227962 0.2377000.2474170.257108 0.266770 0.276400 0.285998 0.2955640.305097 0.3146000.324072 0.333513 0.342923 0.352300 0.3616420.370943 0.3801990.389400第四章源程序：REAL A,B,ANS,DDIMENSION A(4,4),B(4),ANS(4),D(4,5)DATA A/4197,6.8,88.6,1.45,305,71.3,76.4,5.9,-206,&&-47.4,-10.8,6.13,-840,52,802,36.5/DATA B/136,11.5,25.7,6.6/CALL LGAUS(A,B,4,D)CALL GAUSQ(D,4,ANS)WRITE(*,*) 'ANS=',ANSEND结果：ANS=3.226190E-02 0.322386 0.291290 -1.684955E-02第五章源程序：REAL A1,B1,A2,B2,E1,E2A1=0;A2=0B1=3.141592654/2;B2=1E1=0.000001;E2=0.001CALL SIMPB1(A1,B1,E1,S)WRITE(*,*) '第一题simpson:',SCALL SIMPB2(A2,B2,E1,S)WRITE(*,*) '第二题simpson:',SCALL ROMBERG1(A1,B1,E2,R)WRITE(*,*) '第一题ROMBERG',RCALL ROMBERG2(A2,B2,E2,R)WRITE(*,*) '第二题ROMBERG',RENDSUBROUTINE SIMPB1(A,B,E,S)F(X)=SIN(X)/XH=(B-A)/2S2=0;N=1S0=1+F(B)S1=F(A+H)S=(S0+4.0*S1)*H/3.060 N=2*N;H=H/2.0S2=S2+S1S1=0.0X=A+HDO I=1,NS1=S1+F(X)X=X+H+HEND DOS2N=(S0+2.0*S2+4.0*S1)*H/3.0IF(ABS(S2N-S)>E)THENS=S2NGOTO 60END IFEND SUBROUTINESUBROUTINE SIMPB2(A,B,E,S)F(X)=LOG(1+X)/XH=(B-A)/2S2=0;N=1S0=1+F(B)S1=F(A+H)S=(S0+4.0*S1)*H/3.060 N=2*N;H=H/2.0S2=S2+S1S1=0.0X=A+HDO I=1,NS1=S1+F(X)X=X+H+HEND DOS2N=(S0+2.0*S2+4.0*S1)*H/3.0IF(ABS(S2N-S)>E)THENS=S2NGOTO 60END IFEND SUBROUTINESUBROUTINE ROMBERG1(A,B,E,R)F(X)=SIN(X)/XREAL S,DS,TINTEGER K,UDIMENSION S(100,100)S(1,1)=(B-A)*(F(B)+1)/2K=1T=0100 K=K+1S(K,1)= 1/2*S((K-1),1)U=2**(K-2)-1DO J=0,UT=T+F(A+(2*J+1)*(B-A)/(2**(K-1)))END DOS(K,1)=S(K,1)+T*(B-A)/2**(K-1)DO J=2,KS((K-J+1),J)=(4**(J-1)*S((K-J+2),(J-1))+S((K-J+1),(J-1)))/&&(4**(J-1)-1)END DODS=S(1,K)-S(1,K-1)IF(ABS(DS).GT.E)THENR=S(1,K)ELSEGOTO 100END IFEND SUBROUTINESUBROUTINE ROMBERG2(A,B,E,R)F(X)=LOG(X+1)/XREAL S,DS,TINTEGER K,UDIMENSION S(100,100)S(1,1)=(B-A)*(F(B)+1)/2K=1T=0100 K=K+1S(K,1)= 1/2*S((K-1),1)U=2**(K-2)-1DO J=0,UT=T+F(A+(2*J+1)*(B-A)/(2**(K-1)))END DOS(K,1)=S(K,1)+T*(B-A)/2**(K-1)DO J=2,KS((K-J+1),J)=(4**(J-1)*S((K-J+2),(J-1))+S((K-J+1),(J-1)))/&&(4**(J-1)-1)END DODS=S(1,K)-S(1,K-1)IF(ABS(DS).GT.E)THENR=S(1,K)ELSEGOTO 100END IFEND SUBROUTINE结果：第一题simpson: 1.37076第二题simpson: 0.822467第一题ROMBERG 1.37128第二题ROMBERG 0.822811第六章第一题源程序：REAL X0,E,X1,X2,XXXX=1.55X0=1.5E=0.00005CALL SUPPO(X0,E,X1)WRITE(*,*) "X1=",X1CALL SUPPO1(X0,E,X2)WRITE(*,*) "X2=",X2CALL NTSP1(X0,XX,E,X4)WRITE(*,*) "X4=",X4CALL NTQX(X0,E,X3)WRITE(*,*) "X3=",X3END!简单迭代SUBROUTINE SUPPO(X0,E,X1)G(X)=1+1/x**2I=1X1=X0100 X2=X1X1=G(X1)I=I+1Y=X1-X2IF(ABS(Y)>E) THENGOTO 100ENDIFEND SUBROUTINE SUPPO!加速迭代SUBROUTINE SUPPO1(X0,E,X2) G(X)=1+1/X**2G1(X)=-2/X**3I=1X=X0200 X2=G(X)Q=G1(X2)X2=(X2+Q*X)/(1-Q)I=I+1Y=X2**3-X2**2-1X=X2IF(ABS(Y)>E) THENGOTO 200ENDIFEND SUBROUTINE SUPPO1!NEWTONSUBROUTINE NTQX(X0,E,X2)G(X)=X**3-X**2-1G1(X)=3*X**2-2*XX1=X0300 F=G(X1)F1=G1(X1)X2=X1-F/F1X1=X2I=I+1U=X2**3-X2**2-1IF(ABS(U)>E)THENGOTO 300END IFEND SUBROUTINE NTQXSUBROUTINE NTSP1(X0,XX,E,X4)G(X)=X**3-X**2-1I=1X1=XXF1=G(X0)400 F=G(X1)X4=X1-F*(X1-X0)/(F-F1)I=I+1Y=X4**3-X4**2-1X1=X4IF(ABS(Y)>E) THENGOTO 400END IFEND SUBROUTINE NTSP1结果：X1= 1.46556X2= 1.46558X3= 1.46667X4= 1.46557第二题源程序：REAL A,B,E,XDIMENSION A(3,3),B(3),X(3)DATA A/-8,1,1,1,-5,1,1,1,-4/DATA B/1,16,7/E=0.0001CALL YACOBI(A,B,3,E,X)WRITE(*,*) "方程的解为：",XCALL GAUSEI(A,B,3,E,X)WRITE(*,*) "方程的解为：",XENDSUBROUTINE YACOBI(A,B,N,E,X)REAL A,B,X,YDIMENSION A(N,N),B(N),X(N),Y(N)DO I=1,NX(I)=0END DO60 D=0DO I=1,NY(I)=B(I)DO J=1,NIF(I.NE.J)THENY(I)=Y(I)-A(I,J)*X(J)END IFEND DOY(I)=Y(I)/A(I,I)IF(ABS(X(I)-Y(I))>D)THEND=ABS(X(I)-Y(I))ENDIFEND DODO I=1,NX(I)=Y(I)ENDDOIF(D>E)THENGO TO 60END IFEND SUBROUTINE YACOBISUBROUTINE GAUSEI(A,B,N,E,X)DIMENSION A(N,N),B(N),X(N)DO I=1,NX(I)=0END DO100 D=0DO I=1,NY=B(I)DO J=1,NIF(I.NE.J) THENY=Y-A(I,J)*X(J)ENDIFEND DOY=Y/A(I,I)IF(ABS(X(I)-Y)>D)THEND=ABS(X(I)-Y)END IFX(I)=YEND DOIF(D>E)THENGOTO 100END IFEND SUBROUTINE GAUSEI结果：方程的解为：-0.999985 -3.99998 -2.99998 方程的解为：-0.999988 -3.99999 -2.99999第七章第一题源程序：REAL A,B,Y0,YY,Y,XINTEGER N1,N2,N3A=0.0;B=1.0Y0=1.0N1=10;N2=100;N3=1000CALL RKUTTA(A,B,N1,Y0,YY)WRITE(*,*) '步长为0.1，y(1)=',YYCALL RKUTTA(A,B,N2,Y0,YY)WRITE(*,*) '步长为0.01，y(1)=',YYCALL RKUTTA(A,B,N3,Y0,YY)WRITE(*,*) '步长为0.001，y(1)=',YYENDSUBROUTINE RKUTTA(A,B,M,Y0,YY)F(X,Y)=EXP(X)*Y**2-2*YREAL X0,A,B,Y0,YY,HHDIMENSION HH(5)X0=A;Y=Y0H=(B-A)/MHH(1)=H/2;HH(2)=HH(1);HH(3)=H;HH(4)=H;HH(5)=H/2DO I=1,MX1=X0;Y1=Y;YY=YDO J=1,4X=X1+HH(J)AA=F(X,Y)YY=YY+AA*HH(J+1)/3Y=Y1+AA*HH(J)END DOY=YYEND DOEND SUBROUTINE RKUTTA结果：步长为0.1，y(1)= 0.253325步长为0.01，y(1)= 0.239780步长为0.001，y(1)= 0.238542第二题源程序：REAL A,B,Y,YYINTEGER N,MDIMENSION Y(3),YY(3)DATA Y/1,-1,0/N=5/0.25A=0.0;B=5.0M=3CALL HRKUTTA(A,B,N,M,Y,YY)WRITE(*,*) '方程的解为：',YYENDSUBROUTINE HRKUTTA(A,B,N,M,Y,YY)REAL X0DIMENSION HH(5),Y(M),Y1(M),F(M),YY(M)X0=AH=(B-A)/NHH(1)=H/2;HH(2)=HH(1);HH(3)=H;HH(4)=H;HH(5)=H/2DO I=1,NX1=X0+(I-1)*HX=X1DO J=1,MY1(J)=Y(J);YY(J)=Y(J)END DODO K=1,4CALL FUN(X,Y,F)X=X1+HH(K)DO L=1,MYY(L)=YY(L)+F(L)*HH(K+1)/3Y(L)=Y1(L)+F(L)*HH(K)END DOEND DODO J=1,MY(J)=YY(J)END DOEND DOEND SUBROUTINE HRKUTTASUBROUTINE FUN(X,Y,F)REAL X,Y,FDIMENSION Y(3),F(3)F(1)=-Y(1)+2*Y(2)+6*Y(3)F(2)=-1*Y(2)+3*Y(3)+2*SIN(X)F(3)=-1*Y(3)+X**2*EXP(-1*X)+COS(X)END SUBROUTINE FUN结果：方程的解为：-1.42051 -1.67879 -6.023698E-02。

数值分析上机作业2实验一：（1） ①用不动点迭代法求()013=--=x x x f 的根，至少设计两种迭代格式，一个收敛一个发散，1210-=ε。

（2） ②对迭代格式使用Aitken 加速，观察敛散性变化。

1取递推公式31)1(+=x x ，可以得到收敛时的迭代结果为：x=(2)^(1/3); t=1; while(1)if(abs(x-(x+1)^(1/3))<10^-12) break; endx=(x+1)^(1/3);t=t+1;end t xtemp=x^3-x-1 %带回来验证下 t = 16 x =1.324717957243755 temp =-4.225952920933196e-12 加速后代码如下 x=1;x=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)) t=1; while(1)if(abs(x-(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)))<10^-10) break; endx=(x*((x+1)^(1/3)+1)^(1/3)-(x+1)^(2/3))/(x-2*(x+1)^(1/3)+((x+1)^(1/3)+1)^(1/3)); t=t+1; if (t>100000)break; %防止运算次数过多 end endfprintf('需要%d 次',t);输出需要670次>>此处取10^10-=ε，若到10^-12次方则可能需要运行更多次取13-=x x 则迭代发散。

使用aitken 加速计算结果如下 x=2; t=1; while(1)if( abs(x-((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1))<10^-10) break; end t=t+1;x=((x*((x^3-1)^3-1))-(x^3-1)^2)/(x-2*(x^3-1)+(x^3-1)^3-1); if (t>100000)break; %防止运算次数过多 end end t xt = 108x = 1.324717956244172由此可见经过aitken 加速以后，原来发散的迭代格式收敛了。