数值分析上机作业(总).doc

1、数值分析上机实验一、解线性方程组直接法（教材 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

2、)*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

3、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.0894

4、e-013当n=20时c =3.4205e+022d =8.1226e-024备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。三、解线性方程组的迭代法（教材 74 页 14 题）（1）用 Jacobi 迭代法求：Jacobi 迭代法程序如下：function x,n=jacobi(A,b,x0,eps,varargin)if nargin=3eps= 1.0e-6;M = 200;elseif nargin=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp(Warning: 迭代次数太多，可能不收敛 );return;endend本题主程序如

5、下：function yakebidiedaiA=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;b=12;-27;14;-17;12;x0=0;0;0;0;0;x,n=jacobi(A,b,x0)计算结果：x =1.0000-2.00003.0000-2.00001.0000n =67经过 67 次迭代，得到最终结果（2）用 Gauss-Seidel 迭代法求：Gauss-Seidel 迭代法程序如下：function x,n=gauseidel(A,b,x0,eps,M)if nargin=3eps= 1.0e-6;

6、M = 200;elseif nargin = 4M = 200;elseif nargin=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp(Warning: 迭代次数太多，可能不收敛 );return;endend本题主程序如下：function gaosidiedaiA=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;b=12;-27;14;-17;12;x0=0;0;0;0;0;x,n=gauseidel(A,b,x0)计算结果：x =1.0000-2.00003.0000-2.00001.0

7、000n =38经过 38 次迭代，得到最终结果。四、矩阵特征值与特征向量的计算（教材 100 页 13 题）幂法求最大特征值的程序：function l,v,s=pmethod(A,x0,eps)if nargin=2 eps = 1.0e-6;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)epsx0=r;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx; tol=norm

8、(r-x0);n=n+1;if(n100000) disp(迭代步数太多，可能不收敛);return;endend本题解决方案如下：首先，绘制此方程的图形，大概确定其与 X 轴的交点位置。由于 ，可以得出5130.6x71.3.x因此绘制程序如下：ezplot(log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918),-772,772,-10,10);grid on得到图形如下图所示：-600 -400 -200 0 200 400 600-10-8-6-4-20246810xlog(513+0.6651 x)/(513-0.6651 x)-x/(140

9、0 0.0918)经过放大后，发现图形与 x 轴的交点接近 处。75计算非零根：令 为牛顿法接非线性方程的初值。765x程序如下：syms xf=log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918);x0=-765r,n=mulNewton(f,x0)解得：r = -767.3861x0=765r,n=mulNewton(f,x0)解得：r =767.3861结论：此方程的两个非零根分别为： 12 76.38;76.381;xx八、常微分方程数值解法（教材 266 页 19 题）本题分别采用四阶 ADAMS 预测校正算法和经典 RK 法进行求解：四阶

10、 ADAMS 预测校正算法如下：function y = DEYCJZ_yds (f, h,a,b,y0,varvec) format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;y(2) = y0+h*Funval(f,varvec,x(1) y(1);y(3) = y(2)+h*Funval(f,varvec,x(2) y(2);y(4) = y(3)+h*Funval(f,varvec,x(3) y(3);for i=5:N+1 v1 = Funval(f,varvec,x(i-4) y(i-4);v2 = Funval(f

11、,varvec,x(i-3) y(i-3);v3 = Funval(f,varvec,x(i-2) y(i-2);v4 = Funval(f,varvec,x(i-1) y(i-1);t = y(i-1) + h*(55*v4 - 59*v3 + 37*v2 - 9*v1)/24;ft = Funval(f,varvec,x(i) t);y(i) = y(i-1)+h*(9*ft+19*v4-5*v3+v2)/24;end经典 RK 算法程序如下：function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y

12、 = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1K1 = Funval(f,varvec,x(i-1) y(i-1);K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2);K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2);K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3);y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;end其中FUNVAL函数程序如下：function f

13、v = Funval(f,varvec,varval)var = findsym(f);if length(var) 4if var(1) = varvec(1)fv = subs(f,varvec(1),varval(1);elsefv = subs(f,varvec(2),varval(2);endelsefv = subs(f,varvec,varval);end本题解决方案如下（步长 h=0.1）：程序：function changweifensyms x y;z = -y+2*cos(x);yy1 = DELGKT4_lungkuta(z,0.1,0,pi,1,x y);yy2 =

14、DEYCJZ_yds(z,0.1,0,pi,1,x y);a=0:0.1:pi;yy3 = cos(a)+sin(a);plot(a,yy1,r*);grid on;hold on;plot(a,yy2,b+);plot(a,yy3,-go);legend(RK,Adams,准确解)整体图形：0 0.5 1 1.5 2 2.5 3 3.5-1-0.500.511.5RKAdams乘乘乘局部放大图形：2.798 2.799 2.8 2.801 2.802 2.803-0.6095-0.609-0.6085-0.608-0.6075-0.607-0.6065-0.606-0.6055-0.605R

15、KAdams乘 乘 乘继续放大：2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8-0.6072-0.6072-0.6072-0.6072-0.6072-0.6072RKAdams乘 乘 乘数据结果：yy1(RK) yy2(ADAMS) yy3(准确解) yy1-yy3 yy2-yy31 1 1 0 01.094837464 1.1 1.094837582 -1.18389E-07 0.0051624181.178735678 1.189000833 1.178735909 -2.30293E-07 0.0102649241.250856361 1.266114065 1.2508

16、56696 -3.351E-07 0.015257371.310478904 1.324215642 1.310479336 -4.32217E-07 0.0137363051.357007579 1.369446642 1.3570081 -5.21089E-07 0.0124385421.389977487 1.401242287 1.389978088 -6.012E-07 0.0112641991.409059202 1.419252042 1.409059875 -6.72089E-07 0.0101921671.414062067 1.423285143 1.4140628 -7.

17、33353E-07 0.0092223431.404936093 1.41328163 1.404936878 -7.84658E-07 0.0083447531.381772465 1.389323897 1.381773291 -8.25739E-07 0.0075506071.344802625 1.351635461 1.344803481 -8.56415E-07 0.006831981.294395964 1.300578527 1.29439684 -8.76584E-07 0.0061816861.231056128 1.236650238 1.231057014 -8.862

18、29E-07 0.0055932241.155415987 1.160477584 1.155416873 -8.85423E-07 0.0050607111.068231314 1.072811013 1.068232188 -8.74325E-07 0.0045788250.970373228 0.974516833 0.970374081 -8.53183E-07 0.0041427520.862819494 0.866568452 0.862820316 -8.22334E-07 0.0037481360.746644754 0.75003657 0.746645536 -7.8219

19、9E-07 0.0033910340.623009788 0.626078402 0.623010521 -7.33279E-07 0.0030678820.493149914 0.495926042 0.49315059 -6.76156E-07 0.0027754520.358362651 0.360874088 0.358363262 -6.11484E-07 0.0025108260.219994747 0.222266651 0.219995287 -5.39985E-07 0.0022713640.079428728 0.081483866 0.079429191 -4.62442

20、E-07 0.002054676-0.061930915 -0.060071936 -0.061930535 -3.7969E-07 0.001858599-0.202671764 -0.200990293 -0.202671471 -2.92614E-07 0.001681178-0.341387584 -0.339866738 -0.341387382 -2.02132E-07 0.001520644-0.476692371 -0.475316868 -0.476692262 -1.09196E-07 0.001375394-0.607234205 -0.605990211 -0.6072

21、34191 -1.47751E-08 0.00124398-0.731708756 -0.730583746 -0.731708836 8.01498E-08 0.00112509-0.848872314 -0.847854952 -0.848872489 1.74596E-07 0.001017537-0.95755422 -0.95663424 -0.957554488 2.6759E-07 0.000920247结论：通过图形和数据结果可以看出，在本题中利用经典 RK 方法获得解的精确度比 4 阶 ADAMS 方法要高很多。数值计算方法上机实验姓名：刘天博学号：2120100204学院：机电学院