列主元高斯消去法的实现

《数值分析课程设计》

报告

专业：

学号：

学生姓名：

指导教师：

一、题目

列主元guess消去法求方程的解

二、理论

列主元高斯消去法是在高斯消去法的基础上而得到的一种比较快速合理的解线性方程组的方法。它的基本思想是每次在所在列对角线及以下元素中选择绝对值最大的元素作为主元进行消元计算。使用列主元消去法相对于高斯消去法更能减少舍入误差的影响。

三、方法、算法与程序设计

求解Ax=b

第一步：写出增广矩阵[A| b];

第二步：判断增广矩阵的秩r[A|b]与A的秩r[A]的关系：

若r[A|b]= r[A]，线性方程组有唯一解；

若r[A|b]>r[A]，线性方程组没有解；

若r[A|b]

第三步：若|A|≠0，方程组有唯一解

法一：求出A-1，x=A-1b

法二：利用初等行变换将[A| b]中A化为对角阵

计算矩阵行列式：直接调用Det[]函数计算

计算矩阵条件数：第一步：求出矩阵的逆矩阵

第二步：分别计算矩阵极其逆的无穷范数，一范数和二

范数

第三步：求解矩阵的条件数

Cond(A)∞=||A-1||∞*||A||∞

Cond(A)1=||A-1||1*||A||1

Cond(A)2=||A-1||2*||A||2=(λmax(A的转置*A)/

λmin(A的转置*A))1/2

计算机求解

第一步：消元

对ｋ＝１，２，３，……ｎ，进行：

步骤１：选主元（第ｋ列中第ｋ个至第ｎ个元素中绝对值较大者）

步骤２：将主元所在行与第ｋ行交换

步骤３：消元

第二步：回代求解

流程图如下：k=1,2,……,n

i=k+1,k+2,……,n

l

ik =a

ik

/a

kk

得到a

ik

j=k+1,k+2,……,n+1 a ij-a ik*a kj得到a ij

三、算例、应用实例

用列主元消去法解线性方程组Ax=b

⑴ 3.10x1+ 6.03x2+1.99x3=1

1.27x1+ 4.16x2-1.23x3= 1 ;

0.983x1-4.81x2+ 9.34x3=1

⑵ 3.00x1+ 6.03x2+ 1.99x3=1

1.27x1+ 4.16x2-1.23 x3 = 1.

0.990x1 -4.81x2+9.34 x3=1

分别输出A ,b ,detA,解向量x,⑴中A的条件数。分析比较⑴、⑵的计算结果。

输出结果为：

⑴

A =

3.100000000000000 6.030000000000000 1.990000000000000

1.270000000000000 4.160000000000000 -1.230000000000000

0.983000000000000 -4.810000000000000 9.340000000000000

b =

1

1

1

X=

-16.234490957625489

6.763269451277317

5.298697074088838

ans =

(矩阵的行列式的值)

2.9967

ans =

（A的条件数）

314.3810

⑵

A =

3.000000000000000 6.030000000000000 1.990000000000000

1.270000000000000 4.160000000000000 -1.230000000000000

0.990000000000000 -4.810000000000000 9.340000000000000

b =

1

1

1

X=

1.0e+002 *

1.195273381259593

-0.471426044312964

-0.368402561091259

ans =

-0.4070

虽然两题中A的数只差两个数，但是结果完全不同

五、参考文献

数值计算方法与算法（第二版）科学出版社

数值分析（第五版）清华大学出版社

六、附录

%高斯列主元消元法求解线性方程组Ax=b

%A为输入矩阵系数,b为方程组右端系数

%方程组的解保存在x变量中

format long;%设置为长格式显示,显示15位小数

A=[？？？]

det(A);

cond(A);

b=[？？？]'

[m,n]=size(A);

%先检查系数正确性

if m~=n

error('矩阵A的行数和列数必须相同');

return;

end

if m~=size(b)

error('b的大小必须和A的行数或A的列数相同');

return;

end

%再检查方程是否存在唯一解

if rank(A)~=rank([A,b])

error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');

return;

end

c=n+1;

A(:,c)=b; %(增广)

for k=1:n-1

[r,m]=max(abs(A(k:n,k))); %选主元

m=m+k-1; %修正操作行的值

if(A(m,k)~=0)

if(m~=k)

A([k m],:)=A([m k],:); %换行

end

A(k+1:n,k:c)=A(k+1:n,k:c)-(A(k+1:n,k)/A(k,k))*A(k, k:c); %消去end

end

x=zeros(length(b),1); %回代求解

x(n)=A(n,c)/A(n,n);

for k=n-1:-1:1

x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

disp('X=');

disp(x);

format short;%设置为默认格式显示,显示5位