高斯消元法.docx
《高斯消元法.docx》由会员分享，可在线阅读，更多相关《高斯消元法.docx（14页珍藏版）》请在冰豆网上搜索。
高斯消元法
求解线性方程组的直接解法
5.1 Gauss消去法
1三角方程组
先举一个简单的例子来说明消去法的基本思想.
例1.用消去法解方程组
解第一步.将方程
（1）乘上-2加到方程（3）上去,消去（3）中的未知数,得到
第二步.将方程
（2）加到方程（4）上去,消去方程（4）中的未知数,得到与原方程组等价的三角形方程组
显然,方程组（5）是容易求解的,解为
上述过程相当于
其中用表示矩阵的第行.
下面我们讨论求解一般线性方程组的高斯消去法.
一般地
当a11a22…ann≠0时,可解出
xn=bn/ann
fork=n-1:
1
xk=（b1-ak,k+1xk+1-…-aknxn）/akk
end
注:
可用同一组单元.并可解出一个未知数即代入其它方程消去该未知数
Gauss消元法的流程图为:
流程图中，分别为线性方程组的系数矩阵和常数向量；
k是循环次数。
②顺序消去法
一般地,k=1对n阶方程组消去第k个元（akk≠0）:
这里各行变换:
i行-k行×mik,其中mik=aik/akk,i=k+1,…,n而后,k=2对n-1阶方程组消第2个元…我们有如下顺序消元算法:
fork=1:
n-1
ifakk≠0
fori=k+1:
n
mik=aik/akki行=i行-k行×mik
end
elsestopend
end
（每行包括右端项!
）可细化,也可存储mik于aik得:
fork=1:
n-1
ifakk≠0
fori=k+1:
n
aik=aik/akk
forj=k+1:
n
aik=aik-aik×akj
end
bi=bi-aik×bk
end
elsestopend
end
顺序消元过程和回代过程连起来就可得精确解.顺序消元算法也可将系数矩阵和右端项分开:
fork=1:
n-1
ifakk≠0
fori=k+1:
n
aik=aik/akk
forj=k+1:
n
aik=aik-aik×akj
end
end
elsestopend
end
（注意mik在aik）
fork=1:
n-1
fori=k+1:
n
bi=bi-aik×bk
end
end
Gauss消去法运算量
消去第k个元素时，对矩阵作加法和乘法运算各（n-k）×（n-k）次,除法（n-k）次.对右端作加法和乘法运算各（n-k）次.分别共12+22+…+（n-1）2=n（n-1/2）（n-1）/3和1+2+…+（n-1）=n（n-1）/2次加法乘法,消元时还有1+2+…+（n-1）=n（n-1）/2次除法.另外回代过程中加法和乘法运算各n（n-1）/2次,除法n次.运算量主要是消元的贡献,加法和乘法运算各约n3/3.
定理1.设Ax=b,其中A
（1）如果则可通过Gauss消去法将Ax=b约化等价的三角形方程组,
且计算公式为:
（a）消元过程
设,对计算
（b）回代过程
（2）如果A为非奇异矩阵,则可通过Gauss消去法（及交换两行的初等变换）将方程组Ax=b约化为
.
行列式和逆矩阵
易知顺序消元过程中，行列式不变.因此det（A）=det（U）=u11u22…unn,这里U是顺序消元过程结束时的上三角矩阵.A的顺序主子式。
也是同样情况,即有det（Ak）=det（Uk）=u11u22…ukk.。
顺便指出,消第k个元素时，左上角的元素（称之为主元素）非零时才能进行.所有主元素非零等价于det（Ak）≠0,k=1,2,…,n-1.
求A的逆矩阵相当于解n个方程组,它们的系数矩阵是A,n个右端项是n个单位向量.用增广矩阵表示即（A∣I），其中I是单位矩阵.消元过程可对所有右端项同时进行,消元过程结束再分别回代.也可巧妙安排就地求逆.
③主元素法
如果akk=0顺序消元无法进行,如果akk很小，可以进行。
但将引出很大误差.
例2.用顺序消去法解方程组
解用精确运算:
解得（10000/9999,9998/9999）≈（1.00010001000100,0.99989998999900）
若在十进三位尾数舍入的浮点计算机系统中运算,第二行将是
（0-10000∣-10000）
从而得到解,与实解相去甚远.究其原因,小主元作分母扩大误差所致.如果两个方程（两行）交换次序再消元:
得到解,与实解很近.
下面引进的列主元素法在消元前，先选该列中绝对值最大的做主元（交换两行,每行包括右端项!
）:
fork=1:
n-1
找p:
p行k行
ik=p
ifakk≠0
fori=k+1:
n
mik=aik/akk
i行=i行-k行×mik
end
elsestopend
end
ik=p记录每次选主元（两行交换）的指标.跟前面一样,算法还可细化,把矩阵和右端分开算（这时就要用指标ik）。
列主元素法乘数mik绝对值不大于1,不会增加误差。
列主元素法用来求行列式时，要注意两行交换行列式变号。
还有一种全主元素法,它是在整个右下（n-k）×（n-k）矩阵找绝对值最大的做主元（交换行及列）.这对误差控制有利,但搜索太费时.通常列主元素法误差控制就已可以了.
以下是列主元素法的流程图:
流程图中，分别为线性方程组的系数矩阵和常熟向量；k是循环次数；d是主元素；
二实验部分
本章实验内容：
实验题目：
Gauss消元法,追赶法,范数。
实验内容：
①编制用Gauss消元法求解线性方程组Ax=f的程序。
②编制用追赶法求解线性方程组Ax=f的程序。
③编制向量和矩阵的范数程序。
实验目的：
①了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。
②掌握追赶法，能利用追赶法求解线性方程组。
③理解向量和矩阵范数定义,性质并掌握其计算方法.
编程要求：
利用Gauss消元法,追赶法解线性方程组。
分析误差。
计算算法:
①Gauss消元法:
1.消元过程
设,对计算
⒉回代过程
②追赶法：
1.分解Ax=f:
（）
2.解Lg=f,求g:
（）
3．解Ux=g,求x:
（）
③范数：
常用向量范数有：
（令x=（x1,x2,…,xn））
1-范数：
║x║1=│x1│+│x2│+…+│xn│
2-范数：
║x║2=（│x1│2+│x2│2+…+│xn│2）1/2
∞-范数：
║x║=max（│x1│,│x2│,…,│xn│）
常用的三种向量范数导出的矩阵范数是：
1-范数:
║A║1=max{║Ax║1/║x║1=1}=
2-范数:
║A║2=max{║Ax║2/║x║2=1}=,λ1是ATA的最大特征值.
∞-范数:
║A║=max{║Ax║/║x║=1}=
实验例题⑴：
用Gauss消元法解方程组
实验例题⑵:
用追赶法解三对角方程组Ax=b,其中
实验例题⑶：
设
，
计算A的行列范数.
程序①:
Gauss消元法
functionx=Gauss（A,b）
%A是线性方程组的系数矩阵,b为自由项.
n=length（A）
fork=1:
n-1
m（k+1:
n,k）=A（k+1:
n,k）/A（k,k）;
A（k+1:
n,k+1:
n）=A（k+1:
n,k+1:
n）-m（k+1:
n,k）*A（k,k+1:
n）;
b（k+1:
n）=b（k+1:
n）-m（k+1:
n,k）*b（k）;
end
x=zeros（n,1）;
x（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1
x（k）=（b（k）-A（k,k+1:
n）*x（k+1:
n））/A（k,k）;
end
x=x';
disp（sprintf（'kx（k）'））;
fori=0:
n
disp（sprintf（'%d%f',i,x（i+1）））;
end
数值结果：
x=Gauss（A,b）
n=3
kx（k）
01.111111
10.777778
22.555556
程序②:
追赶法
function[x,y,beta]=zhuiganfa（a,b,c,f）
%a,b,c是三对角阵的对角线上的元素,f是自由项.
n=length（b）;
beta
（1）=c
（1）/b
（1）;
fori=2:
n
beta（i）=c（i）/（b（i）-a（i）*beta（i-1））;
end
y
（1）=f
（1）/b
（1）;
fori=2:
n
y（i）=（f（i）-a（i）*y（i-1））/（b（i）-a（i）*beta（i-1））;
end
x（n）=y（n）;
fori=n-1:
-1:
1
x（i）=y（i）-beta（i）*x（i+1）;
end
disp（sprintf（'kx（k）y（k）beta（k）'））;
fori=0:
n
disp（sprintf（'%d%f',i,x（i+1）,y（i+1）,beta（i+1）））;
end
数值结果:
a=[0-1-1-1-1]';
b=[22222]';
c=[-1-1-1-10]';
f=[10000]';
[x,y,beta]=zhuiganfa（a,b,c,f）
kx（k）y（k）beta（k）
00.8333335.000000e-001-0.500000
10.6666673.333333e-001-0.666667
20.5000002.500000e-001-0.750000
30.3333332.000000e-001-0.800000
40.1666671.666667e-0010.000000
程序③：
1.列范数:
functionfan=lie（A）
%A为已知矩阵
n=length（A）
forj=1:
n
x（j）=0
fori=1:
n
x（j）=x（j）+abs（A（i,j））;
end
end
fan=max（x）
disp（sprintf（'nx（n）'））;
fori=0:
n
disp（sprintf（'%d%f',i,x（i+1）））;
end
数值结果:
fan=lie（A）
fan=0.8000
nx（n）
10.700000
20.800000
2.行范数:
functionfan=hang（A）
%A为已知矩阵
n=length（A）
fori=1:
n
x（i）=0
forj=1:
n
x（i）=x（i）+abs（A（i,j））;
end
end
fan=max（x）
disp（sprintf（'nx（n）'））;
fori=0:
n
disp（sprintf（'%d%f',i,x（i+1）））;
end
数值结果:
fan=hang（A）
fan=1.1000
nx（n）
11.100000
20.400000
总结:
从代数上看,直接分解法和Gauss消去法本质上一样,但如果我们采用”双精度累加”计算,那么直接三角分解法的精度要
升级会员 