解線性方程組的直接解法

第九講解線性方程組的直接解法內(nèi)容提要引言Gauss消元法列主元素消元法誤差分析MATLAB的線性方程組求解函數(shù)1小結(jié)1、引言例：小行星軌道問題： 天文學家要確定一小行星的軌道,在軌道平面建立以太陽為原點的直角坐標系.在坐標軸上取天文測量單位(一天文單位為地球到太陽的平均距離：9300萬哩),對小行星作5次觀察,測得軌道上5個點的坐標數(shù)據(jù)如下：

x5.76406.28606.75907.16807.4800

y0.64801.20201.82302.52603.3600

橢圓的一般方程:a1x2+a2xy+a3y2+a4x+a5y+1=0

將數(shù)據(jù)逐個代入,可得五個方程的方程組，求解該線性方程組即可得行星軌道方程。對一般線性方程組:Ax=b,其中當且僅當矩陣A行列式不為0時，即A非奇異時，方程組存在唯一解，可根據(jù)克萊姆法則求解，其算法設(shè)計如下：(1)

輸入系數(shù)矩陣A和右端向量b；(2)計算系數(shù)矩陣A的行列式值D，如果D=0，則輸出錯誤信息，結(jié)束，否則進行第(3)步；(3)

對k=1,2,···,n，用b替換A的第k列數(shù)據(jù)，并計算替換后矩陣的行列式值Dk；(4)

計算并輸出x1=D1/D，x2=D2/D，····,xn=Dn/D，結(jié)束。克萊姆法則只適用于低階方程組，高階方程組工作量太大，故一般用數(shù)值方法求解。數(shù)值方法分兩類：直接解法：在計算過程中，如果所有運算都是精確的，在理論上，經(jīng)過有限次運算就可以得到精確解，適用于變量較少的方程組。迭代解法：近似解法，運算次數(shù)因要求的計算精度而變化。迭代法也稱輾轉(zhuǎn)法，是一種不斷用變量的舊值遞推新值的過程。迭代算法是用計算機解決問題的一種基本方法。它利用計算機運算速度快、適合做重復(fù)性操作的特點，讓計算機對一組指令（或一定步驟）進行重復(fù)執(zhí)行，在每次執(zhí)行這組指令（或這些步驟）時，都從變量的原值推出它的一個新值。迭代方法的使用：確定迭代變量。建立迭代關(guān)系式。給定初值對迭代過程進行控制。對迭代結(jié)果判定2、Gauss消元法2.1基本思想：逐步消去未知元，將方程組化為與其等價的上三角方程組求解。2.2分兩步：第一步:消元過程，將方程組消元化為等價的上三角形方程組;第二步:回代過程，解上三角形方程組,得原方程組的解。2.3Gauss消元的目的a11x1+a12x2+····+a1nxn=b1a21x1+a22x2+····+a2nxn=b2·········································an1x1+an2x2+····+annxn=bn原始方程組約化方程組2.4.1消元過程(化一般方程組為上三角方程組)以四階為例：其系數(shù)增廣矩陣為：第一輪消元：計算3個數(shù):[m21

m31

m41]T=[a21

a31

a41]T/a11

用-m21乘矩陣第一行后加到矩陣第二行;

用-m31乘矩陣第一行后加到矩陣第三行;

用-m41乘矩陣第一行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋旱诙喯河嬎?個數(shù)：[m32

m42]T=[a32(1)

a42(1)]T

/a22(1)

用-m32乘矩陣第二行后加到矩陣第三行;

用-m42乘矩陣第二行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋旱谌喯河嬎?m43=a43(2)/a33(2)用-m43乘矩陣第三行后加到矩陣第四行;其系數(shù)增廣矩陣變?yōu)椋浩鋵?yīng)的上三角方程組為

若對于一般的線性方程組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消去法求解線性方程組%輸入?yún)?shù)：%---A：線性方程組的系數(shù)矩陣%---b：線性方程組的右端項%輸出參數(shù)：%---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的系數(shù)矩陣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的三角分解算法。重復(fù)該過程，最后得記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如果已經(jīng)有A=LU

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

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

U

是上三角矩陣,用回代算法記UX=Y

，LY=b

，則求解方程組分兩步進行：3、Gauss列主元素消元法例：求解方程組（用四位浮點數(shù)計算，精確解舍入到4位有效數(shù)字為：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)顯然，此解并不準確。方法二：交換行，避免絕對值小的主元作除數(shù)。（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基本思想：設(shè)Ax＝b的增廣矩陣為在A的第一列中選絕對值最大的元素作主元，設(shè)該元素所在行為第i1行，交換第一行與第i1行，進行第一次消元；再在第2－n行的第二列中選絕對值最大的元素作主元，設(shè)該元素所在行為第i2行，交換第二行與第i2行，進行第二次消元，……直到消元過程完成為止。例：用列主元法解解：第一列中絕對值最大是10，取10為主元第二列的后兩個數(shù)中選出主元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消去法求解線性方程組%輸入?yún)?shù)：%---A：線性方程組的系數(shù)矩陣%---b：線性方程組的右端項%輸出參數(shù)：%---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為上三角陣。矩陣分解關(guān)系為3.3全主元素消元法定義：則稱為全主元素。經(jīng)過行列互換，使得位于經(jīng)交換行和列后的等價方程組中的位置，然后再實施消元。

全主元素消元法的基本思想：若注：全主元素消元法有可能改變未知數(shù)的順序。4、誤差分析例記方程組（1）為Ax＝b，其精確解為：x1*=2，x2*=0現(xiàn)考察方程組（2）可將其表示為：A(x+x)=b+b，其中b=(0,0.0001)T

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

結(jié)論：（1）的常數(shù)項b的第二個分量只有1/1000的微小變化，