线性方程组求解直接法

解线性方程组的直接法引言消元法：Gauss消元法消元法：列主元素消元法LU分解：求解三对角方程组的追赶法误差分析2.1 引言小行星轨道问题：

天文学家要确定一小行星的轨道,在轨道平面建立以太阳为原点的直角坐标系.在坐标轴上取天文测量单位(一天文单位为地球到太阳的平均距离：9300万哩),对小行星作5次观察,测得轨道上5个点的坐标数据如下：

x

5.76406.28606.75907.16807.4800

y0.64801.20201.82302.52603.3600椭圆的一般方程:a1x2+a2xy+a3y2+a4x+a5y+1=0将数据逐个代入,可得五个方程的方程组，求解该线性方程组即可得行星轨道方程。对一般线性方程组:Ax=b,其中2.1 引言由以前所学内容知，当且仅当矩阵A行列式不为0时，即A非奇异时，方程组存在唯一解，可根据克莱姆法则求解，其算法设计如下：(1)

输入系数矩阵A和右端向量b；(2)计算系数矩阵A的行列式值D，如果D=0，则输出错误信息，结束，否则进行第(3)步；(3)

对k=1,2,···,n，用b替换A的第k列数据，并计算替换后矩阵的行列式值Dk；(4)

计算并输出x1=D1/D，x2=D2/D，····,xn=Dn/D，结束。但克莱姆法则只适用于低阶方程组，高阶方程组工作量太大，故一般用数值方法求解。数值方法分两类：1.直接法2.迭代法2.2Gauss消元法基本思想：逐步消去未知元，将方程组化为与其等价的上三角方程组求解。分两步：第一步:消元过程，将方程组消元化为等价的上三角形方程组;第二步:回代过程，解上三角形方程组,得原方程组的解。Gauss消元的目的a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bn原始方程组约化方程组注意：引入符号的上标含义消元过程(化一般方程组为上三角方程组)以四阶为例：其系数增广矩阵为：第一轮消元：计算3个数:[m21

m31

m41]T=[a21

a31

a41]T/a11用-m21乘矩阵第一行后加到矩阵第二行;用-m31乘矩阵第一行后加到矩阵第三行;用-m41乘矩阵第一行后加到矩阵第四行;其系数增广矩阵变为：消元过程(化一般方程组为上三角方程组)第二轮消元：计算2个数：[m32

m42]T=[a32(1)

a42(1)]T

/a22(1)

用-m32乘矩阵第二行后加到矩阵第三行;用-m42乘矩阵第二行后加到矩阵第四行;其系数增广矩阵变为：消元过程(化一般方程组为上三角方程组)第三轮消元：计算:m43=a43(2)/a33(2)用-m43乘矩阵第三行后加到矩阵第四行;其系数增广矩阵变为：消元过程(化一般方程组为上三角方程组)其对应的上三角方程组为消元过程(化一般方程组为上三角方程组)若对于一般的线性方程组Ax=b，其消元过程的计算公式为：（k=1,2,…,n-1）消元过程(化一般方程组为上三角方程组)回代过程(解上三角方程组)上三角方程组的一般形式为：其中a11…ann≠0回代过程的计算公式：回代过程(解上三角方程组)工作量计算：消去过程：“÷”：第k步，n-k次，共

(n-1)+(n-2)+……+1=n(n-1)/2“×”：第k步，(n-k)(n-k+1)次，共

(n-1)n+(n-2)(n-1)+……+1×2=(n3-n)/3总工作量：

s1=n(n-1)/2+(n3-n)/3回代过程：“÷”：n“×”：1+2+……+(n-1)=n(n-1)/2总工作量：s2=n+n(n-1)/2=n(n+1)/2故Gauss消元法的总工作量为： s=s1+s2=n2+(n3-n)/3克莱姆法则求解的工作量为：“×”：（n+1个n阶行列式的值）(n+1)(n-1)n!“÷”：n故总工作量为：[(n+1)(n-1)]n!+n当n=6时,Gauss消元法工作量为106;而克莱姆法则求解工作量为25206。定理：

约化的主元素ak+1,k+1(k)≠0(k=0,1,···,n-1)的充分必要条件是矩阵A的各阶顺序主子式不为零。即注：对角线上的元素ak+1,k+1(k)在Gauss消去法中作用突出，称约化的主元素。推论:

如果A的顺序主子式Dk

≠0(k=1,···,n-1)，则Gauss消元法中的约化主元可以表示为例

用高斯消元法求解方程组x1=-13,x2=8,

x3=2m21=3/2m31=4/2m32=-3/0.5矩阵的三角分解：对线性方程组Ax=b的系数矩阵A施行初等行变换相当于用初等矩阵左乘A，故第一次消元后方程组化为A(1)x=b(1)，即L1A=A(1)x，L1b=b(1)，其中同理LkA(k-1)=A(k)

Lkb(k-1)=b(k)其中将A

分解为单位下三角矩阵L和上三角矩阵U的乘积的算法称为矩阵A的三角分解算法。重复该过程，最后得记U=A(n-1)，则其中定理：设A为n阶矩阵，若A的顺序主子式Di≠0（i=1,2,…n-1）,则A可分解为一个单位下三角矩阵L和一个上三角矩阵U的乘积，且这种分解是唯一的。由Gauss消元过程可推得

L＝

U即为Gauss消元后所得的上三角方程的系数矩阵。例

对矩阵A＝作LU分解。解由Gauss消去法可得， m21=0，m31=2，m32=－1。 故

A＝＝LU如果已经有A=LU

则AX=b=>LUX=b，（1）求解方程组LY=b得向量Y的值；

L是下三角矩阵,用顺代算法（2）求解方程组UX=Y得向量X的值。 U

是上三角矩阵,用回代算法记UX=Y

，LY=b

，则求解方程组分两步进行：基本思想：Gauss消元法中，若主元akk(k)

太小会使误差增大，故应避免采用绝对值小的元素作主元。最好每一步选取系数矩阵中（或消元后的低阶矩阵中）绝对值最大的元素作主元，以具较好的数值稳定性。2.3Gauss列主元素消元法例：求解方程组（用四位浮点数计算，精确解舍入到4位有效数字为：x1*=-0.4904,x2*=-0.05104,x3*=0.3675）解：《方法一》Gauss消去法（A/b）=其中，

m21=-1.000/0.001=-1000

m31=-2.000/0.001=-2000

m32=4001/2004=1.997 解为x1=-0.4000,x2=-0.09980,x3=0.4000(x1*=-0.4904,x2*=-0.05104,x3*=0.3675)显然，此解并不准确。《方法二》交换行，避免绝对值小的主元作除数。（A/b）=其中，m21=0.5000

m31=-0.0005

m32=0.6300 解为x1=-0.4900,x2=-0.05113,x3=0.3678(x1*=-0.4904,x2*=-0.05104,x3*=0.3675)与方法一相比，此解显然要精确得多。基本思想：设Ax＝b的增广矩阵为2.3Gauss列主元素消元法在A的第一列中选绝对值最大的元素作主元，设该元素所在行为第i1行，交换第一行与第i1行，进行第一次消元；再在第2－n行的第二列中选绝对值最大的元素作主元，设该元素所在行为第i2行，交换第二行与第i2行，进行第二次消元，……直到消元过程完成为止。例：用列主元法解解：第一列中绝对值最大是10，取10为主元第二列的后两个数中选出主元2.5x3=6.2/6.2=1x2=(2.5-5x3)/2.5=-1x1=(7+7x2-0x3)/10=0x1=0x2=-1x3=1列主元矩阵的三角分解：解：交换行变换例：对矩阵A做列主元三角分解：则列主元的Gauss变换可记为A(2)=F2·P2·F1·P1·A记U=A(2)=F2·(P2F1P2)·P2P1·A（因P2·P2=I）P=P2·P1则有对于一般的n阶矩阵的列主元三角分解，通常令定理：（列主元素的三角分解定理）若A非奇异，则存在排列阵P使PA＝LU，其中L为下三角阵，U为上三角阵。矩阵分解关系为经过行列互换，使得位于经交换行和列后的等价方程组中的位置，然后再实施消元。全主元素消元法定义：则称为全主元素。全主元素消元法的基本思想：若注：全主元素消元法有可能改变未知数的顺序。2.4求解三对角方程组的追赶法

且aici≠0，(i=2，3，…，n-1)

其中三对角方程组形式如下：2.4求解三对角方程组的追赶法基本思想：每一轮消元只对增广矩阵中某两行进行——将主对角元化为1，主对角元下方元素化为零。 即将其系数增广矩阵化为如下形式：解三对角方程组Ax＝b的追赶法分“追”和“赶”两个环节：追的过程：（消元过程）按顺序计算系数β1->β2->……->βn-1和y1->y2->……->yn;赶的过程：（回代过程）按逆序求出xn->xn-1->……->x2->x1。2.4求解三对角方程组的追赶法其系数增广矩阵为：

β1=c1/b1

y1=f1/b1

βi=ci/(bi–aiβi–1),(i=2,3,…,n-1)

yi=(fi–aiyi–1)/(bi–aiβi–1),(i=2,3,…,n)总工作量为：5n-42.5误差分析例

记方程组（1）为Ax＝b，其精确解为：x1*=2，x2*=0现考察方程组（2）可将其表示为：A(x+x)=b+b，其中

b=(0,0.0001)T

，x为（1）的解。显然（2）的解为：x+x=(1,1)T

结论：（1）的常数项b的第二个分量只有1/1000的微小变化，方程组的解变化却很大。定义

若矩阵A或常数项b的微小变化引起方程组Ax=b的解的巨大变化，则称此方程组为病态方程组，A为病态矩