《数值分析报告上机作业总.doc》由会员分享，可在线阅读，更多相关《数值分析报告上机作业总.doc（15页珍藏版）》请在三一办公上搜索。
1、word数值分析上机实验一、解线性方程组直接法教材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
2、,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
3、0 0 0 0 -2 5;b=220/27;0;0;0;0;0;0;0;x=followup(A,b)计算结果：x =二、解线性方程组直接法教材49页15题程序如下：function tiaojianshu(n)A=zeros(n);for j=1:1:nfor 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 =当n=10时c = 8.6823e+011d =当n=20时c = 3.4205e+d =备注：对于病态矩阵A来说，d为接近0的数；对于非病态矩阵A来说，d为接近1的数。三、解线性
4、方程组的迭代法教材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
5、;-27;14;-17;12;x0=0;0;0;0;0;x,n=jacobi(A,b,x0)计算结果：x =n =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此题主程
6、序如下：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 =n =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
7、=1:M) y = A*v; m = max(y); v = y/m;if(abs(m - l)eps) l = m; s = k; return;elseif(k=M) disp(收敛速度过慢); l = m; s = M;else l = m;endendend求解此题程序如下：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 =v = 0s = 114结论：经过114次迭代，求得此矩阵的最大特征值为343.0000，
8、与其对应特征向量为-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;endA = 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);
9、 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)计算
10、结果：b =a =绘制图形：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(迭代步数太多，可能不收敛);
11、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)解得：x0=765r,n=mul
12、Newton(f,x0)解得：结论：此方程的两个非零根分别为：八、常微分方程数值解法教材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*Funva
13、l(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
14、经典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,varv
15、ec,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) 4if 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：程序：fu
16、nction changweifensyms xy;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-yy311100结论：通过图形和数据结果可以看出，在此题中利用经典RK方法获得解的准确度比4阶ADAMS方法要高很多。15 / 15
