[理學]第三章matlab線性方程組-文庫吧在線文庫

　　

【正文】 ??Ax b? (0)x01???02???(0)x(0)x( 1 ) ( 1 ) 1 ( 1 ) ( )( 1 ) ( )k k k kx x D b L x U x??? ? ? ?? ? ? ? ?法 2： SOR法的矩陣迭代格式： ( 1 ) ( )kkx B x f??? ? ?11( ) [ ( 1 ) ] , ( )B D L D U f D L b? ? ? ? ? ???? ? ? ? ? ?1( 1 ) ( ) ( 1 ) ( )111( 1 )ink k k ki i i i j j i j jj j iiix x b a x a xa?????? ? ???? ? ? ? ? ???????編寫實現 SOR迭代法的函數 : function s=SOR(a,b,x0,w,err) % seidel迭代法求解線性方程組， a為系數矩陣， b為%ax=b右邊的 %矩陣 b,x0為迭代初值， w為松弛因子，%err為迭代誤差 if nargin==4 err=。 s=B*x0+f。 x0=[0。 2 , 2 ] \ [ 1 。 病態問題 有許多線性方程組理論上是可解的，但實際計算中由于受到舍入誤差的干擾而無法得到精確解，此類問題稱為 病態問題 。 A\bA\b1 for k=3:8 H=hilb(k)。( 3 ) * ( )4 ( * ( ) ) 1 * ( )AAinv inv A AA inv Ac o n d A inv A A inv A?（ ） 的 計 算 值 與 有 較 大 偏 離 ；（ ） 的 計 算 值 明 顯 不 等 于的 計 算 值 與 單 位 矩 陣 有 較 大 偏 離 ；（ ） 與 的 偏 離 可 作 為與 單 位 陣 的 偏 差 的 精 確 度 量 。 1 , 1 。 w=。1 1 15]。%構造嚴格上三角陣 C=inv(D+w*L)。0]。 end n 1 2 31 2 31 2 39710 815 13x x xx x xx x x? ? ???? ? ? ??? ? ? ? ??A=[9 1 1。%構造對角陣 D L=tril(a,1)。8。 n=1。 定理 定理 對于簡單迭代法，若迭代矩陣 ( 1 ) ( ) ( 0 ) 1 kkx B x g xB?? ????對 任 意 初 始 向 量 都 收 斂（ ）設有方程組 ( 其中 ) Ax = b， 即 0iia ?11 1 12 2 1 121 1 22 2 2 21 1 2 2nnnnn n nn n na x a x a x ba x a x a x ba x a x a x b? ? ? ???? ? ? ??????? ? ? ???LLLLLLLLLLLLLLLL(2) 作等價變形 1 1 1 12 2 1112 2 21 1 2 2221 1 2 21()1()1()00 0nnnnn n n n nnnx b x a x a xax b a x x a xax b a x a x xa?? ? ? ? ????? ? ? ? ???????? ? ? ? ????(3) Jacobi迭代法 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 L(k=0,1,2,…) (4) ( 1 ) ( )1nkki i i j j i ij j ix b a x a????????????法 1： 法 2： Ax = b A=D+(L+U) Dx = (L+U ) x + b = x = D1(L+U)x + D1 b L U = x = BJ x+ f = 迭代公式： x(k+1) = BJ x(k)+ f , (k=0,1,2…) BJ= D1（ L+U） , f =D1b D 編寫實現 Jacobo迭代法的函數 : function s=jacobi(a,b,x0,err) % jacobi迭代法求解線性方程組， a為系數矩陣， b為%ax=b右邊的 %矩陣 b,x0為迭代初值， err為迭代誤差 if nargin==3 err=。 A=[1,1,4。 S是要建立的稀疏矩陣的非 0元素； u(i)、 v(i)分別是 S(i)的行標和列標； 該函數生成一個 max(u)行、 max(v)列并以 S為稀疏元素的稀疏矩陣 u=[1,1,4]。1 1 1 1,d d b b??Step2： 1 2 3() nx x x x? ? ?下面開始“趕”： Step3：先求 ( 1 , , 2 , 1 , )in??Step4：再求其他的 :nxnnndxb???1i i iiid c xxb?? ???1 , 1 ,nxx?三對角方程組的求解程序 如下： % tri_diag(a,b,c,d,n) solves a tridiagonal equation. function x = tri_diag(a,b,c,d,n) for i=2:n r=a(i)/b(i1)。1 1 1 1,d d b b??39。1， 3， 0， 2。 x=R\(Q\b) 或采用 QR分解的第 2種格式，命令如下： [Q， R， E]=qr(A)。2 6 9 18。 定理 n 階陣 有唯一 Doolittle分解 (Crout） 的所有 順序主子式不為 0. 三角分解不唯一 ,為此引入 定義 五．利用矩陣三角分解法求解方程組 為什么要討論三角分解？若在消元法進行前能實現三角分解 ， 則 容易回代求解 在 Gauss消去法中，選主元改變了行的次序，盡管對于 Gauss消去法來說，這種次序的變換無法事先知道，但是這個變化的影響卻可以用一個算子 P表示，其中 P是一個置換矩陣。 ]。 a % 第二個主元標準化 a = 0 for i=1:3 if i~=2, a(i,:)=a(i,:)a(i,2)*a(2,:)。 a(1,:)=tempo。 Gauss－ Jordan消去法的優點之一是用它來算逆矩陣的算法非常容易解釋 。 a a = 0 tempo = a(3,:)。 （ ) 例 1：在 MATLAB上，用 Gauss列主元消去法求解方程組： 123 4 4 2 3 6 6 2 1 4 4 8 0xxx? ??? ? ? ???? ? ? ??? ??? ? ? ?? ? ? ?????? ? ? ???clear a = [ 3。 [x,y]=line_solution(A,b) 分別顯示其求解結果 。4,3,1,2。 x=A\b。)。/ 39。2 0 5]。 r=rank(D) t=trace(A) 結果： r = 4 t= 6 九、求矩陣特征多項式、特征值、特征向量 設 A是一個 n n 矩陣， ( ) d e t( )f I A????為 A的特征多項式。 x1=det(D1)/DD。6。 B=pinv(A) 例 求矩陣 A的偽逆 A=[0,0,0。 A=sym(A) inv(A) 結果： ??? Error using == Error, (in inverse) singular matrix 矩陣的偽逆 pinv(A) 若 A不是一個方陣，或 A是一個非滿秩的方陣時，矩陣A沒有逆矩陣，但可以找到一個與 A的轉置矩陣 A’同型的矩陣 B，使得： A norm(double(H*inv(H))eye(size(H))) 說明 3： 對于奇異 矩陣用數值方法的 inv( )函數，會給出錯誤的結果，但符號工具箱中的 inv( )會給出正確的結論 例 A=[16,2,3,13。 x=inv(A)*b 一般情況下，用左除 x=A\b比求矩陣的逆的方法更有效 (因為 A奇異或接近奇異時，用 inv(A)可能出錯 ) 例： Hilbert矩陣（非常有名的病態矩陣）： 驗證從 5 5到 14 14的 Hilbert矩陣病態特征 clear format long。4,0,13,0,22。由楊輝三角形表組成的矩陣稱為帕斯卡 (Pascal)矩陣。 例 A=vander([1。對于 n階魔方陣，其元素由 1,2,3,…,n 2共 n2 個整數組成 . magic(n) 生成一個 n階魔方陣 ,各行各列及兩條 對角線和為 (1+2+3+...+n2 )/n 例 magic(5) ans = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 例：將 101~125等 25個數填入一個 5行 5列的表格中，使其每行每列及對角線的和均為 565。這里 x, y均為向量，兩者不必等長。 命令如下： V=[1,1/2,1]。1,4,9。 fprintf(39。 5,11,10,8。 例 求 A的偽逆，并將結果送 B A=[3,1,1,1。8,5,3,4。 %用方程組的右端向量置換 D的第 2列 D3=[D(:,1:2),b,D(:,4:4)]。4,3,1,2。 ?[V,D]=eig(A,B) 計算廣義特征值和特征向量，使 AV=BVD。[ 1 , 1 ] 39。 if norm(b)0 % b非零， 非齊次方程組 if rank(A)==rank([A,b]) %方程組相容 （ 有解 ） if rank(A)==n %有 唯一解 x=A\b。 end else %方程組不相容 （無解） ，給出最小二乘解 disp(39。r39。3,5,2,2。11 1 12 1 1 , 1222 2 2 2 , 1 2 2 , 1