清华大学数值分析实验报告

数值分析实验报告

一、 实验3.1 题目：

考虑线性方程组b Ax =，n n R A ?∈，n R b ∈，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss 消去过程。

（1）取矩阵??????????

?

??

???=681681681

6 A ，??

???

??

?

????????=1415157 b ，则方程有解()T x 1,,1,1*?=。取10=n 计算矩阵的条件数。分别用顺序Gauss 消元、列主元Gauss 消元和完全选

主元Gauss 消元方法求解，结果如何？

（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。

1. 算法介绍

首先，分析各种算法消去过程的计算公式， 顺序高斯消去法：

第k 步消去中，设增广矩阵B 中的元素()

0k kk a ≠（若等于零则可以判定系数

矩阵为奇异矩阵，停止计算），则对k 行以下各行计算()

(),1,2,,k

ik

ik k kk

a l i k k n a ==++ ，

分别用ik l -乘以增广矩阵B 的第k 行并加到第1,2,,k k n ++ 行，则可将增广矩阵

B 中第k 列中()k kk a 以下的元素消为零；重复此方法，从第1步进行到第n-1步，

则可以得到最终的增广矩阵，即()()(),n n n B A b ??=?

?； 列主元高斯消去法：

第k 步消去中，在增广矩阵B 中的子方阵()()()()k k kk kn

k k nk

nn a a a a ????

????

?? 中，选取()k k i k a 使得

()(k)

max k k

i k ik k i n

a a ≤≤=，当k i k ≠时，对B 中第k 行与第k i 行交换，然后按照和顺序消去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即()

()()111,n n n B

A b ??=?

?； 完全主元高斯消去法：

第k 步消去中，在增广矩阵B 中对应的子方阵()()()()k k kk kn

k k nk

nn a a a a ????????

?? 中，

选取()k k k

i j a 使得()(k)

max k k k i j ij k i n

k j n

a a ≤≤≤≤=，若k i k ≠或k j k ≠，则对B 中第k 行与第k i 行、第k 列与第k

j 列交换，然后按照和顺序消去法相同的步骤进行即可。重复此方法，从第1步进行到第n-1步，就可以得到最终的增广矩阵，即()

()()222,n n n B

A b ??=?

?； 接下来，分析回代过程求解的公式，容易看出，对上述任一种消元法，均有以下计算公式：

()()()

()

1;

;1,2,,1

n

n

n n nn

n

n n k kj

j

j k k n kk

b x a b a

x x k n n a

=+=-=

=--∑

2. 实验程序的设计

一、输入实验要求及初始条件；

二、计算系数矩阵A 的条件数及方程组的理论解； 三、对各不同方法编程计算，并输出最终计算结果。

3. 计算结果及分析

（1）

先计算系数矩阵的条件数，结果如下，

12()2557.500000,()1727.556025,()2557.50000cond A cond A cond A ∞===

可知系数矩阵的条件数较大，故此问题属于病态问题，b 或A 的扰动都可能引起解的较大误差；

采用顺序高斯消去法，计算结果为：

最终解为x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028)T

使用无穷范数衡量误差，得到*X X ε∞=-=2.842170943040401e-14

，可以

发现，采用顺序高斯消元法求得的解与精确解之间误差较小。通过进一步观察，可以发现，按照顺序高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此顺序高斯消去法方式并没有对结果造成特别大的影响。 若采用列主元高斯消元法，则结果为：

最终解为x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000)T

同样使用无穷范数衡量误差，有*X X ε∞=-=0； 若使用完全主元高斯消元法，则结果为

最终解x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000)T

同样使用无穷范数衡量误差，有*X X ε∞=-=0； （2）

若每步都选取模最小或尽可能小的元素为主元，则计算结果为

最终解x=(1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000012 0.999999999999979 1.000000000000028)T

使用无穷范数衡量误差，有*X X ε∞=-为2.842170943040401e-14；而完全主元消去法的误差为*X X ε∞=-=0。

从（1）和（2）的实验结果可以发现，列主元消去法和完全主元消去法都得到了精确解，而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了一定的影响，但由于方程组的维度较低，并且元素之间相差不大，所以误差仍比较小。

为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下：

从上表可以发现，对这个方程组而言，顺序高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小，所以计算过程中的累积误差也较小，最终4种方法的输出结果均较为精确。

在这里，具体解释一下顺序法与模最小法的计算结果完全一致的原因。该矩阵在消元过程中，每次选取主元的一列只有两个非零元素，对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。在这样的情况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一致的。

理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特殊性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。

（3）

n=时，重复上述实验过程，各种方法的计算结果如下所示，在这里，仍20

采用无穷范数衡量绝对误差ε。

可以看出，此时列主元和完全主元的计算结果仍为精确值，而顺序高斯消去和模尽可能小方法仍然产生了一定的误差，并且两者的误差一致。与n=10时候的误差比相比，n=20时的误差增长了大约1000倍，这是由于计算过程中舍入误差的不断累积所致。所以，如果进一步增加矩阵的维数，应该可以看出更明显的现象。（4）

不同矩阵维度下的误差如下，在这里，为方便起见，选取2-条件数对不同维度的系数矩阵进行比较。

从上表可以看出，随着维度的增加，不同方法对计算误差的影响逐渐体现，并且

增长较快，这是由于舍入误差逐步累计而造成的。不过，方法二与方法三在维度小于40的情况下都得到了精确解，这两种方法的累计误差远比方法一和方法四慢；同样地，出于与前面相同的原因，方法一与方法四的计算结果保持一致，方法二与方法三的计算结果保持一致。

4. 结论

本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此，不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显，四种方法都足够精确。

对比四种方法，可以发现采用列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确。不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。

另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发现列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。

附录：程序代码

clear

clc;

format long;

%方法选择

n=input('矩阵A阶数:n=');

disp('选取求解方式');

disp('1 顺序Gauss消元法，2 列主元Gauss消元法，3 完全选主元Gauss消元法，4 模最小或近可能小的元素作为主元');

a=input('求解方式序号：');

%赋值A和b

A=zeros(n,n);

b=zeros(n,1);

for i=1:n

A(i,i)=6;

if i>1

A(i,i-1)=8;

end

if i

A(i,i+1)=1;

end

end

for i=1:n

for j=1:n

b(i)=b(i)+A(i,j);

end

end

disp('给定系数矩阵为：');

A

disp('右端向量为：');

b

%求条件数及理论解

disp('线性方程组的精确解：');

X=(A\b)'

fprintf('矩阵A的1-条件数: %f \n',cond(A,1)); fprintf('矩阵A的2-条件数: %f \n',cond(A)); fprintf('矩阵A的无穷-条件数: %f \n',cond(A,inf));

if a==1

A1=A;b1=b;

for k=1:n

if A1(k,k)==0

disp('主元为零，顺序Gauss消元法无法进行');

break

end

fprintf('第%d次消元所选取的主元：%g\n',k,A1(k,k))

%disp('此次消元后系数矩阵为：');

%A1

for p=k+1:n

l=A1(p,k)/A1(k,k);

A1(p,k:n)=A1(p,k:n)-l*A1(k,k:n);

b1(p)=b1(p)-l*b1(k);

end

end

x1(n)=b1(n)/A1(n,n);

for k=n-1:-1:1

for w=k+1:n

b1(k)=b1(k)-A1(k,w)*x1(w);

end

x1(k)=b1(k)/A1(k,k);

end

disp('顺序Gauss消元法解为:');

disp(x1);

disp('所求解与精确解之差的无穷-范数为');

norm(x1-X,inf)

end

if a==2

A2=A;b2=b;

for k=1:n

[max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k))));

if max_i~=k

A2_change=A2(k,:);

A2(k,:)=A2(max_i,:);

A2(max_i,:)=A2_change;

b2_change=b2(k);

b2(k)=b2(max_i);

b2(max_i)=b2_change;

end

if A2(k,k)==0

disp('主元为零，列主元Gauss消元法无法进行');

break

end

fprintf('第%d次消元所选取的主元：%g\n',k,A2(k,k))

%disp('此次消元后系数矩阵为：');

%A2

for p=k+1:n

l=A2(p,k)/A2(k,k);

A2(p,k:n)=A2(p,k:n)-l*A2(k,k:n);

b2(p)=b2(p)-l*b2(k);

end

end

x2(n)=b2(n)/A2(n,n);

for k=n-1:-1:1

for w=k+1:n

b2(k)=b2(k)-A2(k,w)*x2(w);

end

x2(k)=b2(k)/A2(k,k);

end

disp('列主元Gauss消元法解为:');

disp(x2);

disp('所求解与精确解之差的无穷-范数为');

norm(x2-X,inf)

end

%完全选主元Gauss消元法

if a==3

A3=A;b3=b;

for k=1:n

VV=eye(n);

[max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n)))));

if numel(max_i)==0

[max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n)))));

end

W=eye(n);

W(max_i(1)+k-1,max_i(1)+k-1)=0;

W(k,k)=0;

W(max_i(1)+k-1,k)=1;

W(k,max_i(1)+k-1)=1;

V=eye(n);

V(k,k)=0;

V(max_j(1)+k-1,max_j(1)+k-1)=0;

V(k,max_j(1)+k-1)=1;

V(max_j(1)+k-1,k)=1;

A3=W*A3*V;

b3=W*b3;

VV=VV*V;

if A3(k,k)==0

disp('主元为零，完全选主元Gauss消元法无法进行');

break

end

fprintf('第%d次消元所选取的主元：%g\n',k,A3(k,k))

%disp('此次消元后系数矩阵为：');

%A3

for p=k+1:n

l=A3(p,k)/A3(k,k);

A3(p,k:n)=A3(p,k:n)-l*A3(k,k:n);

b3(p)=b3(p)-l*b3(k);

end

end

x3(n)=b3(n)/A3(n,n);

for k=n-1:-1:1

for w=k+1:n

b3(k)=b3(k)-A3(k,w)*x3(w);

end

x3(k)=b3(k)/A3(k,k);

end

disp('完全选主元Gauss消元法解为:');

disp(x3);

disp('所求解与精确解之差的无穷-范数为');

norm(x3-X,inf)

end

%模最小或近可能小的元素作为主元

if a==4

A4=A;b4=b;

for k=1:n

AA=A4;

AA(AA==0)=NaN;

[min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k))));

if numel(min_i)==0

[min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n))));

end

W=eye(n);

W(min_i(1)+k-1,min_i(1)+k-1)=0;

W(k,k)=0;

W(min_i(1)+k-1,k)=1;

W(k,min_i(1)+k-1)=1;

A4=W*A4;

b4=W*b4;

if A4(k,k)==0

disp('主元为零，模最小Gauss消元法无法进行');

break

end

fprintf('第%d次消元所选取的主元：%g\n',k,A4(k,k))

%A4

for p=k+1:n

l=A4(p,k)/A4(k,k);

A4(p,k:n)=A4(p,k:n)-l*A4(k,k:n);

b4(p)=b4(p)-l*b4(k);

end

end

x4(n)=b4(n)/A4(n,n);

for k=n-1:-1:1

for w=k+1:n

b4(k)=b4(k)-A4(k,w)*x4(w); end

x4(k)=b4(k)/A4(k,k); end

disp('模最小Gauss 消元法解为:'); disp(x4);

disp('所求解与精确解之差的无穷-范数为'); norm(x4-X,inf) end

二、实验3.3 题目：

考虑方程组Hx b =的解，其中系数矩阵H 为Hilbert 矩阵：

()

,,1

,,,1,2,,n

1i j i j n n

H h h i j i j ?==

=+-

这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端的办法给出确定的问题。

（1）选择问题的维数为6，分别用Gauss 消去法（即LU 分解）、J 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何。

（2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何？计算的结果说明的什么？

（3）讨论病态问题求解的算法。

1. 算法设计

对任意线性方程组Ax b =，分析各种方法的计算公式如下， （1）Gauss 消去法：

首先对系数矩阵进行LU 分解，有A LU =，则原方程转化为LUx b =，令

Ux y =，则原方程可以分为两步回代求解：

Ux y Ly b =??

=?

具体方法这里不再赘述。 （2）J 迭代法：

首先分解A D L U =--，再构造迭代矩阵()()1k k X BX f +=+，其中

()11,B D L U f D b --=+=，进行迭代计算，直到误差满足要求。 （3）GS 迭代法：

首先分解A D L U =--，再构造迭代矩阵()()1k k X GX f +=+，其中

()()1

1

,G D L U f D L b --=-=-，进行迭代计算，直到误差满足要求。

（4）SOR 迭代法：

首先分解A D L U =--，再构造迭代矩阵()()1k k X L X f ω+=+，其中

()()()1

1

1,L D L D U f D L b ωωωωωω--=--+=-????，进行迭代计算，直到误差满

足要求。

2. 实验过程

一、根据维度n 确定矩阵H 的各个元素和b 的各个分量值；

二、选择计算方法（ Gauss 消去法，J 迭代法，GS 迭代法，SOR 迭代法），对迭代法设定初值，此外SOR 方法还需要设定松弛因子；

三、进行计算，直至满足误差要求（对迭代法，设定相邻两次迭代结果之差的无穷范数小于0.0001；对SOR 方法，设定为输出迭代100次之后的结果及误差值），输出实验结果。

3. 计算结果及分析

（1）6n =时，问题可以具体定义为

计算结果如下，

Gauss 消去法

第1次消元所选取的主元是：1 第2次消元所选取的主元是：0.0833333 第3次消元所选取的主元是：0.00555556 第4次消元所选取的主元是：0.000357143 第5次消元所选取的主元是：2.26757e-05 第6次消元所选取的主元是：1.43155e-06

解得X=(0.999999999999228 1.000000000021937 0.999999999851792 1.000000000385369 0.999999999574584 1.000000000167680)T

使用无穷范数衡量误差，可得*X X ε∞=-=4.254160357319847e-10；

J 迭代法

设定迭代初值为零，计算得到J 法的迭代矩阵B 的谱半径为4.30853＞1，所以J 法不收敛；

GS 迭代法

设定迭代初值为零，计算得到GS 法的迭代矩阵G 的谱半径为：0.999998＜1，故GS 法收敛，经过541次迭代计算后，结果为X=(1.001178105812706 0.9991440826518600.9689290939849021.0470455699891621.0273231583702810.

954352032784608)T

使用无穷范数衡量误差，有*X X ε∞=-=0.047045569989162；

SOR 迭代法

设定迭代初值为零向量，并设定0.5ω=，计算得到SOR 法迭代矩阵谱半径为0.999999433815223，经过100次迭代后的计算结果为

X=(1.0033806141450780.9624202974584231.0318570231345591.0618149012898811.0140378158271640.917673642493527)T ；

使用无穷范数衡量误差，有*X X ε∞=-=0.082326357506473； 对SOR 方法，ω可变，改变ω值，计算结果可以列表如下

可以发现，松弛因子的取值对迭代速度造成了不同的影响，上述四种方法中，松弛因子=0.5时，收敛相对较快。

综上，四种算法的结果列表如下：

计算可得，矩阵H 的条件数为71.5010?>>1，所以这是一个病态问题。由上表可以看出，四种方法的求解都存在一定的误差。下面分析误差的来源：

LU 分解方法的误差存在主要是由于Hilbert 矩阵各元素由分数形式转换为小数形式时，不能除尽情况下会出现舍入误差，在进行LU 分解时也存在这个问题，所以最后得到的结果不是方程的精确解，但结果显示该方法的误差非常小；

Jacobi 迭代矩阵的谱半径为4.30853，故此迭代法不收敛；

GS 迭代法在迭代次数为541

次时得到了方程的近似解，其误差约为0.05 ，比较大。GS 迭代矩阵的谱半径为0.999998，很接近1，所以GS 迭代法收敛速度较慢；

SOR 迭代法在迭代次数为100次时误差约为0.08，误差较大。SOR 迭代矩阵的谱半径为0.999999，也很接近1，所以0.5ω=时SOR 迭代法收敛速度不是很快，但是相比于GS 法，在迭代速度方面已经有了明显的提高；另外，对不同的ω，SOR 方法的迭代速度会相应有变化，如果选用最佳松弛因子，可以实现更快的收敛； （2）

考虑不同维度的情况，8n =时，

11n =时，

15n =时