西南交大数值分析上机实习报告(共15页).docx

《西南交大数值分析上机实习报告(共15页).docx》由会员分享，可在线阅读，更多相关《西南交大数值分析上机实习报告(共15页).docx（15页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上数值分析上机实习报告（20152016学年第一学期）姓 名： xxxxx xx 学 号： xxxxxxxxxx 专 业： 岩土工程 指导教师： 徐跃良 联系电话： xxxxxxxxxxx 实习成绩： xxxxxxxxx2015年 12月 10日目录一 序言3二 正文3题目33原理33结果34分析35题目46原理46结果47分析47题目57原理57结果58分析59三 总结9四 附录9附录1雅格比迭代法程序代码9附录2高斯赛德尔迭代法程序10附录3求解题目3程序代码11附录4 SOR法程序代码12附录5求解题目4程序代码13附录6 Runge-Kutta 4阶算法程序代码

2、13附录7求解题目5程序代码14一 序言 MATLAB的M语言，一种演算纸方式的编程语言。通过这种语言，用户可以用类似于数学公式的方式来编写算法，大大降低了编程所需的难度并节省了时间，从而让用户把主要的精力集中在算法的构思而不是编程上。为便于检验结果，本上机实习全部使用M语言编程，然后用内置函数求解进行对比。二 正文题目3用雅格比法与高斯赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。(1) (2) (3) 原理： 雅格比迭代法： Jacobi迭代也可看成简单迭代的一种，故对简单迭代的所有性质也成立。从上可知：如果矩阵

3、A的主对角元不为零，则其Jacobi迭代是唯一的。如用矩阵形式表示：则迭代矩阵：B=I 其中：g= b，D=diag(a11,，ann)Jacobi迭代收敛的充要条件是r（）。Gauss-Seidel迭代法 我们称它为方程组Ax=b的Gauss-Seidel迭代式，如写成矩阵形式为： x(k)= D-1 (L x(k)+Ux(k-1)+ D-1b x(k)= (D-L)-1U x(k-1)+ (D-L)-1b其中：L= D=diag(a11,，ann)Gauss-Seidel迭代法的迭代矩阵为(DLU，常数项为(DLb，收敛的充要条件是r（(DLU)1结果3取(1)bJacobi迭代次数kJa

4、cobi解xGS迭代次数kGS解xb131（-0.7273，0.8081，0.2525）T19(-0.7273,0.8081,0.2525)Tb237(36.3636,-2.0707,114.0404)T24(36.3636,-2.0707,114.0404)T(2)bJacobi迭代次数kJacobi解xGS迭代次数kGS解xb1发散无法求出55(4.2308,-0.7692,-0.7692)Tb2发散无法求出65(32.6923,7.6923,-42.3077)T(3)bJacobi迭代次数kJacobi解xGS迭代次数kGS解xb1发散无法求出发散无法求出分析3GS迭代收敛速度一般比Ja

5、cobi迭代收敛速度快，右端项对迭代是否收敛没有影响，但有时对迭代次数会产生较大的影响。题目4松弛因子对SOR法收敛速度的影响。用SOR法求解方程组Ax=b,其中要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, .,1.9进行迭代，记录近似解x(k)达到|x(k)-x(k-1)|0，为可选择的松弛因子，又（1）公式构造一个迭代法，其迭代矩阵为从而得到解的逐次超松弛迭代法。其中：由此，解的SOR-迭代法的计算公式为 （2）观察（2）式，可得结论：（1） 当w=1时，SOR-迭代法为J-迭代法。（2）当w1时，称为超松弛迭代法，当w1时，称为低松弛迭代结果4取时所用的迭代次数k列

6、表如下：w n 1.11.21.31.41.51.61.71.81.9-101234811141822314266138出错出错10发散出错6912151824304367142出错出错12发散出错81113161922334267140出错出错12发散出错101215172024304568140出错出错13发散出错分析4松弛因子的选取会对迭代次数及和是否收敛产生较大影响，松弛因子w应该满足0wemg for i=1:n sum=0; for j=1:n if i=j sum=sum+A(i,j)*x1(j); end end x2(i)=(b(i)-sum)/A(i,i); end r=m

7、ax(abs(x2-x1); x1=x2; k=k+1; if kN disp(迭代失败，返回) return; end endx=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 remg fo

8、r i=1:n sum=0; for j=1:n if ji sum=sum+A(i,j)*x1(j); elseif jN disp(迭代失败，返回) return; end endx=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.7273 0.8081 0.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

9、)x = 36.3636 -2.0707 114.0404k =37A=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,

10、-2;-3,1,4; b=-3;2;4;x0=0;0;0; x,k=Gaussmethod(A,b,x0,100,10(-8)x = -0.7273 0.8081 0.2525k =19A=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.0707 114.0404k =24A=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.23

11、08 -0.7692 -0.7692k =55A=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.6923 7.6923 -42.3077k =65A=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:迭代次数

12、上界，若迭代次数大于 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 remg for i=1:n sum=0; for j=1:n if j=i sum=sum+A(i,j)*x1(j); elseif jN disp(迭代失败，返回) return; end endx=x1; 附录5求解题目4程序代码（只有n=4,w=1.1的程序，其他类似）A=-4 1 0 0;1 -4 1 0;0 1 -4 1;0

13、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.0000 1.0000 1.0000 1.0000k = 8附录6 Runge-Kutta 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)=

14、ya; for i=1:n k1=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.0000 0.1000 0.3333 0.2000 0.1111 0.3000 0.0370 0.4000 0.0123 0.5000 0.0041 0.6000 0.0014 0.7000 0.0005 0.8000 0.0002 0.9000 0.00011.0000 0.0000f=(x,y)-20*y;R=Rungkuta4(f, 0, 1, 5, 1)R = 1.0e+003 * 0 0.0010 0.0002 0.0050 0.0004 0.0250 0.0006 0.1250 0.0008 0.6250 0.0010 3.1250专心-专注-专业