数值分析实验作业-gauss消去法的数值稳定性分析

实验3.1 Gauss 消去法的数值稳定性试验

实验目的:

观察和理解Gauss 消元过程中出现小主元(即)

(k kk

a 很小)时引起的方程

组解的数值不稳定性。 实验内容：

求解方程组b Ax =，其中

（1）⎥⎥⎥

⎥⎥⎦

⎤⎢⎢

⎢⎢

⎢

⎣

⎡⨯=11212592.1121

-130.6-291.51314.59103.015-1A ，⎥

⎥⎥⎥⎦⎤

⎢⎢⎢

⎢⎣⎡=2178.4617.591b ； （2）

⎥⎥⎥⎥⎦⎤⎢

⎢⎢

⎢⎣⎡----=2010

151526990999999999.2310710

2A ,

⎥⎥⎥⎥⎦⎤

⎢⎢⎢⎢⎣⎡=15

019000000000.582b . 实验要求：

（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。

（2） 用Gauss 列主元消去法求得L 和U 及解向量4

21,R x x ∈.

（3） 用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.

（4） 观察小主元并分析其对计算结果的影响. 程序如下:计算矩阵条件数及Gauss 列主元消去法: format longeng

A1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];

b1=[59.17;46.78;1;2]; n=4;

k2=cond(A1) %k2为矩阵的条件数；

for k=1:n-1

a=max(abs(A1(k:n,k))); [p,k]=find(A1==a); B=A1(k,:);c=b1(k);

A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)~=0

A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);

A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else

break end end

L1=tril(A1,0); for i=1:n

L1(i,i)=1; end L=L1

U=triu(A1,0) for j=1:n-1

b1(j)=b1(j)/L(j,j);

b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); end

b1(n)=b1(n)/L(n,n); for j=n:-1:2

b1(j)=b1(j)/U(j,j);

b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); end

b1(1)=b1(1)/U(1,1); x1=b1

运行结果如下： K2=68.43；

⎥

⎥

⎥⎥⎦⎤

⎢

⎢⎢

⎢⎣⎡

---⨯=-14929.00202.00893

.0011755.04724.00011079.2600

118L ⎥⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592

.11U 1x =[18.9882;3.3378;-34.747;-33.9865]

不选主元的Gauss 消去法程序： clear

format longeng

A1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;

for k=1:n-1

A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);

A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); end

L1=tril(A1,0); for i=1:n

L1(i,i)=1; end L=L1

U=triu(A1,0) for j=1:n-1

b1(j)=b1(j)/L(j,j);

b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); end

b1(n)=b1(n)/L(n,n); for j=n:-1:2

b1(j)=b1(j)/U(j,j);

b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); end

b1(1)=b1(1)/U(1,1); x1=b1

程序运行结果如下：

⎥

⎥

⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=

10189.010

333.3011168.21033.3700110637.17000

1~151515L ⎥⎥

⎥⎥

⎥

⎦

⎤⎢⎢

⎢⎢

⎢

⎣

⎡--⨯-⨯-⨯-⨯=

-5.00

08

16

010637.171091.5210043.101

314.59103.0~1515

18

15U ]0;0;0005.1;6848.23[~1=x

同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;