上机作业1.docx

1、上机作业1 2015级研究生计算方法作业 2015年11月上机作业1数值试验33-1试验目的：考察不动点迭代法的局部收敛性试验内容： 至少采用3种迭代法，迭代100次，考察收敛性，改变初值符号，再做迭代。分析收敛与发散的原因。(1)迭代原理：若实数满足，称为函数的一个不动点，迭代称为不动点迭代，称为迭代函数。由不动点方程建立迭代法，其中称为初值，需要预先给定。 方程分别对应下列不同形式的不动点方程： 取，按迭代，并分析收敛性。(2)不动点迭代法代码编制函数文件Iteratepro.mfunction y=iteratepro(x)x1=g(x);n=1;while (norm(x1-x) =1

2、.0e-5) &(n=1.0e-5）&（n=1.0e-5 x0=x1; x1=x0-fc(x0)/df(x0);n=n+1;end不动点迭代法编制函数文件Iteratepro.mfunction y=iteratepro(x)x1=g(x);n=1;while norm(x1-x)=1.0e-6 x=x1; x1=g(x);n=n+1;end函数文件g.mfunction y=g(x)y=sin(x);运行结果近似解初值x=1迭代次数n牛顿法1.9449e-0421不动点法0.0182k=410(3)为了提高牛顿法求重根的收敛阶采用以下俩种方法：方法一：程序：functionp1,err,y=

3、newtonroot(f,df,p0,eps,max1)p0for k=1:max1 p1=p0-2*（feval(f,p0)/feval(df,p0)）； err=abs(p1-p0); p0=p1; p1,err,k,y=feval(f,p1) if(erreps)|(y=0) break;end取初值p=0.5,当迭代次数k=9满足精度要求,近似解为2.4925e-05方法二：程序：functionp1,err,y=newtonroot(f,df,ddf,p0,eps,max1)p0for k=1:max1p1=p0-(feval(f,p0)*feval(df,p0)/(feval(df

4、,p0)2-feval(f,p0)*feval(ddf,p0)err=abs(p1-p0); p0=p1; p1,err,k,y=feval(f,p1) if(erreps)|(y=0) break; endend取初值p=0.5,k =3执行结果为：p1 =1.5022e-08 err =2.2741e-08 ans = 1.5022e-08以上四种方法中简单迭代法的收敛速度最慢 最后一种收敛速度最快。3-4 试验目的：体验Steffensens method加速技巧试验内容： 先用Newton法求解方程，再用Steffensens method求解，比较迭代步数，精确到0.00001。（1

5、）简单原理：其基本思想是：对给定的初值p0(0),首先应用不动点迭代法pn+1=g(pn)计算俩步，然后用Atiken公式加速。（2）程序代码Newton法迭代程序：functionp1,err,y=newtonroot(f,df,p0,eps,max1)p0for k=1:max1 p1=p0-feval(f,p0)/feval(df,p0); err=abs(p1-p0); p0=p1; p1,err,k,y=feval(f,p1) if(erreps)|(y=0) break; endEnd定义迭代函数：存储为f.mfunction y=f(x)y=x-tan(x)end定义迭代函数：存

6、储为df.mfunction g=df(x)g=1-1/(cos(x)2endAtiken公式加速程序：function n,err,p1,p2=steffen(g,p0,tol,max1)p0 n=0,p(1)=p0;while n=max1 for k=2:3 p(k)=feval(g,p(k-1); endp1=p(1)-(p(2)-p(1)2/(p(3)-2*p(2)+p(1);err=abs(p1-p(1);n=n+1;p(1)=p1n,err; if(err=1.0e-5 x0=y; y=B*x0+f;n=n+1;endynSeidel迭代法 D=diag(diag(A); D0=

7、inv(D); L=tril(A)-D; U=triu(A)-D; M=L+D; M0=inv(M); B=-M0*U; d=M0*b; x0=0 0 0 0 0; x=B*x0+d; while norm(x-x0)1e-5*norm(x)x0=x;x=B*x0+d;i=i+1;endSor迭代式求解方程组程序代码：建立文件名为sor.m的M文件，具体如下：function y=sor(a,b,w,x0)D=diag(diag(a);U=-triu(a,1);L=-tril(a,-1);B=(D-w*L)(1-w)*D+w*U);f=(D-w*L)b*w;y=B*x0+f;n=1;while

8、 norm(y-x0)=1.0e-5 x0=y; y=B*x0+f; n=n+1;Endyn2程序计算结果三种迭代法求解方程组（1）结果列表SorJacobiSeidel0.81.11.21.31.41.50.39060.39060.39060.39060.3906NaN0.39060.39060.16780.16780.16780.16780.1678NaN0.16780.16780.09960.09960.09960.09960.0996NaN0.09960.09960.17540.17540.17540.17540.1754NaN0.17540.17540.04470.04470.044

9、70.04470.0447-Inf0.04470.0447迭代次数10111323598357169三种迭代法求解方程组（2）结果列表SorJacobiSeidel0.81.11.21.31.41.50.47200.47200.47200.4720NaNNaN0.47200.47200.13300.13300.13300.1330NaNNaN0.13300.1330-0.1302-0.1302-0.1302-0.1302NaNNaN-0.1302-0.13020.16280.16280.16280.1628NaNNaN0.16280.16280.17010.17010.17010.1701Na

10、NNaN0.17010.1701迭代次数101727598058226024123结果分析：（1）Seidel迭代法收敛速度最快，迭代次数最少。（2）当松弛因子越大，其收敛速度越慢，迭代次数越多（3） 当松弛因子选的合适时（在1.0附近），Sor迭代法的收敛速度要比Jacobi迭代法快点。上机作业3数值试验5-2实验目的：掌握Newton法与最速下降法求解非线性方程组，观察各自的优势。实验内容：（1）分别用Newton法和最速下降法求解下面非线性方程组，初值精度要求 （2）改变初值，再用如上两种方法求解，得到什么结果？（3）采用初值，先用最速下降法求解3步，再用Newton迭代法求解，得到什么

11、结果？对以上运算结果作分析。(1)简单原理：Newton法:当方程组的解 可表示为称为求解非线性方程组的newton法迭代格式。最速下降法简单原理：在包含的某个区域内任取一初始点，然后在此点沿着下降最快的负梯度方向，取一个合适的步长，求得下一个，依次循环，直到来终止计算，此时就是满足要求的近似解。(2)MATLAB程序牛顿法解非线性方程组程序syms x1 x2 x3 f1=3*x1-cos(x2*x3)-(1/2);f2=(x1)2-81*(x2+0.1)2+sin(x3)+1.06;f3=exp(-x1*x2)+20*x3+(10*pi-3)/3);g1=diff(f1,x1);g2=di

12、ff(f1,x2);g3=diff(f1,x3);g4=diff(f2,x1);g5=diff(f2,x2);g6=diff(f2,x3);g7=diff(f3,x1);g8=diff(f3,x2);g9=diff(f3,x3);x0=0.1;0.1;-0.1;Tol=10(-5);% x0=20;20;20;F2=subs(f1,x1,x2,x3,x0(1),x0(2),x0(3); subs(f2,x1,x2,x3,x0(1),x0(2),x0(3); subs(f3,x1,x2,x3,x0(1),x0(2),x0(3);F3=subs(g1,x1,x2,x3,x0(1),x0(2),x0

13、(3) subs(g2,x1,x2,x3,x0(1),x0(2),x0(3) subs(g3,x1,x2,x3,x0(1),x0(2),x0(3); subs(g4,x1,x2,x3,x0(1),x0(2),x0(3) subs(g5,x1,x2,x3,x0(1),x0(2),x0(3) subs(g6,x1,x2,x3,x0(1),x0(2),x0(3); subs(g7,x1,x2,x3,x0(1),x0(2),x0(3) subs(g8,x1,x2,x3,x0(1),x0(2),x0(3) subs(g9,x1,x2,x3,x0(1),x0(2),x0(3);k=1;z=-inv(F3)

14、;y=z*F2;while k=20 x=x0+y; p=norm(x,inf);q=norm(x0,inf);if (p-q)/pTol) breakelse x0=x; F2=subs(f1,x1,x2,x3,x0(1),x0(2),x0(3); subs(f2,x1,x2,x3,x0(1),x0(2),x0(3); subs(f3,x1,x2,x3,x0(1),x0(2),x0(3); F3=subs(g1,x1,x2,x3,x0(1),x0(2),x0(3) subs(g2,x1,x2,x3,x0(1),x0(2),x0(3) subs(g3,x1,x2,x3,x0(1),x0(2),

15、x0(3); subs(g4,x1,x2,x3,x0(1),x0(2),x0(3) subs(g5,x1,x2,x3,x0(1),x0(2),x0(3) subs(g6,x1,x2,x3,x0(1),x0(2),x0(3); subs(g7,x1,x2,x3,x0(1),x0(2),x0(3) subs(g8,x1,x2,x3,x0(1),x0(2),x0(3) subs(g9,x1,x2,x3,x0(1),x0(2),x0(3); z=-inv(F3); y=z*F2; x=x0+y; k=k+1;endenddisp(x);disp(k);最速下降法解非线性方程程序syms x1 x2 x

16、3 tf1=3*x1-cos(x2*x3)-(1/2);f2=(x1)2-81*(x2+0.1)2+sin(x3)+1.06;f3=exp(-x1*x2)+20*x3+(10*pi-3)/3);F=(f1)2+(f2)2+(f3)2;g1=diff(F,x1);g2=diff(F,x2);g3=diff(F,x3);G11=g1;g2;g3;G1=-G11;x0=0.1;0.1;-0.1;G0=subs(G1,x1,x2,x3,x0(1),x0(2),x0(3);z0=x0+(t*G0);F1=subs(F,x1,x2,x3,z0(1),z0(2),z0(3);FD=diff(F1,t);xs

17、0=solve(FD,t);FD_xs0=vpa(subs(FD,t,xs0),6);F1_xs0=vpa(subs(F1,t,xs0),6);t=xs0;k=1;Tol=10(-5);while k=20 x=x0+(t*G0); p=norm(x,inf);q=norm(x0,inf);if (p-q)/pTol) breakelse x0=x; G0=subs(G1,x1,x2,x3,x0(1),x0(2),x0(3); x=x0+t*G0; k=k+1;endenddisp(x);disp(k);(3)程序运行的结果 （1）将初值定为分别运行上述两个程序，所得结果分别为： 牛顿法： 0

18、.5000 0.0000 -0.5236 4 最速下降法：0.12928443961888410941500498462202.0770* -0.52527958444225176826882692504685 3 （2）改变初值为，分别运行上述两个程序，所得结果分别为： 牛顿法： NaN NaN NaN 21 最速下降法： 20.219127745313790277417091321903 2.158*27343723864664313333 19.999648528574337249240101667477 2 （3）采用初值，先用最速下降法求解3步，再用牛顿迭代法求解。首先将第一个程序运

19、行3次后，将结果作为第二个程序的初值再做迭代，运行结果如下： 0.4274 -0.0318 -0.4736 2结果分析 Newton法的收敛速度较快，并且用相邻两次迭代的无穷范数相对误差来控制迭代次数，由程序运行的结果可知，改变初值后Newton法迭代的结果发散，由此可见Newton法对初值的要求较高。而最速下降法，对初值的要求不高，它对任意的初值都是收敛的，因此，为了提高迭代的效率，应先用最速下降法迭代几次，来确定初值，再用牛顿迭代法进行迭代。上机作业41.用规范的幂法与反幂法求矩阵A的按模最大，最小特征值与对应的特征向量。 A=,e=10-5（1）幂法算法步骤 (1)取初始向量（例如取）,

20、置精度要求，置。 (2)计算, (3)若，则停止计算（作为绝对值最大特征值，作为相应的特征向量）否则置，转(2)。程序代码function l,v,s=mifa(A,x0,eps)if nargin=2 eps=1.0e-5;endv=x0; M=5000; m=0;l=0;for k=1:M y=A*v; m=max(y); v=y/m; if(abs(m-l) A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2; x0=1 1 1 1; l,v,s=mifa(A,x0)结果如下：特征向量，特征值为 。（2）反幂法 计算的按模最小的特征值的问题就是计算的按模最大的特征值问

21、题。对于应用幂法迭代（又称反幂法），可求得矩阵的主特征值，从而求得的按模最小的特征值。反幂法的迭代公式为：任取初始向量，构造向量序列 迭代向量可以通过解线性方程组求得。程序代码functionl,v=fanmifa(A,x0,eps)if nargin=2 eps=1.0e-5;endv=x0; M=500; m=0;l=0;for k=1:M y=A*v;m=max(y); v=y/m;if(abs(m-l) A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2; x0=1 1 1 1; l,v=fanmifa(A,x0)结果如下：特征向量，特征值为 。（3）结论 幂法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。幂法的缺点是开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择