数值分析上机作业(总)(13页).doc
-数值分析上机作业(总)-第 13 页数值分析上机实验一、解线性方程组直接法（教材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 =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 =备注：对于病态矩阵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<3 error returnelseif nargin =5 M = varargin1;end D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); B=D(L+U);f=Db;x=B*x0+f;n=1; while norm(x-x0)>=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 =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<3 error return;end D=diag(diag(A); L=-tril(A,-1); U=-triu(A,1); G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1; while norm(x-x0)>=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 =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=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 =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(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 =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(n<5 && nargin = 4) AK = 0; XK = 0;else XK1=(b-a)/2)*XK+(a+b)/2); I=(b-a)/2)*sum(AK.*subs(sym(f),findsym(f),XK1);Endta = (b-a)/2;tb = (a+b)/2;switch n case 0, I=2*ta*subs(sym(f),findsym(sym(f),tb); case 1, I=ta*(subs(sym(f),findsym(sym(f),ta*0.5773503+tb)+. subs(sym(f),findsym(sym(f),-ta*0.5773503+tb); case 2, I=ta*(0.55555556*subs(sym(f),findsym(sym(f),ta*0.7745967+tb)+. 0.55555556*subs(sym(f),findsym(sym(f),-ta*0.7745967+tb)+. 0.88888889*subs(sym(f),findsym(sym(f),tb); case 3, I=ta*(0.3478548*subs(sym(f),findsym(sym(f),ta*0.8611363+tb)+. 0.3478548*subs(sym(f),findsym(sym(f),-ta*0.8611363+tb)+. 0.6521452*subs(sym(f),findsym(sym(f),ta*0.3398810+tb). +0.6521452*subs(sym(f),findsym(sym(f),-ta*0.3398810+tb); case 4, I=ta*(0.2369269*subs(sym(f),findsym(sym(f),ta*0.9061793+tb)+. 0.2369269*subs(sym(f),findsym(sym(f),-ta*0.9061793+tb)+. 0.4786287*subs(sym(f),findsym(sym(f),ta*0.5384693+tb). +0.4786287*subs(sym(f),findsym(sym(f),-ta*0.5384693+tb)+. 0.5688889*subs(sym(f),findsym(sym(f),tb);end本题计算程序如下（采用4个节点）：function shuzhijifena=1;b=1;for s=-5:0.05:5 q1(1,a)=IntGauss('cos(1/2*x2)',0,s,4); q2(1,b)=IntGauss('sin(1/2*x2)',0,s,4); a=a+1; b=b+1;endplot(q1,q2);绘制图形：七、非线性方程及非线性方程组的解法本题采用牛顿法进行求解：牛顿法解非线性方程程序如下：function r,n=mulNewton(F,x0,eps)if nargin=2 eps=1.0e-4;endx0 = transpose(x0);Fx = subs(F,findsym(F),x0);var = findsym(F);dF = Jacobian(F,var);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx;n=1;tol=1;while tol>eps 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(n>100000) 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)解得：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,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-yy311100结论：通过图形和数据结果可以看出，在本题中利用经典RK方法获得解的精确度比4阶ADAMS方法要高很多。数值计算方法上机实验 姓名：刘天博 学号：2120100204 学院：机电学院
