【正文】
方程組化為 2 2222 2( 2 ) m a x 2 .ir inr a a r????找 , 使 , 對 調(diào) 行22 222222220 ( 3 , 4 , , ) :2 3 , 4 , 2 , 3 , 1iiiij j ija a i naiaaa a a i n j na???? ? ???????? ? ? ? ? ????? 消 元 : 用 把 消 為 第 行 第 行 , 則( , ; , )1 1 1 1 2 2 1 3 3 1 , 112 2 2 2 3 3 2 2 , 1 3 3 3 , 133 3 nnnn n nn n nna x a x a x a x aa x a x a x aa x a x aa???? ? ? ? ?? ? ? ?? ? ? LLLLLLLLLLLLLLLLLL,13 n n n n nx a x a????????? ? ? ??L(上三角方程組 ) ( ) (n1) 原方程組化為 以上為消元過程。 0] %先定義增廣矩陣 x = [0,0,0]‘。 a %第一次選主元(第一行與第二行交換) a(2,:) = a(2,:) a(1,:)*a(2,1)/a(1,1)。 a(2,:)=tempo。 x(1) = (a(1,4) a(1,2:3)*x(2:3))/a(1,1)。 例 2:用 Gauss- Jordan消去法求解上例中的矩陣 的逆矩陣。a tempo = a(2,:)。 a(i,:)=a(i,:) a(i,1)*a(1,:)。 a(2,:)=tempo。 if i~=3, a(i,:)=a(i,:)a(i,3)*a(3,:)。 1 , nxxL四 、 Gauss消元法的應(yīng)用 求行列式: det(A) 1 1 1 1 1 1111. . . . . .. . . . . . . . . . . . ( 1 ) . . . . . . ( 1 ) . . ....nnmmnnn n n n na a b bbba a b? ? ? ? ?1( ) . . . ( )A E E A ??? 求逆矩陣: inv( A) 或 A^(1) rref(A,eye(size(A))) 將 (A,E)化為行最簡形 , 其實就是 GaussJordon消去法 ( 3) 求解方程組 Ax=b ?求逆矩陣的思想 : x=inv(A)*b或 x=A^(1)*b ?左除法 ( 原理上是運用高斯消元法求解 ,但 Matlab在實際執(zhí)行過程中是通過 LU分解法進行的 ) : x=A\b ?符號矩陣法 ( 此法最接近精確值 , 但運算速度慢 )x=sym(A)\sym(b) ?化為行最簡形: C=[A,b] rref(C) 例 1:在 MATLAB上,用 Gauss列主元消去法求解方程組: 123 4 4 2 3 6 6 2 1 4 4 8 0xxx? ??? ? ? ???? ? ? ??? ??? ? ? ?? ? ? ?????? ? ? ???A= [ 。1。 下三角 , 若 L為單位下三角陣(對角元全為 1), U為上三角陣,則稱 A=LU 為 Doolittle分解 。 實際上,將選主元 Gauss消去法里的行交換同樣作用于單位矩陣,所得矩陣即為 P。6 15 18 40]%先錄入系數(shù)矩陣 b=[9 23 22 47]’ [L,U]=lu(A) %對 A矩陣做 LU分解 y=L\b %求解方程組 Ly=b x=U\y %求解方程組 Ux=y得到方程組的最終解 故方程組的最終解為 : x =(,)’ 或 A=[2 4 2 6。 A=[3, 1, 4, 1。 b=[12, 6, 4, 0]39。R=X 若 X為非對稱正定,則輸出一個出錯信息 ? [R, p]=chol(X):這個命令格式將不輸出出錯信息。 【 例 】 用 Cholesky分解求解線性方程組。1, 6, 1, 3]。解法如下: 由第一個方程 1 1 1 2 1b x c x d?? 1 1 211d c xxb???令 39。1d c xxb???將其代入第二個方程 ,得: 2 1 2 2 2 3 2a x b x c x d? ? ?22 1 2 3122211ad d c xbxabcb???????21,arb??39。3d c xxb???再令 3 3 2 ,d d r d????3 3 2b b r c? ???以此類推, …… 令 1,iiarb ???1i i iiid c xxb?? ????1i i ib b r c ?? ??1 ,i i id d r d ??????( 2, 3, , 1 )in??( 1 , 2, 3, , 1 )in??將 1nx ?代入最后一個方程 1n n n n na x b x d? ??令 nnndxb????1 ,n n nd d r d ?????1,nnarb ???以上過程稱為 “追”; 1n n nb b r c ?? ???總結(jié)“追”的算法: Step1:(初始化變量) 1,iiarb ???1i i ib b r c ?? ??1 ,i i id d r d ??????( 2 , 3 , , )in?39。 b(i)=b(i)r*c(i1)。 迭代法適用于解大型稀疏方程組 (萬階以上的方程組 ,系數(shù)矩陣中零元素占很大比例 ,而非零元按某種模式分布 ) 背景 : 電路分析、邊值問題的數(shù)值解和數(shù)學(xué)物理方程 問題 : (1)如何構(gòu)造迭代格式? (2)迭代格式是否收斂? (3)收斂速度如何? (4)如何進行誤差估計? 167。 s=[1,3,7]。 (i,3) 第 i個非 0元素值的實部。 1,3,7]。 xxx( 0 ):1 Bx? ?注 意 ( ) 時 不 能 說 “ 對 任 意 都 不 收 斂 ” 。%構(gòu)造對角陣 D D1=inv(D)。%求出迭代矩陣 f=D1*b。 s=B*x0+f。1 1 15]。 x0=[0。 s=jacobi(A,b,x0,err) 結(jié)果: n 1 2 3 4 5 6 7 例 三、 GaussSeidel 迭代法 對于 Ax = b 22( 1 ) ( ) ( ) ( )1 1 1 1 2 2 111( 1 ) ( ) ( ) ( )2 2 2 1 1 2 2( 1 ) ( ) ( ) ( )1 1 2 21( 0 ...... .... )1( 0 ...... .... )1( ...... .... 0 )k k k knnk k k knnk k k kn n n n nnnx b x a x a xax b a x x a xax b a x a x xa????? ? ? ? ????? ? ? ? ??????? ? ? ? ???L L L L L L L L L L L L LJacobi迭代公式為 (k=0,1,2,…) (4) ( 1 ) ( )1nkki i i j j i ij j ix b a x a????????????( 1 ) ( ) ( ) ( )1 1 1 1 2 2 111( 1 ) ( 1 ) ( ) ( )2 2 2 1 1 2 222( 1 ) ( 1 ) ( 1 ) ( )1 1 2 21 01 0 1 0k k k knnk k k knnk k k kn n n n nnnx b x a x a xax b a x x a xax b a x a x xa???? ? ????? ? ? ? ???????? ? ? ? ??????? ? ? ? ??????????與 ( 4) 對應(yīng)的迭代公式為: —— GaussSeidel迭代法 1( 1 ) ( 1 ) ( )11ink k ki i i j j i j j i ij j ix b a x a x a???? ? ???? ? ???????(5) 法 1: 法 2: Ax = b A=D+L+U (D+L) x = U x + b = 迭代公式: x(k+1) = BGS x(k)+ f , (k=0,1,2…) 11( ) ( )x D L U x D L b??? ? ? ? ? ?( 1 ) 1 ( ) 1( ) ( )kkx D L U x D L b? ? ?? ? ? ? ? ?11( ) , ( )GSB D L U f D L b??? ? ? ? ? ? ?L U D 編寫實現(xiàn) Seidel迭代法的函數(shù) : function s=seidel(a,b,x0,err) % seidel迭代法求解線性方程組, a為系數(shù)矩陣, b為%ax=b右邊的 %矩陣 b,x0為迭代初值, err為迭代誤差 if nargin==3 err=。%構(gòu)造嚴格上三角陣 C=inv(D+L)。 s=B*x0+f。1 1 15]。 x0=[0。 s=seidel(A,b,x0,err) 結(jié)果: n 1 2 2 3 4 例 發(fā)現(xiàn): seidel收斂速度比 jacobi快 — 松弛因子, =1 即 Seidel方法 (5) (6)是一種加權(quán)平均。%構(gòu)造對角陣 D L=tril(a,1)。 f=w*C*b。 end n 1 2 31 2 31 2 39710 815 13x x xx x xx x x? ? ???? ? ? ??? ? ? ? ??A=[9 1 1。8。0]。 不可解問題 線性方程組并不總是數(shù)值可解的,考慮如下三 個方程組。 1 , 1 ] \ [ 1 。 1 。 MATLAB中計算條件數(shù)的命令是: cond(A) 對于病態(tài)矩陣,逆矩陣和行列式的計算都會變得不精確。 b=[1 2 2