数值分析上机实验——解线性方程组

实验报告

哈尔滨工程大学教务处制

实验四 解线性方程组

一．解线性方程组的基本思想 1．直接三角分解法：

将系数矩阵A 转变成等价两个矩阵L 和U 的乘积 ，其中L 和U 分别是下三角和上三角矩阵。当A 的所有顺序主子式都不为0时，矩阵A 可以分解为A=LU ，且分解唯一。其中L 是单位下三角矩阵，U 是上三角矩阵。 2．平方根法：

如果矩阵A 为n 阶对称正定矩阵，则存在一个对角元素为正数的下三角实矩阵L ，使得：A=LL^T 。当限定L 的对角元素为正时，这种分解是唯一的，称为平方根法（Cholesky ）分解。 3．追赶法：

设系数矩阵为三对角矩阵

1122233111000000000

000000

n n n n

n b c a b c a b A a b c a b ---⎛⎫ ⎪ ⎪ ⎪=

⎪ ⎪

⎪

⎪ ⎪⎝

⎭

则方程组Ax=f 称为三对角方程组。

设矩阵A 非奇异，A 有Crout 分解A=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记

1122

233

1

1000010

000

0001000

000100,00000000

00

0001n n n

n b L U γαβγββγβ--⎛⎫⎛⎫ ⎪

⎪ ⎪ ⎪ ⎪ ⎪∂==

⎪

⎪ ⎪

⎪ ⎪ ⎪

⎪ ⎪

⎪

⎪

∂⎝

⎭

⎝

⎭ 可先依次求出L ，U 中的元素后，令Ux=y ，先求解下三角方程组Ly=f 得出y ，再求解上三

角方程组Ux=y 。

4．雅克比迭代法：

首先将方程组中的系数矩阵A 分解成三部分，即：A = L+D+U ，如图1所示，其中D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。

之后确定迭代格式，X )1(+k = BX )(k +f ，如图2所示，其中B 称为迭代矩阵，雅克比迭代法中一般记为J 。（k = 0,1，......）再选取初始迭代向量X )0(，开始逐次迭代。

5．超松弛迭代法（SOR ）

它是在GS 法基础上为提高收敛速度，采用加权平均而得到的新算法。 选取分裂矩阵M 为带参数的下三角矩阵

M ＝

ω

1

（D －L ω）， 其中ω＞0 为可选择的松弛因子，一般当1<ω<2时称为超松弛。 二.实验题目及实验目的

1．(第五章习题8)用直接三角分解（杜利特尔（Doolittle ）分解）求线性方程组

141x +251x +36

1

x = 9， 131

x +241x +351x = 8，

12

1

x + 2x +32x = 8 的解。

2．(第五章习题9)用追赶法解三对角方程组Ax=b ，其中

A=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------2100012100012100012100012，b=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛00001. 3．（第五章习题10）用改进的平方根法解线性方程组

⎪⎪⎪⎭⎫ ⎝⎛---131321112⎪⎪⎪⎭⎫ ⎝⎛321x x x = ⎪⎪⎪⎭

⎫ ⎝⎛654 4．（第六章习题7）用SOR 方法解线性方程组（分别取松弛因子ω=1.03，ω=1，

ω=1.1）

41x - 2x = 1， -1x +42x - 3x = 4，

-2x +43x = -3.

精确解x *=（

21，1，-2

1

）T .要求当)(*k x x -∞<5×106-时迭代终止，并且对每一个ω值确定迭代次数.

5.（第六章习题8）用SOR 方法解线性方程组（取ω=0.9） 51x -22x + 3x = -12， -1x +42x - 23x = 20， 21x -32x +103x = 3.

要求当)

()1(k k x x -+∞

<104-时迭代终止.

6．（第六章习题9）设有线性方程组Ax=b ，其中A 为对称正定阵，迭代公式

)()1(k k x x =++ω(b- A )(k x )，k=0，1，2…，

试证明当0<ω<

β

2

时上述迭代法收敛（其中0<α≤λ(A)≤β）. 7．（第六章计算实习题1）给出线性方程组H n x=b ，其中系数矩阵H n 为希尔伯特矩阵：

H n x=(h ij )∈R n n ⨯, h ij =

1

1

-+j i ，i ，j=1，2，…，n.

假设x *=（1，1，…，1）T ∈R n ，b= H n x *.若取n=6，8，10，分别雅克比迭代法及SOR 迭代（ω=1，1.25，1.5）求解.比较计算结果. 三．实验手段：

指操作环境和平台:win7系统下MATLAB R2009a

程序语言：一种类似C 语言的程序语言，但比C 语言要宽松得多，非常方便。 四.程序

1.

①直接三角分解（文件ZJsanjiao.m ） function x=ZJsanjiao(A,b) [m,n]=size(A); [l u]=lu(A); s=inv(l)*[A,b]; x=ones(m,1);

for i=m:-1:1

h=s(i,m+1);

for j=m:-1:1;

if j~=i

h=h-x(j)*s(i,j);

end

end

x(i)=h/s(i,i);

end

②控制台输入代码：

>> A=[1/4,1/5,1/6;1/3,1/4,1/5;1/2,1,2]; >> b=[9;8;8];

>> x=ZJsanjiao(A,b)

2.

①追赶法（文件ZG_SDJ.m）

function x=ZG_SDJ(a,b,c,f)

%aÊǶԽÇÏßÔªËØ

%bÊǶԽÇÏßÉÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö

%cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö

%fÊdz£ÊýÏîb

N=length(a);

b=[b,0];

c=[0,c];

a1=zeros(N,1);

b1=zeros(N,1);

y=zeros(N,1);

x=zeros(N,1);

a1(1)=a(1);

b1(1)=b(1)/a1(1);

y(1)=f(1)/a1(1);

for j1=2:N

a1(j1)=a(j1)-c(j1)*b1(j1-1);

b1(j1)=b(j1)/a1(j1);

temp1=f(j1)-c(j1)*y(j1-1);

y(j1)=temp1/a1(j1);

end

j1=N;

x(j1)=y(j1);

for j1=N-1:-1:1

x(j1)=y(j1)-b1(j1)*x(j1+1);