西南交大数值分析上机实习报告

数值分析上机实习报告

（2015～2016学年第一学期）

姓 名：学 号：专 业：指导教师：联系电话：实习成绩：

xxxxx xx xxxxxxxxxx 岩土工程 *** xxxxxxxxxxx

xxxxxxxxx

2015年 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

1

附录3求解题目3程序代码 ............................................................... 11 附录4 SOR法程序代码 ....................................................................... 12 附录5求解题目4程序代码 ............................................................... 13 附录6 Runge-Kutta 4阶算法程序代码 .............................................. 13 附录7求解题目5程序代码 ............................................................... 14

2

一 序言

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

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

题目3用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究

其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。

621-3100142, b2, b-200(1) A12 314434510.80.8350.810.8, b2, b0(2) A12 0.80.811-10(3) A, b 716原理： 雅格比迭代法：

(k)x11341a111a221annk1)(a12x(2(k1)k1)a13x3a1nx(b1)n

k)x(2(k1)(a21x1(k1)a23x3a2nx(nk1)b2)

k)x(n(k1)k1)(k1)k1)(an1x1an2x(an3x3ann1x(2n1bn)3

Jacobi迭代也可看成简单迭代的一种，故对简单迭代的所有性质也成立。从上可知：如果矩阵A的主对角元不为零，则其Jacobi迭代是唯一的。如用矩阵形式表示：则迭代矩阵：B=I－D1A

其中：g= b，D=diag(a11,…，ann)

Jacobi迭代收敛的充要条件是（Ｉ－D1Ａ）<１。 Gauss-Seidel迭代法

x1(k)1(k1)(a12x2a111(a21x1(k)a22(k1)a13x3(k1)a1nxnb1)

(k)x2(k1)a23x3(k1)a2nxnb2)

(k)xn1(k)(k)(an1x1(k)an2x2an3x3ann(k)ann1xn1bn)我们称它为方程组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

0,U0

a120

a1n

0a其中：L=－21an10an2ann1 a(n1)n0 D=diag(a11,…，ann)

Gauss-Seidel迭代法的迭代矩阵为(D－L)1U，常数项为(D－L)1b，收敛的充要条件是（(D－L)1U)<1 结果3

4

取(1)

b b1 b2 Jacobi迭代次数k 31 37 x(k)x(k1)10^(8)

Jacobi解x （-0.7273，0.8081，0.2525）T (36.3636,-2.0707,114.0404)T (2) GS迭代次数k 19 24 GS解x (-0.7273,0.8081,0.2525)T (36.3636,-2.0707,114.0404)T b b1 b2 Jacobi迭代次数k 发散 发散 Jacobi解x 无法求出 无法求出 GS迭代次数k 55 65 (3) GS解x (4.2308,-0.7692,-0.7692)T (32.6923,7.6923,-42.3077)T b b1 Jacobi迭代次数k 发散 Jacobi解x 无法求出 GS迭代次数k 发散 GS解x 无法求出 分析3

GS迭代收敛速度一般比Jacobi迭代收敛速度快，右端项对迭代是否收敛没有影响，但有时对迭代次数会产生较大的影响。 题目4松弛因子对SOR法收敛速度的影响。

用SOR法求解方程组Ax=b,其中

41-3141-2-2...A, b

....-2141-314要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代，记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当

5

w0或w2会有什么影响？ 原理：

逐次超松弛迭代法（SOR-迭代法）：

选取矩阵A的下三角矩阵分量并赋予参数w，将之作为矩阵M，

M1(DwL)，其中，w>0，为可选择的松弛因子，又（1）公式构造一个迭w代法，其迭代矩阵为BsIw(DwL)1A从而得到解A*Xb的逐次超松弛迭代法。

(0)X (k1)(k)BsXf k0,1,2,....X其中：

Bs(DwL)1((1w)DwU)fw(DwL)b1

由此，解A*Xb的SOR-迭代法的计算公式为

(0)(0)TX(0)(X1(0),X2,....Xn)i1n w(k1)(k)(k1)(k1)Xi(biaijXjaijXj)Xiaiij1ji1 （2）

观察（2）式，可得结论：

（1）当w=1时，SOR-迭代法为J-迭代法。

（2）当w>1时，称为超松弛迭代法，当w<1时，称为低松弛迭代 结果4

取w n 4 8 11 14 18 22 31 42 66 138 出错 出10 发出错 散 错 x(k)x(k1)10^(6)时所用的迭代次数

1.5 1.6 1.7 1.8 1.9 -1 k列表如下：

0 1 2 3 1.1 1.2 1.3 1.4 6

6 9 12 15 18 24 30 43 67 142 出错 出错 出错 出12 错 出12 错 出13 错 发散 发散 发散 出错 出错 出错 8 11 13 16 19 22 33 42 67 140 10 12 15 17 20 24 30 45 68 140 分析4

松弛因子的选取会对迭代次数及和是否收敛产生较大影响，松弛因子w应该满足0步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。 注：此方程的精确解为：y=e-20x 原理：

一阶常微分方程初值问题的数值解法是近似计算中很重要的部分。

dyf(x,y) （5.1） dxy(x)y00 常微分方程初值问题的数值解法是求方程（5.1）的解在点列

xnxn1hn(n0,1,)上的近似值yn，这里hn是xn1到xn的步长，一般略去下标记为h。

常微分方程初值问题的数值解法一般分为两大类：

（1）单步法：这类方法在计算yn时，只用到xn1、xn和yn，即前一步的值。因此，在有了初值以后就可以逐步往下计算。典型方法如龙格–库塔(RK)方法。

（2）多步法：这类方法在计算yn1时，除用到xn1、xn和yn以外，还要用ynp(p1,2,,k;k0)，即前面k步的值。典型方法如Adams方法。

7

经典的RK方法是一个四阶的方法，它的计算公式是：

hyy(K12K22K3K4)nn16K1f(xn,yn)hh （6.2） Kf(x,yK1)2nn22hhKf(x,yK2)nn322K4f(xnh,ynhK3)单步法、精度高，计算过程便于改变步长，RK方法的优点是：

缺点是计算量较大，每前进一步需要计算四次函数值f。 结果5

x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 h=0.1对应y 0.3333 0.1111 0.0370 0.0123 0.0041 0.0014 0.0005 0.0002 0.0001 0.0000 h=0.2对应y 5 25 125 625 3125 真值y 0.1353 0.0183 0.0025 3.3546e-004 4.5400e-005 6.1442e-006 8.3153e-007 1.1254e-007 1.5230e-008 2.0612e-009 分析5

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

实际工程中的数学问题往往比较复杂，需要借助于计算机这一工具进行计算。由于计算机实质上只会做加减乘除运算，所以研究怎样通过计算机所能执行的基本计算，求得数学问题的有效数值解是数值分析的最终任务。数值分析计算原理一般都不会太难，但是运算量很大，尤其是问题中存在大型矩阵时，情况更是如此。同时，运算方法

8

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

附录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>emg

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=max(abs(x2-x1)); x1=x2; k=k+1; if k>N

disp('迭代失败，返回') return; end end x=x1;

附录2高斯－赛德尔迭代法程序

function [x,k]=Gaussmethod(A,b,x0,N,emg) % A:线性方程组左端矩阵 % b:线性方程组右端向量

9

% 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>emg

for i=1:n sum=0; for j=1:n if j>i

sum=sum+A(i,j)*x1(j); elseif j sum=sum+A(i,j)*x2(j); end end x2(i)=(b(i)-sum)/A(i,i); end r=max(abs(x2-x1)); x1=x2; k=k+1;

if k>N

disp('迭代失败，返回') return; end end x=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.2525 k =

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

10

x = 36.3636 -2.0707 114.0404 k =

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.7273 0.8081 0.2525 k =

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.0707 114.0404 k =

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.7692

11

k =

55

>>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.6923 7.6923 -42.3077 k =

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>emg

for i=1:n sum=0; for j=1:n if j>=i

sum=sum+A(i,j)*x1(j); elseif jx2(i)=x1(i)+w*(b(i)-sum)/A(i,i); end r=max(abs(x2-x1)); x1=x2;

12

k=k+1; if k>N

disp('迭代失败，返回') return; end end

x=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.0000 1.0000 1.0000 1.0000 k = 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)=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; end R=[x',y'];

13

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.0001

1.0000 0.0000

f=@(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

14

附录