北科计算方法上机作业

1、计算方法上机作业姓 名： 班 级： 学 号： 单 位： 授课教师： 目 录上机作业11上机作业25上机作业39上机作业411上机作业514上机作业616上机作业11、分别用不动点迭代与Newton法求解方程3x-5x+4=0的正根与负根。解：a）不动点迭代：.迭代公式pk+1=log5(3pk+4) 初值p0=1程序以及结果：p0=1;N=2000;Tol=10(-4); k=1;while k=Np=log(3*p0+4)/log(5);if abs(p-p0)Tol break endk=k+1;p0=p;enddisp(p);disp(k) 结果：近似解以及迭代次数 迭代公式pk+1=

2、5pk-43 初值p0=-1程序以及结果：p0=-1;N=2000;Tol=10(-4);k=1;while k=Np=(5p0-4)/3;if abs(p-p0)Tol break endk=k+1;p0=p; enddisp(p);disp(k) 结果：近似解以及迭代次数b）Newton法迭代：迭代公式pk+1=pk-(3pk-5pk+4)/(3-5pk*ln5) 初值p0=1程序以及结果： p0=1;N=2000;Tol=10(-4);for k=1:Np1=p0-(3*p0-5p0+4)/(3-5p0*log(5);if abs(p1-p0)Tolbreakendp0=p1;enddi

3、sp(p1);disp(k) 结果：近似解以及迭代次数 当初值p0取-1时，迭代公式不变，仍是pk+1=pk-(3pk-5pk+4)/(3-5pk*ln5)，程序不变（只是p0=-1）结果如下：近似解以及迭代次数结果分析：当初值p0=1时，设不动点每步迭代的结果为p牛顿迭代法每步迭代的结果为p1，则n01234567p11.20911.26241.27531.27841.27911.27931.2793p111.39631.29221.27951.27931.2793当初值p0=-1时，设不动点每步迭代的结果为p，牛顿迭代法每步迭代的结果为p1，则n012345p-1-1.2667-1.289

4、9-1.2915-1.2916-1.2916p1-1-1.2987-1.2916-1.2916由以上的迭代结果可以看出，Newton法迭代相对于不动点迭代法，收敛速度更快。2、用Newton法与重根计算法求解方程 x-sinx0 的根.再用 Steffensens method加速Newton法收敛，比较结果。解：A) Newton法求解：迭代公式 pk+1=pk-pk-sin(pk)/(1-cospk)程序以及结果： p0=1;N=2000;Tol=10(-4);for k=1:Np1=p0-(p0-sin(p0)/(1-cos(p0);if abs(p1-p0)Tolbreakendp0=

5、p1;enddisp(p1);disp(k) 近似解以及迭代次数：B) 重根法求解：设f(x)= x-sinx，则有f(0)=f(0)=f(0)，f(3)(0)=1，则p=0是方程的三重根。于是有迭代公式 pk+1=pk-3pk-sin(pk)/(1-cospk)程序以及结果：p0=1;N=2000;Tol=10(-4);for k=1:Np1=p0-3*(p0-sin(p0)/(1-cos(p0);if abs(p1-p0)Tolbreakendp0=p1;enddisp(p1);disp(k) 近似解以及迭代次数： C) Steffensens method加速Newton法收敛：程序以及

6、结果：p0=1;p1=0.6551;N=2000;Tol=10(-4);n=0;p(1)=p0;while n=Np(2)=p(1)-(p(1)-sin(p(1)/(1-cos(p(1);p(3)=p(2)-(p(2)-sin(p(2)/(1-cos(p(2);p1=p(1)-(p(2)-p(1)2/(p(3)-2*p(2)+p(1);f0=p1-sin(p1);if abs(f0) A=-14*eye(10)+ones(10,10);b=-1*ones(10,1);x0=zeros(10,1);e=10(-5);N=2000;x=x0;for k=1:Nfor i=1:10x(i)=(b(i

7、)-A(i,:)*x0)/A(i,i)+x0(i);endif norm(x-x0)=Ndisp(发散)endend 结果：近似解以及迭代次数2、Seidel法：迭代格式 xi(k+1)=1aii(bi-),i=1,2,10,k=1,2,程序以及结果： A=-14*eye(10)+ones(10,10);b=-1*ones(10,1);x0=zeros(10,1);e=10(-5);N=2000;x=x0;for k=1:Nfor i=1:10x(i)=(b(i)-A(i,:)*x)/A(i,i)+x(i); endif norm(x-x0)=Ndisp(发散)endend 结果： 近似解以及

8、迭代次数3、Sor迭代法：迭代格式 xi(k+1)=(1-)xi(k)+aii(bi-), i=1,2,10,k=1,2,程序以及结果：（=0.8时） A=-14*eye(10)+ones(10,10);b=-1*ones(10,1);x0=zeros(10,1);e=10(-5);N=2000;x=x0;w=0.8;for k=1:Nfor i=1:10x(i)=w*(b(i)-A(i,:)*x)/A(i,i)+x(i);endif norm(x-x0)=Ndisp(发散)endend 结果：近似解以及迭代次数当分别取1.1,1.2,1.3,1.4,1.5时，程序不变，只是改变赋值条件中w的

9、值。当=1.1时， =1.2时， =1.3时， =1.4时， =1.5时 结果分析：Jacobi法的迭代列表k12345678910x1(k)0.07690.07690.07690.07690.07690.07690.07690.07690.07690.0769x2(k)0.13020.13020.13020.13020.13020.13020.13020.13020.13020.1302x3(k)0.1670.1670.1670.1670.1670.1670.1670.1670.1670.167x4(k)0.19260.19260.19260.19260.19260.19260.19260.

10、19260.19260.1926一直计算到29步得到：x(29)=(0.25,0.25)T。Seidel法的迭代列表如下所示k12345678910x1(k)0.07690.08280.08920.09610.10350.11140.120.12920.13920.1499x2(k)0.15550.16110.16660.1720.17730.18240.18720.19160.19570.1992x3(k)0.20250.20570.20870.21160.21420.21660.21890.2210.2230.2248x4(k)0.22650.22810.22960.2310.23230.

11、23350.23460.23560.23660.2375一直计算到17步得到：x(17)=(0.25,0.25)T。Sor迭代法的收敛速度取决于的值，取值恰当，便可以加速收敛，在本题中当=1.2时，收敛速度最快，其迭代列表如下k12345678910x1(k)0.09230.10080.11010.12030.13140.14350.15680.17130.18710.2043x2(k)0.19620.20330.21010.21630.22190.22670.23050.23310.23420.2335x3(k)0.23860.24040.24190.2430.24380.24440.244

12、90.24550.24630.2476x4(k)0.24750.24780.2480.24830.24850.24880.2490.24920.24930.2492一直计算到9步得到：x(19)=(0.25,0.25)T。比较以上三种方法可以看出Jacobi的收敛速度慢于Seidel，而Sor则在的值选取合适的时候，才能加速收敛。上机作业31、用Newton法与最速下降法求方程组x2+2y2-2=0 x2=y 在(0.8 , 0.7) 附近的根。解:1、Newton法：原理：由题目可得F(x)=,F(x)=,应用牛顿迭代法Fxk,ykz=-F(xk,y(k)x(k+1)y(k+1)=xky(k

13、)+z求解得到近似解。程序以及结果：x0=0.8;0.7;e=10(-5);N=2000;for k=1:Ny=-2*x0(1),4*x0(2); 2*x0(1),-1x0(1)2+2*x0(2)2-2; x0(1)2-x0(2);x=x0+y;if norm(y)=N disp(发散)endend 结果：近似解以及迭代次数 2、最速下降法：原理:设(x)=f12(x)+f22(x)=(x2+2y2-2)2+(x2-y)2，（1）求解出其在x(0)点的负梯度方向，G0 =grad(x(0)=-()T=-4x(0)1 *(x(0)12+2 x(0)22-2)+2 x(0)1*( x(0)12-

14、x(0)2); 8 x(0)2*( x(0)12+2 x(0)22-2)-2 (x(0)12- x(0)2);（2）然后求正数0, 使得(x(0)+ 0G0)= min(x(0)+0G0)( 00); (3) 取x(1) = x(0) + 0G0, 则满足(x(1) ) (x(0) );(4) x(1) x(0),重复上述过程，直到满足(x(k) )。程序以及结果：x0=0.8;0.7;e=10(-5);N=2000;sym afor k=1:NQ=-4*x0(1)*( x0(1)2+2*x0(2)2-2)+2*x0(1)*( x0(1)2-x0(2);8*x0(2)*( x0(1)2+2*x

15、0(2)2-2)-2* ( x0(1)2-x0(2);y=inline(x12+2*x22-2)2+(x12-x2)2,x1,x2);u= fminbnd(a)y(x0(1)+a*Q(1),x0(2)+a*Q(2),0,100);x=x0+u*Q;if (y(x(1),x(2)=N disp(发散)endend 结果：近似解以及迭代次数 上机作业4 1. 用幂法与反幂法求矩阵A的按模最大、最小特征值与对应的特征向量。A=解：a）幂法（求模最大特征值与对应特征向量）程序以及结果：A=4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2;v=1;1;1;1;eps=10(-5);N=

16、2000;lamda=0;err=1;k=1;while(keps)u=A*v;m j=max(abs(u);dc=abs(lamda-m);u=u/m;dv=norm(u-v);err=max(dc,dv);v=u;lamda=m;k=k+1; enddisp(lamda);disp(v); disp(k) 结果：模最大的特征值和对应的特征向量及迭代次数 b) 反幂法（求模最小特征值与对应特征向量）程序以及结果：A=4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2;v0=1;1;1;1;eps=10(-5);N=2000;m,i=max(abs(v0);lam0=m;u0=

17、v0/lam0;for k=1:Nv1=Au0;m1,i=max(abs(v1);lam1=m1;u1=v1/lam1;if abs(1/lam1-1/lam0)eps break endu0=u1;v0=v1;lam0=lam1;endlam=1/lam1;disp(lam);disp(u0);disp(k) 结果：模最大的特征值和对应的特征向量及迭代次数2. 用Householder变换求矩阵A的QR分解，并用QR方法做3次迭代A=解：程序以及结果：A=4,1,1,1;1,3,-1,1;1,-1,2,0;1,1,0,2;Q=eye(4);for k=1:3c=norm(A(k:4,k);w=zeros(4,1);w(k+1:4)=A(k+1:4,k);w(k)=A(k,k)+sign(A(k,k)*c;P=eye(4)-2*w*w/norm(w)2;Q=Q*P;A=P*A;enddisp(Q);disp(A)结果：矩阵Q以及R（A）QR法做3次迭代