数值分析大作业

数值分析大作业

学 号：*********

专 业：机械工程

学生姓名：***

2014年10月

摘要：在自然科学与工程技术中，很多问题的解决常常归结为求解线性方程组Ax=b。随着计算机的发展，利用计算机这个强有力的计算工具去求解线性方程组是一个非常实用的问题。在求解大型线性方程组时，直接法在多次消元，回代的过程中，四则运算的误差累计与传播无法控制，致使计算结果的精度就无法保证，特别是求解大型稀松矩阵时，还要对系数矩阵进行分解。而迭代法相对于直接法而言，具有保持迭代矩阵不变的特点，计算程序一般也比较简单，且对于许多问题收敛速度比较快。比较常用的迭代法有雅克比迭代法、高斯一塞德尔迭代法和逐次超松弛迭代法等，本次研究目的是通过求解一个线性方程组来比较它们的迭代效果，验证一些已有的结论。

1．数学原理

1.1雅可比迭代法

将线性方程组的系数矩A?(aij)?Rn?n分解为A=D+L+U，其中D是由A的主对角元素构成的对角矩阵，L是由A的严格下三角部分构成的严格下三角矩阵， U是由A的严格上三角部分构成的严格上三角矩阵，即

?a11?D??

??????L??

????

0a21a31?an1

a22

???, ??

?ann?

??0?????,U??????

?0???

a12

a13a23?

0??

a1n?a2n????. ?an?1,n?0??

a32?an2

0?

?an,n?1

若系数矩阵A的对角元素aii?0(i?1,2,?,n)，则矩阵D非奇异，取M=D，N=-（L+U)，则x??D?1(L?U)x?D?1b?GJx?gJ， 因而，构造的迭代法为：

x(k?1)?GJx(k)?gJ,GJ??D(L?U),gJ?Db.

1.2高斯-赛得尔迭代法

?1

?1

将线性方程组的系数矩A?(aij)?Rn?n分解为A=D+L+U。若系数矩阵A的对角元素不等于0，则矩阵D非奇异，取M=L+D，N=-U，则

x???L?D?Ux??L?D?b?GGx?gG

?1

?1

因而，构造的迭代法为：

x(k?1)?GGx(k)?gG,

GG???L?D?U,gG??L?D?b.?1?1

1.3逐次超松弛迭代法

线性方程组的系数矩A?(aij)?Rn?n分解为A=D+L+U。 取M?1

??D??L?,N???1???D??U?,则 ?

x??D??L???1???D??U?x???D??L?b?GSx?gS ?1?11因而，构造的迭代法为：

x(k?1)?GSx(k)?gS,

GG??D??L???1???D??U?,gS???D??L?b.?1?1

2.程序设计

2.1雅可比迭代法求解

雅可比迭代法MATLAB程序如下：

%majacobi.m

function x=majacobi(A,b,x0,ep,N)

%用途：用Jacobi迭代法解线性方程组Ax=b

%格式;x=majacobi(A,b,x0,ep,N) A为系数矩阵，b为右端向量， %x0为初始向量（默认零向量），ep为精度（默认1e-6), %N为最大迭代次数（默认500次），x返回近似解向量 n=length(b);

if nargin<5,N=500;end

if nargin<4,ep=1e-6;end

if nargin<3,x0=zeros(n,1);end

x=zeros(n,1);k=0;

while k<N

for i=1:n

x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); end

if norm(x-x0,inf)<ep,break; end

x0=x ;k=k+1;

end

if k==N,Warning('100');end

disp(['k=',num2str(k)])

在MATLAB命令窗口执行创建的m文件

>> edit majacobi

>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]

A =

10.9000 1.2000 2.1000 0.9000

1.2000 11.2000 1.5000 2.5000

2.1000 1.5000 9.8000 1.3000

0.9000 2.5000 1.3000 12.3000

>> b=[-7.0,5.3,10.3,24.6]';

>> x=majacobi(A,b)

k=17

x =

-0.9986

0.0071

1.0032

1.9656

2.2高斯-赛得尔迭代法求解

高斯-赛得尔迭代法MATLAB程序如下：

%maseidel.m

function x=maseidel (A,b,x0,ep,N)

%用途：用Gauss-Seidel迭代法解线性方程组Ax=b

%格式：x=maseidel （A,b,x0,ep,N) A为系数矩阵，b为右端向量， %x0为初始向量（默认零向量），ep为精度（默认1e-6),

%N为最大迭代次数（默认500次），x返回近似解向量

n=length(b);

if nargin<5,N=500;end

if nargin<4,ep=1e-6;end

if nargin<3,x0=zeros(n,1);end

x=zeros(n,1);k=0;

while k<N

for i=1:n

if i==1

x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);

else if i==n

x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);

else

x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i); end

end

end

if norm(x-x0,inf)<ep,break;end

x0=x;k=k+1;

end

if k==N,Warning('100');end

disp(['k=',num2str(k)])

在MATLAB命令窗口执行创建的m文件：

>> edit maseidel

>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]

A =

10.9000 1.2000 2.1000 0.9000

1.2000 11.2000 1.5000 2.5000

2.1000 1.5000 9.8000 1.3000

0.9000 2.5000 1.3000 12.3000

>> b=[-7.0,5.3,10.3,24.6]';

>> x=maseidel(A,b)

k=7

x =

-0.9986

0.0071

1.0032

1.9656

2.3逐次超松弛迭代法求解

逐次超松弛迭代法MATLAB程序如下：

%masor.m

function x=masor(A,b,omega,x0,ep,N)

%用途：用SOR迭代法解线性方程组Ax=b

%格式：x=maseidel （A,b,x0,ep,N) A为系数矩阵，b为右端向量，

%omega为松弛因子（默认1.5),x0为初始向量（默认零向量），ep为精度（默认1e-6),

%N为最大迭代次数（默认500次），x返回近似解向量

n=length(b);

if nargin<6,N=500;end

if nargin<5,ep=1e-6;end

if nargin<4,x0=zeros(n,1);end

if nargin<3,omega=1.5;end

x=zeros(n,1);k=0;

while k<N

for i=1:n

if i==1

x1(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);

else if i==n

x1(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);