解线性方程组的直接解法

第九讲解线性方程组的直接解法内容提要引言Gauss消元法列主元素消元法误差分析MATLAB的线性方程组求解函数1小结1、引言例：小行星轨道问题： 天文学家要确定一小行星的轨道,在轨道平面建立以太阳为原点的直角坐标系.在坐标轴上取天文测量单位(一天文单位为地球到太阳的平均距离：9300万哩),对小行星作5次观察,测得轨道上5个点的坐标数据如下：

x5.76406.28606.75907.16807.4800

y0.64801.20201.82302.52603.3600

椭圆的一般方程:a1x2+a2xy+a3y2+a4x+a5y+1=0

将数据逐个代入,可得五个方程的方程组，求解该线性方程组即可得行星轨道方程。对一般线性方程组:Ax=b,其中当且仅当矩阵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，结束。克莱姆法则只适用于低阶方程组，高阶方程组工作量太大，故一般用数值方法求解。数值方法分两类：直接解法：在计算过程中，如果所有运算都是精确的，在理论上，经过有限次运算就可以得到精确解，适用于变量较少的方程组。迭代解法：近似解法，运算次数因要求的计算精度而变化。迭代法也称辗转法，是一种不断用变量的旧值递推新值的过程。迭代算法是用计算机解决问题的一种基本方法。它利用计算机运算速度快、适合做重复性操作的特点，让计算机对一组指令（或一定步骤）进行重复执行，在每次执行这组指令（或这些步骤）时，都从变量的原值推出它的一个新值。迭代方法的使用：确定迭代变量。建立迭代关系式。给定初值对迭代过程进行控制。对迭代结果判定2、Gauss消元法2.1基本思想：逐步消去未知元，将方程组化为与其等价的上三角方程组求解。2.2分两步：第一步:消元过程，将方程组消元化为等价的上三角形方程组;第二步:回代过程，解上三角形方程组,得原方程组的解。2.3Gauss消元的目的a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bn原始方程组约化方程组2.4.1消元过程(化一般方程组为上三角方程组)以四阶为例：其系数增广矩阵为：第一轮消元：计算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）2.4.2回代过程(解上三角方程组)上三角方程组的一般形式为：其中a11…ann≠0回代过程的计算公式：2.5工作量计算：消去过程：“÷”：第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。function[U,x]=Gauss(A,b)%顺序Gauss消去法求解线性方程组%输入参数：%---A：线性方程组的系数矩阵%---b：线性方程组的右端项%输出参数：%---U：消元后的上三角方程组的增广矩阵%---x：线性方程组的解n=length(b);%消元过程fork=1:n-1m=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);endU=[A,b];%回代过程x=zeros(n,1);x(n)=b(n)/A(n,n);%求x_nfork=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);%求x_k，k=n-1,n-2,…,1end定理：约化的主元素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)，则其中>>clearall;A=[234;352;4330];b=[6532]';[L,U]=lu(A);x=U\(L\b)x=-13.00008.00002.0000MATLAB处理：例对矩阵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

，则求解方程组分两步进行：3、Gauss列主元素消元法例：求解方程组（用四位浮点数计算，精确解舍入到4位有效数字为：x1*=-0.4904,x2*=-0.05104,x3*=0.3675）解：方法一Gauss消去法（A/b）=其中，

m21=-1.000/0.001=-1000 m31=-2.000/0.001=-2000m32=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.0005m32=0.6300 解为x1=-0.4900,x2=-0.05113,x3=0.3678(x1*=-0.4904,x2*=-0.05104,x3*=0.3675)与方法一相比，此解显然要精确得多。3.1基本思想：设Ax＝b的增广矩阵为在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=1function[U,x]=Gauss_column(A,b)%列主元Gauss消去法求解线性方程组%输入参数：%---A：线性方程组的系数矩阵%---b：线性方程组的右端项%输出参数：%---U：消元后的上三角方程组的增广矩阵%---x：线性方程组的解n=length(b);fork=1:(n-1)%选主元

[Y,I]=max(abs(A(k:n,k)));I=I+k-1;ifI>kt=A(k,:);A(k,:)=A(I,:);A(I,:)=t;t=b(k);b(k)=b(I);b(I)=t;end%消元

m=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);endU=[A,b];%回代x=zeros(n,1);x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);end3.2列主元矩阵的三角分解解：交换行变换例：对矩阵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为上三角阵。矩阵分解关系为3.3全主元素消元法定义：则称为全主元素。经过行列互换，使得位于经交换行和列后的等价方程组中的位置，然后再实施消元。

全主元素消元法的基本思想：若注：全主元素消元法有可能改变未知数的顺序。4、误差分析例记方程组（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的微小变化，