用高斯消元法的MATALAB程序设计.docx
《用高斯消元法的MATALAB程序设计.docx》由会员分享，可在线阅读，更多相关《用高斯消元法的MATALAB程序设计.docx（17页珍藏版）》请在冰豆网上搜索。
用高斯消元法的MATALAB程序设计
毕业论文
序言
高斯消元法，又称高斯消去法，实际上就是我们俗称的加减消元法。
数学上，高斯消去法或称高斯－约当消去法，由高斯和约当得名（很多人将高斯消去作为完整的高斯－约当消去的前半部分），它是线性代数中的一个算法，用于决定线性方程组的解，决定矩阵的秩，以及决定可逆方矩阵的逆。
当用于一个矩阵时，高斯消去产生“行消去梯形形式”。
该方法以数学家高斯命名，但最早出现于中国古籍《九章算术》，成书于约公元前150年。
一．高斯消元法
例如：
一个二元一次方程组，设法对每个等式进行变形，使两个等式中的同一个未知数的系数相等，这两个等式相减，得到一个新的等式，在这个新的等式中，细数相等的未知数就被除去了（系数为0。
同样的也适合多元多次方程组。
什么是线性方程组？
含m个方程和n个未知量的方程组定义为
a（11）x
（1）+a（12）x
（2）+...+a（1n）x（n）=b
（1）
a（21）x
（1）+a（22）x
（2）+...+a（2n）x（n）=b
（2）
...
a（m1）x
（1）+a（m2）x
（2）+...+a（mn）x（n）=b（m）
这个方程组称为m*n线性方程组，其中a（ij）和b（i）为实数，括号中为下标。
这个方程组有多种表示方法。
例如，我们知道m*n矩阵（用大写字母表示）是一个m行n列的数阵，n维向量（用加粗的小写字母表示）是n个数的数组，也就是一个n*1矩阵（列向量。
我们不考虑行向量）。
另外，大家也都知道矩阵乘法。
因此一个m*n线性方程组可以表示为Ax=b，其中A是由系数aij组成的m*n矩阵即系数矩阵，x是n维的未知数向量，b是m维的结果向量。
如果把向量b写到A的右边得到m*（n+1）的矩阵，得到的新矩阵称为这个方程组的增广矩阵。
每一个方程组均对应于一个增广矩阵。
下面介绍一下矩阵的初等行变换:
1交换两行
2用非零实数乘以任一行
3把某一行的倍数加到另一行上
同理可以定义初等列变换。
初等行变换和初等列变换统称初等变换。
定理：
对于一个方程组对应的增广矩阵进行有限次初等行变换，所得矩阵对应的方程组与原方程组等价。
即它们是等价方程组。
高斯消元的过程，就是利用初等行变换将原来不容易求解的方程组转化为容易求解的方程组。
下面我们以求解n*n线性方程组为例。
因为如果m<>n,这个m*n方程组一般不会有唯一的解。
定义：
上三角矩阵是一个n*n的矩阵，其中对于任意i>j的项a（ij）=0。
若i=j时a（ij）<>0则称为严格上三角矩阵。
下面就是一个严格上三角矩阵。
121
011
001
我们发现，如果将系数矩阵化为上三角矩阵，就可以轻而易举地求解。
可以证明，n*n方程组有唯一解等价于它的增广矩阵可以化为严格上三角矩阵。
消元过程如下：
对于增广矩阵A
1fori:
=1ton-1do
2选择第i至第n行中第i个元素绝对值最大的行，与第i行交换//选择主元，由于涉及实数除法运算，选取大绝对值可以减小误差
3forj:
=i+1tondo
4使用第3种行变换使a（j,i）=0
根据得到的矩阵，可以回代求出每个未知数x（i）
这只是对于有唯一解的情况。
那么，其他情况呢？
我们定义在消元过程中保留的那行称为主行，主行第一个非零元素称为主元。
消元的任务就是使主元非0，使主元下方元素变为0。
当消元过程中无法选择非0主元，则此方程组无解或无穷解。
此时，我们选择右边一列操作。
即，若第i行为主行，在第j列无法选择主元，就对第i+1列操作。
这样，最后的结果不再是上三角形，而是行阶梯形：
每一次在竖直方向下降1，在水平方向扩展可能大于1。
在增广矩阵中能够选择主元的列对应的变量称为首变量，跳过的列对应的变量称为自由变量。
定义：
行阶梯形矩阵是一个满足
（1）每一个非零行的第一个非零元为1;
（2）若第i行为非零行，则第i+1行行首的0多于第i行;（3）非零行在全零行之前
如这就是一个行阶梯形矩阵：
1234
0013
0001
当然，对于竞赛解方程没必要把首元素化为1。
定理：
m*n线性方程组有解当且仅当其行阶梯形矩阵不含这样一行：
[00...01]
亦即，对于n*n线性方程组，若没有唯一解，我们就将其化为行阶梯形，若最后某行形如[00...0a]（a<>0）则无解，否则有无穷多的解。
上面我们主要讨论了n*n方程组。
下面讨论m<>n的方程组：
1若m>n则这个方程组称为超定方程组。
超定方程组对应的矩阵行数增加，因此通常无解（但若其行阶梯形矩阵上方为严格上三角，下方全为[00...0]，则有唯一解；若下方全为[00...0]而上方是阶梯形，则有无穷解）
2若m 由此可见，对于任意m*n线性方程组，求解均将其转化为行阶梯形矩阵。
这就是我们所讨论的：
定义利用矩阵行初等变换，将线性方程组的增广矩阵化为行阶梯形的过程称为高斯消元法（Gaussianelimination）。
而求解线性方程组的步骤：
1将其增广矩阵化为行阶梯形
2若最后有形如[00...0a]（a<>0）的行则无解否则
3若含有自由变量则有无穷组解否则
4原方程有唯一解。
采用回代求解。
二．关于矩阵可对角化的问题
有限维线性空间上的线性变换可对角化的问题与矩阵可对角化的问题有密切的联系。
设
是数域P上的n维线性空间V的一个线性变换，若存在V的一个基，似的
关于这个基的矩阵是对角形矩阵，则呈线性变换
可对角化。
类似的，设A是数域P上的一个n阶矩阵，若A能与一个对角型矩阵相似，则称矩阵A可对角化。
由于线性变换
关于线性空间V的不同基的矩阵是相似矩阵，从而有：
若A是线性变换
关于线性变换V的某个基的矩阵，那么
可对角化当且仅当A可对角化。
为简便起见,我们只讨论矩阵可对角化的问题。
设A
A可对角化的充分条件是A在数域P上有n个不同的特征值，充要条件是A在数域P上有n个线性无关的特征向量。
因此，通常判定矩阵A可对角化的方法是：
首先求矩阵A的特征值，若特征多项式
在数域P上有n个不同的特征根，则可判定A可对角化；若特征多项式
在数域P上有重根，则需确定A的属于每个特征值的特征向量，若A有n个线性无关的特征向量,则A可对角化,否则A不可对角化.
定义：
设A
若存在多项式
使得f（A）=0,则称f（x）为矩阵A的化零多项式，其中首项系数为1的最低次化零多项式成为A的最小多项式。
由hamlton-caylay定理可知,矩阵A的特征多项式是A的化零多项式,因此,A的最小多项式是A的特征多项式的一个因式。
性质最小多项式的性质有
1.矩阵A的最小多项式是惟一的
2.设f（x）为矩阵A的最小多项式，则g（x）为A的化零多项式的充要条件是f（x）|g（x）;
3.相似矩阵有相同的最小多项式
求矩阵最小多项式的方法：
1.利用矩阵的特征多项式，因为矩阵的最小多项式是他的特征多项式的因式。
2.利用矩阵的线性相关性质，线性空间
中向量组I，A,
，….
的向量数为
+1,必线性相关，因而，存在一个最小正整数k，使得
，（
3.利用矩阵的特征矩阵的标准型，通常，人们利用矩阵的特征矩阵在初等变换下的标准型来求矩阵的若当标准型，我们也可以利用矩阵的特征矩阵的标准型来求他的最小多项式，矩阵A的特征矩阵
经初等变换化为标准型后，其主对角线上非零元素称为
的不变因子，简称A的不变因子，则最后一个不变因子就是A的最小多项式。
4.利用向量的线性相关性质，将有限维线性空间使我线性变换的最小多项式，我们将它转换为讨论矩阵的最小多项式并加以证明。
设A
，v
在线性空间
上向量组v,Av,
必线性相关，因而存在一个最小整数
使得
令
则
是使得
的最低首次一多项式。
定理一：
设A
且
线性无关。
令
的最低首一多项式，则矩阵A的最小多项式f（x）是
的最小公倍式，即f（x）=
.
证：
设g（x）=
.
则
又
线性无关，
即g（x）是A的化零多项式，由最小多项式的性质，我们有
反之，对于j=1,2，…,n,设
，则有
，但
是使得
的最低首一多项式，从而
.因此
的最小公倍式也整除f（x）,即g（x）|f（x）,又f（x）,g（x）都是首一多项式，故
。
注：
在求矩阵的最小多项式时，向量
通常就取n维单位向量
.
定理2数域P上N阶矩阵A可对角化的充要条件为A的最小多项式是P上互质的一次因式的乘积。
三．高斯消元法的流程头
四．高斯消元法的程序设计
functionx=DelGauss（a,b）
%Gauss消去法
[n,m]=size（a）;
nb=length（b）;
det=1;%存储行列式值
x=zeros（n,1）;
fork=1:
n-1
fori=k+1:
n
ifa（k,k）==0
return
end
m=a（i,k）/a（k,k）;
forj=k+1:
n
a（i,j）=a（i,j）-m*a（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*a（k,k）;%计算行列式
end
det=det*a（n,n）;
fork=n:
-1:
1%回代求解
forj=k+1:
n
b（k）=b（k）-a（k,j）*x（j）;
end
x（k）=b（k）/a（k,k）;
end
例1：
编写高斯消元法MATLAB文件如下：
clear;
A=[521;28-3;1-3-6];
b=[8;21;1;];
[RA,RB,n,X]=gaus（A,b）
运行结果为：
请注意：
因为RA=RB=n，所以此方程组有唯一解.
RA=
3
RB=
3
n=
3
X=
1
2
-1
例2：
>>A=[1.0170-0.00920.0095;-0.00920.99030.0136;0.00950.01360.9898];
>>b=[101]';
>>x=DelGauss（A,b）
x=
0.9739
-0.0047
1.0010
五．由高斯总结的LU矩阵的三角形分解的程序设计及其它分解的程序设计
（一）LU矩阵的三角分解算法
M文件
functiony=LU（A,B）;
n=length（A）;
A=[AB];
fork=1:
n-1;
fori=k:
n;
if（abs（A（i,k））==max（abs（A（k:
n,k））））
P（k）=i;
temp=A（k,:
）;
A（k,:
）=A（i,:
）;
A（i,:
）=temp;
end
end
forj=k+1:
n;
A（j,k）=A（j,k）/A（k,k）;
A（j,k+1:
n+1）=A（j,k+1:
n+1）-A（j,k）*A（k,k+1:
n+1）;
end
end
P（n）=n;
L（1,1）=1;
L（2:
n,1）=A（2:
n,1）;
L（1,2:
n）=0;
U（1,1）=A（1,1）;
U（2:
n,1）=0;
U（1,2:
n）=A（1,2:
n）;
fori=2:
n;
L（i,1:
i-1）=A（i,1:
i-1）;
L（i,i）=1;
L（i,i+1:
n）=0;
U（i,1:
i-1）=0;
U（i,i:
n）=A（i,i:
n）;
end
x（n）=A（n,n+1）/U（n,n）;
fork=n-1:
-1:
1
x（k）=A（k,n+1）;
forp=n:
-1:
k+1;
x（k）=x（k）-U（k,p）*x（p）;
end
x（k）=x（k）/U（k,k）;
end
x
L
U
P
End
Matlab运行结果
（二）龙贝格（Romberg）算法
M文件
function[t]=romberg（f,a,b,e）
t=zeros（15,4）;
t（1,1）=（b-a）/2*（f（a）+f（b））;
fork=2:
4
sum=0;
fori=1:
2^（k-2）
sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;
end
t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;
fori=2:
k
t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;
end
end
fork=5:
15
sum=0;
fori=1:
2^（k-2）
sum=sum+f（a+（2*i-1）*（b-a）/2^（k-1））;
end
t（k,1）=0.5*t（k-1,1）+（b-a）/2^（k-1）*sum;
fori=2:
4
t（k,i）=（4^（i-1）*t（k,i-1）-t（k-1,i-1））/（4^（i-1）-1）;
end
ifk>6
ifabs（t（k,4）-t（k-1,4））disp（['答案',num2str（t（k,4））]）;
break;
end
end
end
ifk>=15
disp（['溢出']）;
end
Matlab运行结果
（三）列主元Gauss消去法：
functionx=detGauss（a,b）
%Gauss列主元消去法
[n,m]=size（a）;
nb=length（b）;
det=1;%存储行列式值
x=zeros（n,1）;
fork=1:
n-1
amax=0;%选主元
fori=k:
n
ifabs（a（i,k））>amax
amax=abs（a（i,k））;r=i;
end
end
ifamax<1e-10
return;
end
ifr>k%交换两行
forj=k:
n
z=a（k,j）;a（k,j）=a（r,j）;a（r,j）=z;
end
z=b（k）;b（k）=b（r）;b（r）=z;det=-det;
end
fori=k+1:
n%进行消元
m=a（i,k）/a（k,k）;
forj=k+1:
n
a（i,j）=a（i,j）-m*a（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*a（k,k）;
end
det=det*a（n,n）;
fork=n:
-1:
1%回代求解
forj=k+1:
n
b（k）=b（k）-a（k,j）*x（j）;
end
x（k）=b（k）/a（k,k）;
end
Example:
>>x=detGauss（A,b）
x=
0.9739
-0.0047
1.0010
升级会员 