power法一次只能求一個特徵值,Rutishauser在1958年提出\(LR\)演算法,可以求出全部的特徵值。
\(LR\)演算法是基於\(LU\)分解,也就是將\(A\)分解成\(L\)和\(U\)兩個矩陣,其中\(L\)為下三角矩陣,\(U\)為上三角矩陣,再將\(U\)乘上\(L\)形成下一個矩陣\(A\),如此反覆計算最後的矩陣\(A\)的對角線就是全部的特徵值。
關於\(LU\)分解請見
https://ccjou.wordpress.com/2010/09/01/lu-分解/
---------------------------
虛擬碼
1.對於\(k=0,1,2,\ldots\),直到\( \lambda(1),\lambda(2),\ldots,\lambda(n) \)收斂為止,計算
2.\(A_k=L_k U_k\),計算\(A_k\)的\(LU\)分解。
3.\(A_{k+1}=U_k L_k\),將\(L_k\)和\(U_k\)對調後相乘得到下一個矩陣\(A_{k+1}\)
4.特徵值\( \lambda(1),\lambda(2),\ldots,\lambda(n) \)在矩陣\(A_{k+1}\)的對角線上。
---------------------------
計算\(A=\left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right] \)全部的特徵值。
[計算過程]
第1輪:
計算矩陣\(A\)的\(LU\)分解
\( \left[ \matrix{-4&0&1 \cr -1&3&1 \cr 1&2&5} \right]=\left[ \matrix{1&0&0 \cr 0.25&1&0 \cr -0.25&0.6667&1} \right] \left[ \matrix{-4&0&1 \cr 0&3.0&0.75 \cr 0&0&4.75} \right] \)
將\(L\)和\(U\)對調後相乘得到下一個矩陣\(A\)
\( \left[ \matrix{-4&0&1 \cr 0&3.0&0.75 \cr 0&0&4.75} \right] \left[ \matrix{1&0&0 \cr 0.25&1&0 \cr -0.25&0.6667&1} \right]=\left[ \matrix{-4.25&0.6667&1 \cr 0.5625&3.5&0.75 \cr -1.1875&3.1667&4.75} \right] \)
特徵值\( \left[-4.25,3.5,4.75 \right] \)
第2輪:
計算矩陣\(A\)的\(LU\)分解
\( \left[ \matrix{-4.25&0.6667&1 \cr 0.5625&3.5&0.75 \cr -1.1875&3.1667&4.75} \right]=\left[ \matrix{1&0&0 \cr -0.1324&1&0 \cr 0.2794&0.8306&1} \right] \left[ \matrix{-4.25&0.6667&1 \cr 0&3.5882&0.8824 \cr 0&0&3.7377} \right] \)
將\(L\)和\(U\)對調後相乘得到下一個矩陣\(A\)
\( \left[ \matrix{-4.25&0.6667&1 \cr 0&3.5882&0.8824 \cr 0&0&3.7377} \right] \left[ \matrix{1&0&0 \cr -0.1324&1&0 \cr 0.2794&0.8306&1} \right]=\left[ \matrix{-4.0588&1.4973&1.0 \cr -0.2284&4.3211&0.8824 \cr 1.0444&3.1045&3.7377} \right] \)
特徵值\( \left[ -4.0588,4.3211,3.7377 \right] \)
第3輪:
計算矩陣\(A\)的\(LU\)分解
\( \left[ \matrix{-4.0588&1.4973&1.0 \cr -0.2284&4.3211&0.8824 \cr 1.0444&3.1045&3.7377} \right]=\left[ \matrix{1&0&0 \cr 0.0563&1&0 \cr -0.2573&0.8237&1} \right] \left[ \matrix{-4.0588&1.4973&1.0 \cr 0&4.2369&0.8261 \cr 0&0&3.3146} \right] \)
將\(L\)和\(U\)對調後相乘得到下一個矩陣\(A\)
\( \left[ \matrix{-4.0588&1.4973&1.0 \cr 0&4.2369&0.8261 \cr 0&0&3.3146} \right] \left[ \matrix{1&0&0 \cr 0.0563&1&0 \cr -0.2573&0.8237&1} \right]=\left[ \matrix{-4.2319&2.3209&1.0 \cr 0.0253&4.9173&0.8261 \cr -0.8529&2.7301&3.3145} \right] \)
特徵值\( \left[ -4.2319 , 4.9173 , 3.3146 \right] \)
反覆計算直到全部特徵值\( \left[ 5.7542,-4.1444,2.3902 \right] \)
設定為數值運算
(%i1) numer:true;
(%o1) true
不顯示小數轉換成分數的過程
例如rat: replaced -4.25 by -17/4 = -4.25
(%i2) ratprint:false;
(%o2) false
LR迭代法副程式
/*print(...);*/將兩邊的/*和*/刪除後可以顯示計算過程
但顯示太多計算過程maxima程式會當掉,只好註解起來不顯示
(%i3)
LR(A,epsilon):=block
([count:0,P,L,U,n:length(A),eig1,eig2],
eig2:makelist(0,n),
do
(count:count+1,
/*print("第",count,"輪"),*/
[P,L,U]:get_lu_factors(lu_factor(A)),
/*print("計算矩陣A的LU分解",A,"="),*/
/*print(L,"*",U),*/
A:U.L,
/*print("將L和U對調後相乘得到下一個矩陣A"),*/
/*print(U,"*",L,"=",A),*/
eig1:create_list(A[i,i],i,1,n),
print("特徵值",eig1),
if lmax(abs(eig1-eig2))<epsilon then
(return(eig1)),
eig2:eig1
)
)$
矩陣A
(%i4)
A:matrix([-4,0,1],
[-1,3,1],
[1,2,5]);
(%o4) \( \left[ \matrix{-4&0&1\cr -1&3&1\cr 1&2&5} \right] \)
計算矩陣A全部的特徵值,誤差10^-5
(%i5) LR(A,10^-5);
0 errors, 0 warnings
特徵值\( \left[-4.25,3.5,4.75\right] \)
特徵值\( \left[-4.058823529411764,4.321118611378978,3.737704918032787\right] \)
特徵值\( \left[-4.231884057971015,4.91729873861925,3.314585319351763\right] \)
特徵值\( \left[-4.044520547945205,5.313248850521489,2.731271697423716\right] \)
特徵值\( \left[-4.27519051651143,5.669125227250055,2.606065289261376\right] \)
特徵值\( \left[-3.97464844523668,5.528967112469076,2.445681332767604\right] \)
特徵值\( \left[-4.386336456049431,5.951388036215438,2.434948419833993\right] \)
特徵值\( \left[-3.82905992615734,5.432559904220157,2.396500021937184\right] \)
特徵值\( \left[-4.616971218176369,6.21703651117396,2.39993470700241\right] \)
特徵值\( \left[-3.556070568397674,5.165962621376724,2.39010794702095\right] \)
特徵值\( \left[-5.09592634683869,6.703414031205933,2.392512315632756\right] \)
特徵值\( \left[-3.07014305149453,4.680374938114763,2.389768113379766\right] \)
特徵值\( \left[-6.15763906410022,7.766852181380823,2.390786882719396\right] \)
特徵值\( \left[-2.263100148442174,3.873135110765942,2.389965037676232\right] \)
特徵值\( \left[-8.927750729051786,10.53740902683053,2.390341702221251\right] \)
特徵值\( \left[-1.061374420236201,2.671290372000583,2.390084048235619\right] \)
特徵值\( \left[-20.85888707850496,22.46867122537871,2.390215853126253\right] \)
特徵值\( \left[0.4665444106715815,1.143322725027585,2.390132864300833\right] \)
特徵值\( \left[52.72558253627871,-51.11576036528454,2.390177829005824\right] \)
特徵值\( \left[2.06213873325895,-0.4522894269100236,2.390150693651076\right] \)
特徵值\( \left[13.17442195234688,-11.56458779576555,2.390165843418669\right] \)
特徵值\( \left[3.419996923785756,-1.810153822408377,2.390156898622622\right] \)
特徵值\( \left[8.58287969592947,-6.97304166698185,2.39016197105238\right] \)
特徵值\( \left[4.388368293639632,-2.77852730349426,2.39015900985463\right] \)
特徵值\( \left[7.044154778522382,-5.434315481286824,2.390160702764444\right] \)
特徵值\( \left[4.995309976952969,-3.385469697199792,2.390159720246825\right] \)
特徵值\( \left[6.383873007186085,-4.774033291505344,2.390160284319262\right] \)
特徵值\( \left[5.34546797814534,-3.735627936085535,2.390159957940198\right] \)
特徵值\( \left[6.071147185395128,-4.46130733112236,2.390160145727236\right] \)
特徵值\( \left[5.537890814471758,-3.928050851714158,2.390160037242404\right] \)
特徵值\( \left[5.916131959742552,-4.306292059474105,2.390160099731557\right] \)
特徵值\( \left[5.640814092845454,-4.030974156506412,2.390160063660963\right] \)
特徵值\( \left[5.8375586037984,-4.227718688248872,2.390160084450476\right] \)
特徵值\( \left[5.695070878935839,-4.08523095139107,2.390160072455235\right] \)
特徵值\( \left[5.797281242393553,-4.187441321764426,2.390160079370877\right] \)
特徵值\( \left[5.723453549499145,-4.113613624880672,2.390160075381531\right] \)
特徵值\( \left[5.776515673768994,-4.166675751450871,2.390160077681882\right] \)
特徵值\( \left[5.738241274640624,-4.128401350995675,2.390160076355056\right] \)
特徵值\( \left[5.765777948950188,-4.15593802607038,2.390160077120198\right] \)
特徵值\( \left[5.745929679925635,-4.136089756604526,2.390160076678896\right] \)
特徵值\( \left[5.760217049729156,-4.150377126662543,2.390160076933392\right] \)
特徵值\( \left[5.749922650695337,-4.140082727481946,2.390160076786614\right] \)
特徵值\( \left[5.757334865818812,-4.14749494269007,2.390160076871262\right] \)
特徵值\( \left[5.751995220747789,-4.142155297570226,2.390160076822442\right] \)
特徵值\( \left[5.7558404324125,-4.146000509263092,2.390160076850597\right] \)
特徵值\( \left[5.753070680591324,-4.143230757425678,2.39016007683436\right] \)
特徵值\( \left[5.755065392906378,-4.145225469750097,2.390160076843724\right] \)
特徵值\( \left[5.753628652968832,-4.14378872980715,2.390160076838324\right] \)
特徵值\( \left[5.754663399407986,-4.144823476249418,2.390160076841438\right] \)
特徵值\( \left[5.753918118401036,-4.144078195240673,2.390160076839642\right] \)
特徵值\( \left[5.754454883559434,-4.144614960400106,2.390160076840678\right] \)
特徵值\( \left[5.75406828136153,-4.144228358201605,2.390160076840081\right] \)
特徵值\( \left[5.754346722233892,-4.144506799074311,2.390160076840425\right] \)
特徵值\( \left[5.754146178167928,-4.144306255008148,2.390160076840227\right] \)
特徵值\( \left[5.754290615935278,-4.144450692775613,2.390160076840341\right] \)
特徵值\( \left[5.754186586569208,-4.144346663409477,2.390160076840275\right] \)
特徵值\( \left[5.754261511799679,-4.144421588639986,2.390160076840314\right] \)
特徵值\( \left[5.754207548012657,-4.144367624852942,2.390160076840291\right] \)
特徵值\( \left[5.754246414489054,-4.144406491329351,2.390160076840304\right] \)
特徵值\( \left[5.75421842151387,-4.144378498354161,2.390160076840297\right] \)
特徵值\( \left[5.754238582980808,-4.144398659821102,2.390160076840301\right] \)
特徵值\( \left[5.75422406200547,-4.144384138845762,2.390160076840299\right] \)
特徵值\( \left[5.75423452049635,-4.144394597336643,2.3901600768403\right] \)
特徵值\( \left[5.75422698793699,-4.144387064777282,2.390160076840299\right] \)
(%o5) \( \left[5.75422698793699,-4.144387064777282,2.390160076840299\right] \)
參考資料
https://books.google.com.tw/book ... e&q&f=false
[
本帖最後由 bugmens 於 2016-9-16 10:41 PM 編輯 ]