数值分析上机实习报告（2015～2016学年第一学期）姓名：xxxxx xx学号：xxxxxxxxxx专业：岩土工程指导教师：***联系电话：xxxxxxxxxxx实习成绩：xxxxxxxxx2015年12月10日目录一序言 (3)二正文 (3)题目3 (3)原理3 (3)结果3 (4)分析3 (5)题目4 (6)原理4 (6)结果4 (7)分析4 (7)题目5 (7)原理5 (7)结果5 (8)分析5 (9)三总结 (9)四附录 (9)附录1雅格比迭代法程序代码 (9)附录2高斯－赛德尔迭代法程序 (10)附录3求解题目3程序代码 (11)附录4 SOR法程序代码 (12)附录5求解题目4程序代码 (13)附录6 Ru n ge-Kutt a 4阶算法程序代码 (13)附录7求解题目5程序代码 (14)一 序言MATLAB 的M 语言，一种演算纸方式的编程语言。

通过这种语言，用户可以用类似于数学公式的方式来编写算法，大大降低了编程所需的难度并节省了时间，从而让用户把主要的精力集中在算法的构思而不是编程上。

为便于检验结果，本上机实习全部使用M 语言编程，然后用内置函数求解进行对比。

二 正文题目3用雅格比法与高斯－赛德尔迭代法解下列方程组Ax =b ，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。

(1) 12621-3100142, b 2, b -2003144345A -⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪=-== ⎪ ⎪ ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ (2) 1210.80.8350.810.8, b 2, b 00.80.811-10A ⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪=== ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭(3) 134, b 716A ⎛⎫⎛⎫== ⎪ ⎪-⎝⎭⎝⎭原理： 雅格比迭代法：)b x a x a x a x a (a 1x )b x a x a x a (a 1x )b x a x a x a (a 1x n )1k (1n 1nn )1k (33n )1k (22n )1k (11n nn)k (n2)1k (n n 2)1k (323)1k (12122)k (21)1k (nn 1)1k (313)1k (21211)k (1-++++-=-+++-=-+++-=------------Jacobi 迭代也可看成简单迭代的一种，故对简单迭代的所有性质也成立。

从上可知：如果矩阵A 的主对角元不为零，则其Jacobi 迭代是唯一的。

如用矩阵形式表示：则迭代矩阵：B=I －A D 1-其中：g= b ，D=diag(a 11,…，a nn )Jacobi 迭代收敛的充要条件是ρ（Ｉ－1D -Ａ）<１。

Gauss-Seidel 迭代法()(1)(1)(1)11221331111()()(1)(1)22112332222()()()()()112233111()1()1()k k k k n n k k k k n n k k k k k n n n n nn n n nnx a x a x a x b a x a x a x a x b a x a x a x a x a x b a -------=-+++-=-+++-=-++++-我们称它为方程组Ax=b 的Gauss-Seidel 迭代式，如写成矩阵形式为： x (k)= D -1 (L x (k)+Ux (k-1))+ D -1bx (k)= (D-L)-1U x (k-1)+ (D-L)-1b其中：L=－12121(1)1210000,00nn n n n nn a a aU a a a a --⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦D=diag(a 11,…，a nn )Gauss-Seidel 迭代法的迭代矩阵为(D －L 1)-U ，常数项为(D －L 1)-b ，收敛的充要条件是ρ（(D －L 1)-U)<1 结果3取()(1)10^(8)k k x x ∞-≤--(1)分析3GS 迭代收敛速度一般比Jacobi 迭代收敛速度快，右端项对迭代是否收敛没有影响，但有时对迭代次数会产生较大的影响。

题目4松弛因子对SOR 法收敛速度的影响。

用SOR 法求解方程组Ax =b ,其中41-3141-2...-2, b ....141-214-3A -⎛⎫⎛⎫⎪ ⎪- ⎪ ⎪ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎪ ⎪- ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭要求程序中不存系数矩阵A ,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代，记录近似解x (k)达到||x (k)-x (k-1)||<10-6时所用的迭代次数k ，观察松弛因子对收敛速度的影响，并观察当w ≤0或w ≥2会有什么影响？ 原理：逐次超松弛迭代法（SOR-迭代法）：选取矩阵A 的下三角矩阵分量并赋予参数w ，将之作为分裂矩阵M ，)(1wL D wM -=，其中，w>0，为可选择的松弛因子，又（1）公式构造一个迭代法，其迭代矩阵为A wL D w I B s 1)(---≡从而得到解b X A =*的逐次超松弛迭代法。

(0)(1)()k 0,1,2,....k k s XXB X f +⎧⎪⎨=+=⎪⎩ 其中：bwL D w f wU D w wL D B s 11)())1(()(---=+--=由此，解b X A =*的SOR-迭代法的计算公式为⎪⎩⎪⎨⎧--+==∑∑-=+=+++)(),....,(111)1()1()()1()0()0(2)0(1)0(i j n i j k j ij k j ij i ii k i k iTn X a X a b a w X X X X X X （2）观察（2）式，可得结论：（1）当w=1时，SOR-迭代法为J-迭代法。

（2）当w>1时，称为超松弛迭代法，当w<1时，称为低松弛迭代 结果4取()(1)10^(6)k k x x ∞-≤--时所用的迭代次数k 列表如下：分析4松弛因子的选取会对迭代次数及和是否收敛产生较大影响，松弛因子w 应该满足0<w<2，否则会出错或者出现发散的情况。

题目5用Ru n ge-Kutt a 4阶算法对初值问题y /=-20*y ，y (0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。

注：此方程的精确解为：y =e -20x 原理：一阶常微分方程初值问题的数值解法是近似计算中很重要的部分。

0(,)()dyf x y dxy x y⎧=⎪⎨⎪=⎩ （5.1） 常微分方程初值问题的数值解法是求方程（5.1）的解在点列1(0,1,)n n n x x h n -=+=上的近似值n y ，这里n h 是1n x -到n x 的步长，一般略去下标记为h 。

常微分方程初值问题的数值解法一般分为两大类：（1）单步法：这类方法在计算n y 时，只用到1n x +、n x 和n y ，即前一步的值。

因此，在有了初值以后就可以逐步往下计算。

典型方法如龙格–库塔()R K -方法。

（2）多步法：这类方法在计算1n y +时，除用到1n x +、n x 和n y 以外，还要用(1,2,,;0)n p y p k k -=>，即前面k 步的值。

典型方法如Adams 方法。

经典的R K -方法是一个四阶的方法，它的计算公式是：112341213243(22)6(,)(,)22(,)22(,)n n n n n nn n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （6.2） R K -方法的优点是：单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值f 。

结果5分析5步长的选取会直接影响能否求出准确结果以及所求结果的准且性。

三 总结实际工程中的数学问题往往比较复杂，需要借助于计算机这一工具进行计算。

由于计算机实质上只会做加减乘除运算，所以研究怎样通过计算机所能执行的基本计算，求得数学问题的有效数值解是数值分析的最终任务。

数值分析计算原理一般都不会太难，但是运算量很大，尤其是问题中存在大型矩阵时，情况更是如此。

同时，运算方法的选择，初始参数的选择都会对能否得到正确结果和得到正确结果的运算量产生较大影响。

四附录附录1雅格比迭代法程序代码function [x,k]=Jacobimethod(A,b,x0,N,emg)% A:线性方程组左端矩阵% b:线性方程组右端向量% x0:迭代初值% N:迭代次数上界，若迭代次数大于n，则迭代失败% emg:精度指标% k:迭代次数% x:用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1);x2=zeros(n,1);x1=x0; k=0;r=max(abs(b-A*x1));while r>emgfor i=1:nsum=0;for j=1:nif i~=jsum=sum+A(i,j)*x1(j);endendx2(i)=(b(i)-sum)/A(i,i);endr=max(abs(x2-x1));x1=x2;k=k+1;if k>Ndisp('迭代失败，返回')return;endendx=x1;附录2高斯－赛德尔迭代法程序function [x,k]=Gaussmethod(A,b,x0,N,emg)% A:线性方程组左端矩阵% b:线性方程组右端向量% x0:迭代初值% N:迭代次数上界，若迭代次数大于n，则迭代失败% emg:精度指标% k:迭代次数% x:用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1);x2=zeros(n,1);x1=x0;r=max(abs(b-A*x1));k=0;while r>emgfor i=1:nsum=0;for j=1:nif j>isum=sum+A(i,j)*x1(j);elseif j<isum=sum+A(i,j)*x2(j);endendx2(i)=(b(i)-sum)/A(i,i);endr=max(abs(x2-x1));x1=x2;k=k+1;if k>Ndisp('迭代失败，返回')return;endendx=x1附录 3 求解题目3程序代码>>A=[6,2,-1;1,4,-2;-3,1,4]; b=[-3;2;4];x0=[0;0;0];>> [x,k]=Jacobimethod(A,b,x0,100,10^(-8))x =-0.72730.80810.2525k =31>> A=[6,2,-1;1,4,-2;-3,1,4]; b=[100;-200;345];x0=[0;0;0];>> [x,k]=Jacobimethod(A,b,x0,100,10^(-8))36.3636-2.0707114.0404k =37>>A=[1,0.8,0.8;0.8,1,0.8;0.8,0.8,1]; b=[3;2;2];x0=[0;0;0]; >> [x,k]=Jacobimethod(A,b,x0,100,10^(-8))迭代失败，返回>>A=[1,0.8,0.8;0.8,1,0.8;0.8,0.8,1]; b=[5;0;-10];x0=[0;0;0]; >> [x,k]=Jacobimethod(A,b,x0,100,10^(-8))迭代失败，返回>>A=[1,3;-7,1]; b=[4;6];x0=[0;0];>> [x,k]=Jacobimethod(A,b,x0,100,10^(-8))迭代失败，返回>>A=[6,2,-1;1,4,-2;-3,1,4]; b=[-3;2;4];x0=[0;0;0];>> [x,k]=Gaussmethod(A,b,x0,100,10^(-8))x =-0.72730.80810.2525k =19>>A=[6,2,-1;1,4,-2;-3,1,4]; b=[100;-200;345];x0=[0;0;0]; >> [x,k]=Gaussmethod(A,b,x0,100,10^(-8))x =36.3636-2.0707114.0404k =24>>A=[1,0.8,0.8;0.8,1,0.8;0.8,0.8,1]; b=[3;2;2];x0=[0;0;0]; >> [x,k]=Gaussmethod(A,b,x0,100,10^(-8))x =4.2308-0.7692-0.769255>>A=[1,0.8,0.8;0.8,1,0.8;0.8,0.8,1]; b=[5;0;-10];x0=[0;0;0];>> [x,k]=Gaussmethod(A,b,x0,100,10^(-8))x =32.69237.6923-42.3077k =65>>A=[1,3;-7,1]; b=[4;6];x0=[0;0];>> [x,k]=Gaussmethod(A,b,x0,100,10^(-8))迭代失败，返回附录4 SOR法程序代码function [x,k]=SORmethod(A,b,x0,N,emg,w)% A:线性方程组左端矩阵% b:线性方程组右端向量% x0:迭代初值% N:迭代次数上界，若迭代次数大于n，则迭代失败% emg:精度指标% w:松弛因子% x:用迭代法求得的线性方程组的近似解n=length(A);x1=zeros(n,1);x2=zeros(n,1);x1=x0;r=max(abs(b-A*x1));k=0;while r>emgfor i=1:nsum=0;for j=1:nif j>=isum=sum+A(i,j)*x1(j);elseif j<isum=sum+A(i,j)*x2(j);endendx2(i)=x1(i)+w*(b(i)-sum)/A(i,i);endr=max(abs(x2-x1));x1=x2;k=k+1;if k>Ndisp('迭代失败，返回')return;endendx=x1;附录5求解题目4程序代码（只有n=4,w=1.1的程序，其他类似）A=[-4 1 0 0;1 -4 1 0;0 1 -4 1;0 0 1 -4];b=[-3;-2;-2;-3];x0=[0;0;0;0];>> [x,k]=SORmethod(A,b,x0,10000,10^(-6),1.1)x =1.00001.00001.00001.0000k =8附录6 Ru n ge-Kutt a 4阶算法程序代码function R=Rungkuta4(f, a, b, n, ya)% f:微分方程右端函数句柄% a,b:自变量取值区间的两个端点% n:区间等分的个数% ya:函数初值y(a)% R=[x',y']:自变量X 和解Y 所组成的矩阵h=(b-a)/n;x=zeros(n+1);y=zeros(1,n+1);x=a:h:b;y(1)=ya;for i=1:nk1=h*feval(f,x(i),y(i));k2=h*feval(f,x(i)+h/2,y(i)+k1/2);k3=h*feval(f,x(i)+h/2,y(i)+k2/2);k4=h*feval(f,x(i)+h,y(i)+k3);y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;endR=[x',y'];附录7求解题目5程序代码f=@(x,y)-20*y;R=Rungkuta4(f, 0, 1, 10, 1)R =0 1.00000.1000 0.33330.2000 0.11110.3000 0.03700.4000 0.01230.5000 0.00410.6000 0.00140.7000 0.00050.8000 0.00020.9000 0.00011.0000 0.0000f=@(x,y)-20*y;R=Rungkuta4(f, 0, 1, 5, 1)R =1.0e+003 *0 0.00100.0002 0.00500.0004 0.02500.0006 0.12500.0008 0.62500.0010 3.1250。