数值分析上机作业(总)

《数值分析上机作业(总)》由会员分享，可在线阅读，更多相关《数值分析上机作业(总)（12页珍藏版）》请在装配图网上搜索。

1、.数值分析上机实验一、解线性方程组直接法教材49页14题追赶法程序如下：function x=followupn = rank;forifA=0 disp;return;endend; d = ones;a = ones;c = ones;for a=A; c=A; d=A;endd = A;for d=d - a/d*c; b=b - a/d*b;endx = b/d;fori=:-1:1 x = b-c*x/d;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、 -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计算结果：x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.11940.0477二、解线性方程组直接法教材49页15题程序如下：function tiaojianshuA=zeros;for j=1:1:nfor i=1:1:n A=;endendc=condd=rcond当n=5时c = 5.3615e+005d = 9.

3、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=jacobiif nargin=3 eps= 1.0e-6; M = 200;elseif nargin3 errorreturnelseif nargin =5 M = varargin1;endD=diagdiag; L=-tril;

4、 U=-triu; B=D;f=Db;x=B*x0+f;n=1; while norm=eps x0=x; x=B*x0+f; n=n+1;if=M disp;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计算结果：x = 1.0000 -2.0000 3.0000 -2.0000 1.0000n =67经过67次迭代,得到最终结果2用Gauss-S

5、eidel迭代法求：Gauss-Seidel迭代法程序如下：function x,n=gauseidelif nargin=3 eps= 1.0e-6; M = 200;elseif nargin = 4 M = 200;elseif nargin3 errorreturn;endD=diagdiag; L=-tril; U=-triu; G=U;f=b;x=G*x0+f;n=1; while norm=eps x0=x; x=G*x0+f; n=n+1;if=M disp;return;endend本题主程序如下：function gaosidiedaiA=10 1 2 3 4;1 9 -1

6、 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计算结果：x = 1.0000 -2.0000 3.0000 -2.0000 1.0000n =38经过38次迭代,得到最终结果。四、矩阵特征值与特征向量的计算教材100页13题幂法求最大特征值的程序：function l,v,s=pmethodif nargin=2 eps = 1.0e-6;endv = x0; M = 5000; m = 0; l = 0;for y = A*v; m = max; v = y/m

7、;ifabs l = m; s = k; return;elseif 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求解结果：l = 343.0000v = -0.6667 -2.0000 0 1.0000s = 114结论：经过114次迭代,求得此矩阵的最大特征值为343.0000,及其对应特征向量为-0.6667;-2.0000;0;1.0000

8、五、函数逼近教材164页16题本题采用最小二乘法进行拟合：线性最小二乘法程序如下：function a,b=LZXECiflength = length n = length; else disp;return;endA = zeros;A = n;B = zeros;for i=1:n A = A + x*x; A = A + x; B = B + x*y; B = B + y;endA = A;s = AB;a = s;b = s;首先利用1/y代替y,1/x代替x并采用线性最小二乘法求出a与b：function zuixiaoerchengx1=2 3 5 6 7 9 10 11 12

9、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计算结果：b = 8.4169e-004a =0.0090绘制图形：fplotx/,2,20grid ontitle六、数值微分与数值积分教材207页26题本题采用高斯勒让德求积公式求解：高斯勒让德求积公式程序如下：function I = IntGaussifn AK = 0; XK = 0;else XK1=

10、/2*XK+/2; I=/2*sumAK.*subssym,findsym,XK1;Endta = /2;tb = /2;switch ncase 0, I=2*ta*subssym,findsymsym,tb;case 1, I=ta*subssym,findsymsym,ta*0.5773503+tb+. subssym,findsymsym,-ta*0.5773503+tb;case 2, I=ta*0.55555556*subssym,findsymsym,ta*0.7745967+tb+. 0.55555556*subssym,findsymsym,-ta*0.7745967+tb+

11、. 0.88888889*subssym,findsymsym,tb;case 3, I=ta*0.3478548*subssym,findsymsym,ta*0.8611363+tb+. 0.3478548*subssym,findsymsym,-ta*0.8611363+tb+. 0.6521452*subssym,findsymsym,ta*0.3398810+tb. +0.6521452*subssym,findsymsym,-ta*0.3398810+tb;case 4, I=ta*0.2369269*subssym,findsymsym,ta*0.9061793+tb+. 0.23

12、69269*subssym,findsymsym,-ta*0.9061793+tb+. 0.4786287*subssym,findsymsym,ta*0.5384693+tb. +0.4786287*subssym,findsymsym,-ta*0.5384693+tb+. 0.5688889*subssym,findsymsym,tb;end本题计算程序如下采用4个节点：function shuzhijifena=1;b=1;for s=-5:0.05:5 q1=IntGausscos,0,s,4; q2=IntGausssin,0,s,4; a=a+1; b=b+1;endplot;绘制

13、图形：七、非线性方程及非线性方程组的解法本题采用牛顿法进行求解：牛顿法解非线性方程程序如下：function r,n=mulNewtonif nargin=2 eps=1.0e-4;endx0 = transpose;Fx = subsF,findsym,x0;var = findsym;dF = Jacobian;dFx = subsdF,findsym,x0;r=x0-inv*Fx;n=1;tol=1;while toleps x0=r; Fx = subsF,findsym,x0; dFx = subsdF,findsym,x0; r=x0-inv*Fx; tol=norm; n=n+1

14、;if100000 disp;return;endend本题解决方案如下：首先,绘制此方程的图形,大概确定其与X轴的交点位置。由于,可以得出因此绘制程序如下：ezplotlog/-x/,-772,772,-10,10;grid on得到图形如下图所示：经过放大后,发现图形与x轴的交点接近处。计算非零根：令为牛顿法接非线性方程的初值。程序如下：syms xf=log/-x/;x0=-765r,n=mulNewton解得：r = -767.3861x0=765r,n=mulNewton解得：r =767.3861结论：此方程的两个非零根分别为：八、常微分方程数值解法教材266页19题本题分别采用四

15、阶ADAMS预测校正算法和经典RK法进行求解：四阶ADAMS预测校正算法如下：function y = DEYCJZ_yds format long;N = /h;y = zeros;x = a:h:b;y = y0;y = y0+h*Funvalf,varvec,x y;y = y+h*Funvalf,varvec,x y;y = y+h*Funvalf,varvec,x y;for i=5:N+1 v1 = Funvalf,varvec,x y; v2 = Funvalf,varvec,x y; v3 = Funvalf,varvec,x y; v4 = Funvalf,varvec,x

16、y; t = y + h*/24; ft = Funvalf,varvec,x t; y = y+h*/24;end经典RK算法程序如下：function y = DELGKT4_lungkutaformat long;N = /h;y = zeros;y = y0;x = a:h:b;var = findsym;for i=2:N+1 K1 = Funvalf,varvec,x y; K2 = Funvalf,varvec,x+h/2 y+K1*h/2; K3 = Funvalf,varvec,x+h/2 y+K2*h/2; K4 = Funvalf,varvec,x+h y+h*K3; y

17、 = y+h*/6;end其中FUNVAL函数程序如下：function fv = Funvalvar = findsym;if length 4if var = varvec fv = subsf,varvec,varval;else fv = subsf,varvec,varval;endelse fv = subs;end本题解决方案如下步长h=0.1：程序：function changweifensyms xy;z = -y+2*cos;yy1 = DELGKT4_lungkuta;yy2 = DEYCJZ_yds;a=0:0.1:pi;yy3 = cos+sin;plot;grid

18、on;hold on;plot;plot;legend整体图形：局部放大图形：继续放大：数据结果：yy1yy2yy3yy1-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.35700757

19、91.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.2

20、5739E-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.9

21、703732280.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.3

22、58363262-6.11484E-070.0025108260.2199947470.2222666510.219995287-5.39985E-070.0022713640.0794287280.0814838660.079429191-4.62442E-070.002054676-0.061930915-0.060071936-0.061930535-3.7969E-070.001858599-0.202671764-0.200990293-0.202671471-2.92614E-070.001681178-0.341387584-0.339866738-0.341387382-2.0

23、2132E-070.001520644-0.476692371-0.475316868-0.476692262-1.09196E-070.001375394-0.607234205-0.605990211-0.607234191-1.47751E-080.00124398-0.731708756-0.730583746-0.7317088368.01498E-080.00112509-0.848872314-0.847854952-0.8488724891.74596E-070.001017537-0.95755422-0.95663424-0.9575544882.6759E-070.000920247结论：通过图形和数据结果可以看出,在本题中利用经典RK方法获得解的精确度比4阶ADAMS方法要高很多。数值计算方法上机实验 _刘天博 _2120100204 学院：机电学院.