数值分析上机作业(总)

1、-作者xxxx-日期xxxx数值分析上机作业(总)【精品文档】数值分析上机实验一、解线性方程组直接法（教材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:

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

3、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 =65二、解线性方程组直接法（教材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 =当n=10时c = 8.6823e+011d =当n=20时c = 3.4205e+022d =

4、备注：对于病态矩阵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

5、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 = 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+

6、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 =n =38经过38次迭代，得到最终结果。四、矩阵特征值与特征向量的计算（教材100页13题）幂法求最大特征值的程序：function l,v,s=pmethod(A,x0,eps)if n

7、argin=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=

8、pmethod(A,x0)求解结果：l =v = 0s = 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(

9、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 1

10、10.72 110.90 110.76 111.10 111.30;x2=1./x1;y2=1./y1;b,a=LZXEC(x2,y2)计算结果：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);

11、 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)

12、/(513-0.6651*x)-x/(1400*0.0918);x0=-765r,n=mulNewton(f,x0)解得：x0=765r,n=mulNewton(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,

13、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)/

14、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

15、)+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