计算方法作业.doc

应用计算方法上机作业应用计算方法学院：自动化学院班级：硕1104班姓名：陈兆辉计算方法2班教师：张晓丹目录上机作业12上机作业25上机作业310上机作业412上机作业516上机作业618上机作业719上机作业11.分别用不动带你迭代法与Newton法求解方程2x- ex +3=0的正根和负根解：（1）用不动点迭代法求方程的正根1）迭代公式：，初值p0=1.02）程序：k=1;Tol=0.00001;p0=1.0;N=100;while k=Np=log(2*p0+3); if abs(p-p0)Tol break; end k=k+1;p0=p;enddisp(p);disp(k)3)显示结果： P= 1.9239K=11（2） 用不动点迭代法求解方程的负根1）迭代公式： ，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k=Np=(exp(p0)-3)/2; if abs(p-p0)Tol break; end k=k+1;p0=p;enddisp(p);disp(k)3）显示结果： P= -1.3734k= 7 (3)用牛顿法求方程的正根1）迭代公式：，初值p0=1.02） 程序：k=1;Tol=0.00001;p0=1.0;N=100;while k=N p=p0-(2*p0-exp(p0)+3)/(2-exp(p0); if abs(p-p0)Tol break; end k=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P=1.9239 K=8（4）用牛顿法求解方程的负根1）迭代公式：，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k=N p=p0-(2*p0-exp(p0)+3)/(2-exp(p0); if abs(p-p0)Tol break; end k=k+1;p0=p;enddisp(p);disp(k)3）显示结果： P= -1.3734 K= 44）结果分析：不动点迭代法是求解非线性方程的主要方法，其中牛顿法应用最广泛，它求解方程的单根时具有二阶收敛性，但是它要求较好的初始近似值，牛顿法比普通的迭代法收敛速度快，能较少的迭代达到理想的值。2、用Newton法求解方程x-sinx=0的根，再用steffensens method加速其收敛。（1）用牛顿法求解方程的根1）Newton迭代方程：，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k=N p=p0-(p0-sin(p0)/(1-cos(p0); if abs(p-p0)Tol break; end k=k+1;p0=p;enddisp(p);disp(k)3）结果显示： P= -1.9449e-004K= 21（2）用steffensens method 加速Newton法收敛1）牛顿迭代方程 ，初值p0=-1.02）steffen加速收敛程序：n=1;Tol=0.00001;p0=-1.0;p(1)=p0;N=100;while n=N for k=1:2 p(k+1)=p(k)-(p(k)-sin(p(k)/(1-cos(p(k); end p1=p(1)-(p(2)-p(1)2/(p(3)-2*p(2)+p(1); f0=p1-sin(p1); if abs(f0)Tol break; end n=n+1;p(1)=p1;enddisp(p1);disp(n)3）结果显示：P1=-0.0355n=14）结果分析：steffen 加速极大地改善了原迭代的收敛性质，它加速了牛顿法的收敛速度，仅迭代一次就超过了原来的迭代结果！上机作业21、分别用jacobi，seidel，sor（=1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b，这里:， (1)、jacobi方法求解方程组Ax=b1）Jacobi不动点方程：Jacobi迭代方程：2）Jacobi程序:A=-12*eye(10)+ones(10);b=-ones(10,1);X0=zeros(10,1);e=1.0e-5;N=100;X=X0;for K=1:N for i=1:10 X(i)=(b(i)-A(i,:)*X0)/A(i,i)+X0(i); end if norm(X-X0)e disp(X);disp(K); return; end X0=X;enddisp(发散);3）程序结果：X=0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 K= 53（2）seidel方法求解方程组Ax=b1）seidel迭代方程：2）seidel迭代程序：A=-12*eye(10)+ones(10);b=-ones(10,1);X0=zeros(10,1);e=1.0e-5;N=100;X=X0;for K=1:N for i=1:10 X(i)=(b(i)-A(i,:)*X)/A(i,i)+X(i); end if norm(X-X0)e disp(X);disp(K); return; end X0=X;enddisp(发散);3）程序结果：X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000K=29（3）sor（=1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b1）sor迭代方程2）sor迭代程序：A=-12*eye(10)+ones(10);b=-ones(10,1);X0=zeros(10,1);e=1.0e-5;N=100;X=X0;w=1.1;for K=1:N for i=1:10 X(i)=w*(b(i)-A(i,:)*X)/A(i,i)+X(i); end if norm(X-X0)e disp(X);disp(K); return; end X0=X;enddisp(发散);3）程序结果：w = 1.1000X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000K= 24w = 1.2000X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 K=19w = 1.3000X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 K=14w = 1.4000X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 K=13w = 1.5000X= 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 K=174）结果分析：从计算结果可以看出，若Jacobi迭代法和Gauss-Seidel迭代法都收敛时，Gauss-Seidel迭代法比用于Jacobi迭代法求解收敛速度快。在该方程组中， Jacobi迭代法次数为53次，而Gauss-Seidel迭代法的迭代次数为29次，大大加快了收敛速度。但对不收敛的情况失效。应用SOR时选取合适的松弛因子是很关键的，从解该方程组程序及结果对照表中可以看出，若松弛因子w选取的好，可以加速其收敛，当w=1.11.4时，随着w的值的增大，迭代收敛速度加快，但当w=1.5时，迭代法需要迭代17次，说明松弛因子不是越大越好。上机作业31、用Newton法与最速下降法求方程组在（0.8,0.7）附近的根（1）Newton法：1）解：由题目所给方程组可以得到取初始向量，应用牛顿迭代法2）牛顿迭代法程序：x0=0.8;0.7;e=1.0e-5;N=100;for k=1:N Fx0=x0(1)*x0(1)+2*x0(2)*x0(2)-2 ; x0(1)*x0(1)-x0(2); FFx0=2*x0(1) 4*x0(2); 2*x0(1) -1; y=-FFx0Fx0; x=x0+y; if norm(y)e disp(x);disp(k); return; end x0=x;enddisp(发散)3）程序结果：X=0.88360.7808K=4（2）最速下降法1）解：由原方程组构造一个模函数其中2）程序：x0=0.8;0.7;e=1.0e-5;N=100;function f=2*x14+4*x24+4*x1*x22-4*x12-7*x22-2*x12*x2+4for k=1:N g=8*x13+8*x1*x22-8*x1-4*x1*x2;16*x23+8*x12*x2-14*x2-2*x12 h(t)=2*(x1-t*g)4+4*(x2-t*g)4+4*(x1-t*g)4*(x2-t*g)2-4*(x1-t*g)2-7*(x2-t*g)2-2*(x1-t*g)2*(x2-t*g)+4 u=fmin(h,0,100) x=x-u*g if fe disp(x);disp(k);return; end x0=x enddisp(发散)3）程序结果：X=0.74130.6470K=74）结果分析：最速下降法计算简便，但收敛速度比牛顿法慢。上机作业41、用幂法与反幂法求矩阵A的按模最大、最小特征值与对应的特征向量（1） 幂法求矩阵A的按模最大特征值与对应的特征向量1）程序： A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2; k=1; N=100; eps=1e-006; err=1; v=1;1;1;1; lamda=0; 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+1end2）程序结果：v =1.0000;0.6180;0.1180;0.5000;lamda=5.2361k=32（2） 反幂法求矩阵A的按模最小特征值与对应的特征向量1）反幂法的计算公式： 2）程序：A=4 1 1 1;1 3 -1 1;1 -1 2 0;1 1 0 2; k=1; N=100; eps=1e-006; err=1; v=1;1;1;1; lamda=0;while (keps) u=(inv(A)*v; m=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(u);disp(k)3）显示结果V= -0.4721； 0.7639；1.0000；-0.2361；K=244）程序结果：lamda=1.3090,由于本程序是对A求逆，所以本程序所求的特征值对应反幂法中最小的特征值，即：反幂法所求的最小特征值 lamda1=1/lamda=1/1.3090=0.7639,但特征向量是相同的。2、用Householder变换求矩阵A的QR分解，并用QR方法做3次迭代1）解：Householder变换公式：利用Householder矩阵求矩阵A的QR分解公式：，2）程序：A=4 1 1 1;1 3 -1 1;-1 -1 2 0;1 1 0 2;n=size(A)Q=eye(n)for k=1:n-1 c=norm(A(k:n,k); w=zeros(n,1); w(k+1:n)=A(k+1:n,k) w(k)=A(k,k)+sign(A(k,k)*c; P=eye(n)-2*w*w/norm(w)2 Q=Q*P A=P*A R=P*Q;end3）程序结果： n = 4 4Q =1 0 0 00 1 0 00 0 1 00 0 0 1w = 0 1 -1 1P = -0.9177 -0.2294 0.2294 -0.2294 -0.2294 0.9726 0.0274 -0.0274 0.2294 0.0274 0.9726 0.0274 -0.2294 -0.0274 0.0274 0.9726Q = -0.9177 -0.2294 0.2294 -0.2294 -0.2294 0.9726 0.0274 -0.0274 0.2294 0.0274 0.9726 0.0274 -0.2294 -0.0274 0.0274 0.9726A = -4.3589 -2.0647 -0.2294 -1.6059 -0.0000 2.6334 -1.1471 0.6882 0.0000 -0.6334 2.1471 0.3118 0 0.6334 -0.1471 1.6882w = 0 0 -0.6334 0.6334P = 1.0000 0 0 0 0 -0.9467 0.2277 -0.2277 0 0.2277 0.9734 0.0266 0 -0.2277 0.0266 0.9734Q = -0.9177 0.3217 0.1650 -0.1650 -0.2294 -0.9083 0.2474 -0.2474 0.2294 0.1892 0.9536 0.0464 -0.2294 -0.1892 0.0464 0.9536A = -4.3589 -2.0647 -0.2294 -1.6059 0.0000 -2.7815 1.6084 -0.9650 0.0000 0.0000 1.8248 0.5051 0.0000 0 0.1752 1.4949w = 0 0 0 0.1752P = 1.0000 0 0 0 0 1.0000 0 0 0 0 -0.9954 -0.0956 0 0 -0.0956 0.9954Q = -0.9177 0.3217 -0.1484 -0.1800 -0.2294 -0.9083 -0.2227 -0.2700 0.2294 0.1892 -0.9537 -0.0450 -0.2294 -0.1892 -0.1373 0.9448A = -4.3589 -2.0647 -0.2294 -1.6059 0.0000 -2.7815 1.6084 -0.9650 -0.0000 -0.0000 -1.8332 -0.64570.0000 -0.0000 0 1.4397上机作业5目的：观察Lagrange插值的Runge现象，对函数 进行lagrange插值取不同的等分数n（n=5,10）；将区间-1, 1n等分，取等距节点，把插值多项式的曲线在同一张图上进行比较。解：lagrange插值程序：%绘制原函数图形figure(1)tt=-1:0.01:1;y=1./(1+25*tt.2);plot(tt,y,k)grid on;hold on; %Lagrange插值n=5;x=-1:2/n:1;%五等分y=1./(1+25*x.2);yy=0;xx=-1:0.1:1;for k=1:n+1 t=1; for i=1:n+1; if i=k; t=t.*(xx-x(i)/(x(k)-x(i); end end yy=yy+t*y(k);endplot(xx,yy,b:.)hold on; n=10;x=-1:2/n:1;%十等分y=1./(1+25*x.2);yy=0;xx=-1:0.1:1;for k=1:n+1 t=1; for i=1:n+1; if i=k; t=t.*(xx-x(i)/(x(k)-x(i); end end yy=yy+t*y(k);endplot(xx,yy,r-*)title(Lagrange插值 -.-5五等分 -*-10十等分 -原函数)结果分析：从图中可以看出，在靠近-1和1时，余项会随n的增加而增加，产生较大的误差，说明靠节点数增多而余项不会随n的增加而趋近于0，需要考虑其他方法。上机作业6分别用Romberg求积公式计算积分，使精度达到解：1）Romberg求积公式：2）Romberg求积分程序：a=1;b=3;eps=1e-006;M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(exp(a)*cos(a)+exp(b)*cos(b)/2; while(erreps) J=J+1; h=h/2; x=a+h:2*h:b-h; R(J+1,1)=R(J,1)/2+h*sum(exp(x).*cos(x); for K=1:min(3,J) R(J+1,K+1)=R(J+1,K)+(R(J+1,K)-R(J,K)/(4K-1); end if(J3) err=abs(R(J+1,4)-R(J,4); end endquad=R(J,4)3）程序结果：quad=-10.4031上机作业71.分别用Euler方法（h=0.025）、改进的Euler方法（h=0.05）、经典的R-K方法（h=0.1）求初值问题的数值解计算 0.1,0.2,0.3,0.4,0.5处解的近似值。并与精确解比较，分析误差。解：（1）Euler方法：1）欧拉公式：2）程序：x0=0;y0=0.5;h=0.025;n=5;x=x0;y=y0;for i=1:n y=y+h*(y-x2+1) x