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

数值分析上机实验一、解线性方程组直接法（教材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.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.11940.0477二、解线性方程组直接法（教材49页15题）程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:n for i=1:1:n A(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的数。三、解线性方程组的迭代法（教材74页14题）（1）用Jacobi迭代法求：Jacobi迭代法程序如下：function x,n=jacobi(A,b,x0,eps,varargin)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可能不收敛); return; endend本题主程序如下：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.0000 3.0000 -2.0000 1.0000n =67经过67次迭代，得到最终结果（2）用Gauss-Seidel迭代法求：Gauss-Seidel迭代法程序如下：function x,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin = 4 M = 200;elseif nargin=eps x0=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.0000 3.0000 -2.0000 1.0000n =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)eps) l = m; s = k; return; else if(k=M) disp(收敛速度过慢); l = m; s = M; else l = m; end endend求解本题程序如下：function mifaA=190 66 -84 30;66 303 42 -36;336 -168 147 -112;30 -36 28 291;x0=0 0 0 1;l,v,s=pmethod(A,x0)求解结果：l = 343.0000v = -0.6667 -2.0000 0 1.0000s = 114结论：经过114次迭代，求得此矩阵的最大特征值为343.0000，及其对应特征向量为-0.6667;-2.0000;0;1.0000五、函数逼近（教材164页16题）本题采用最小二乘法进行拟合：线性最小二乘法程序如下：function a,b=LZXEC(x,y)if(length(x) = length(y) n = length(x); else disp(x和y的维数不相等); return;end A = zeros(2,2);A(2,2) = n;B = zeros(2,1);for i=1:n A(1,1) = A(1,1) + x(i)*x(i); A(1,2) = A(1,2) + x(i); B(1,1) = B(1,1) + x(i)*y(i); B(2,1) = B(2,1) + y(i);endA(2,1) = A(1,2);s = AB;a = s(1);b = s(2);首先利用1/y代替y，1/x代替x并采用线性最小二乘法求出a与b：function zuixiaoerchengx1=2 3 5 6 7 9 10 11 12 14 16 17 19 20;y1=106.42 108.26 109.58 109.50 109.86 110.00 109.93 110.59 110.60 110.72 110.90 110.76 111.10 111.30;x2=1./x1;y2=1./y1;b,a=LZXEC(x2,y2)计算结果：b = 8.4169e-004a =0.0090绘制图形：fplot(x/(0.0090*x+8.4169e-004),2,20)grid ontitle(最小二乘拟合)六、数值微分与数值积分（教材207页26题）本题采用高斯勒让德求积公式求解：高斯勒让德求积公式程序如下：function I = IntGauss(f,a,b,n,AK,XK)if(neps x0=r; Fx = subs(F,findsym(F),x0); dFx = subs(dF,findsym(dF),x0); r=x0-inv(dFx)*Fx; tol=norm(r-x0); n=n+1; if(n100000) disp(迭代步数太多，可能不收敛); return; endend本题解决方案如下：首先，绘制此方程的图形，大概确定其与X轴的交点位置。由于，可以得出因此绘制程序如下：ezplot(log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918),-772,772,-10,10);grid on得到图形如下图所示：经过放大后，发现图形与x轴的交点接近处。计算非零根：令为牛顿法接非线性方程的初值。程序如下：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结论：此方程的两个非零根分别为：八、常微分方程数值解法（教材266页19题）本题分别采用四阶ADAMS预测校正算法和经典RK法进行求解：四阶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,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 = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 K1 = 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 fv = Funval(f,varvec,varval)var = findsym(f);if length(var) 4 if var(1) = varvec(1) fv = subs(f,varvec(1),varval(1); else fv = subs(f,varvec(2),varval(2); endelse fv = 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 = 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,准确解)整体图形：局部放大图形：继续放大：数据结果：yy1(RK)yy2(ADAMS)yy3(准确解)yy1-yy3yy2-yy3111001.0948374641.11.094837582-1.18389E-070.0051624181.1787356781.1890008331.178735909-2.30293E-070.0102649241.2508563611.2661140651.250856696-3.351E-070.015257371.3104789041.3242156421.310479336-4.32217E-070.0137363051.3570075791.3694466421.3570081-5.21089E-070.0124385421.3899774871.4012422871.389978088-6.012E-070.0112641991.4090592021.4192520421.409059875-6.72089E-070.0101921671.4140620671.4232851431.4140628-7.33353E-070.0092223431.4049360931.413281631.404936878-7.84658E-070.0083447531.3817724651.3893238971.381773291-8.25739E-070.0075506071.3448026251.3516354611.344803481-8.56415E-070.006831981.2943959641.3005785271.29439684-8.76584E-070.0061816861.2310561281.2366502381.231057014-8.86229E-070.0055932241.1554159871.1604775841.155416873-8.85423E-070.0050607111.0682313141.0728110131.068232188-8.74325E-070.0045788250.9703732280.9745168330.970374081-8.53183E-070.0041427520.8628194940.8665684520.862820316-8.22334E-070.0037481360.7466447540.750036570.746645536-7.82199E-070.0033910340.6230097880.6260784020.623010521-7.33279E-070.0030678820.4931499140.4959260420.49315059-6.76156E-070.0027754520.3583626510.3608740880.3583632