Doolittle分解法(LU分解)詳細分析以及matlab的實現
一、基本介紹
前面介紹的Gauss消去法實際上做的事情是將係數矩陣A做了一個三角分解,即:
A=LU
式(1)
其中,L為單位下三角陣,U為上三角陣,該分解唯一。若A為非奇異,則U也非奇異。
實際消元過程如下所示:
第1步對應將係數矩陣和右端項左乘一個初等變換陣,即
式(1.1)
其中,有
式(1.2)
式(1.3)
消元的第2步對應為
式(1.4)
其中,有
式(1.5)
式(1.6)
以此類推,第n-1步消元后,有:
式(1.7)
則消去過程對應的矩陣變換為
式(1.8)
由式(1.8)可匯出LU分解式(1),此處L=L1L2...Ln-1仍然是單位下三角陣,U為上三角陣,具體形式為:
,
式(1.9)
在實際程式設計計算的時候,只需要最終的兩個矩陣L和U,此時,原方程Ax=b的求解就轉化為了兩個三角形方程組的求解
式(1.10)
採用回代的方式即可求出中間y和結果x。
式(1.11)
式(1.12)
Doolittle分解可直接通過矩陣乘法匯出計算過程。設A=LU,即
式(1.13)
由矩陣乘法及兩矩陣相等,可得
第一步,應l11=1,故u1j
= a1j(j=1,2,...,n)
一般地,若U的前i-1行及L的前i-1列元素已計算出來,則
第i步,由
式(1.14)
得
式(1.15)
又由
式(1.16)
得
式(1.17)
綜上可得,A的LU分解公式如下
對i = 1,2,...,n
式(1.18)
二、演算法分析
在程式設計實現的時候,真正需要主要的公式並不多,很多看似複雜的推到,實際上只是舉例幾步,使得我們能夠舉一反三,更加清楚地理解推導式的由來。那麼在這個程式碼中,我們需要使用到的公式包括:(1.11)、(1.12)和(1.18)。而其他公式如果無法理解其實並不影響程式碼的實現,只是為了更好地理解該知識點,還是建議大家自己動手推導一兩步,可能會花一些時間,但肯定都是值得的。
程式碼實現具體步驟:
第一步:
初始化u1i,其中i = 1,2,...,n。這裡初始化上三角陣的第一行的原因是在式(1.18)的累加求和中使用到i-1,當i=1時,對於uk0我們並不能找到這樣的一個值。同理也需要對下三角陣的第一列進行初始化。(如果哪位小夥伴有更好的方法,歡迎留言討論)
第二步:
採用式(1.18)遞推得到整個上三角陣和下三角陣。
第三步:
採用式(1.11)回代求解中間矩陣y。
第四步:
採用式(1.12)回代求解最終結果x。
注:或許有些小夥伴對於中間複雜公式怎麼求解,其實這只是一個迭代過程,只不過在迭代地過程中需要理解什麼時候需要使用迴圈,所使用的迴圈初始值是多少,結束是多少,是否需要中間變數。如果理解清楚這幾點,相信對於同類型的問題應該是沒有任何難度問題了。如果還是沒有理解清楚,歡迎私信我。
本文只給出了matlab語言的實現,不同語言實現的有一定的區別,歡迎大家使用更多不同語言來編寫程式。
三、程式碼實現
matlab程式碼實現如下:
function [L_matrix,U_matrix,y_matrix,x_matrix] = LU_separetion(A_matrix, B_matirx)
% LU係數矩陣分解
% 2017-11-09 xh_scu
% inputs:
% A_matrix:輸入的係數矩陣,尺寸為[n,n]
% B_matrix:輸入的乘積矩陣,尺寸為[n,1]
% outputs:
% L_matrix:下三角陣,尺寸為[n,n]
% U_matrix:上三角陣,尺寸為[n,n]
% y_matrix:中間矩陣,尺寸為[n,1]
% x_matrix:結果矩陣,尺寸為[n,1]
%% 第一步:初始化
% 獲取n值
[row_a, col_a] = size(A_matrix);
% 初始化上三角陣的第一行
for j = 1:col_a % for-1
U_matrix(1,j) = A_matrix(1,j);
end % for-1
% 初始化下三角陣的第一列
L_matrix(1,1) = 1;
for i = 2:row_a % for-2-s
L_matrix(i,1) = A_matrix(i,1)/A_matrix(1,1); % 對應式(1.3)
end % for-2-e
%% 第二步:前向分解計算
for i = 2:row_a % for-3-s
for j = i:col_a % for-4-s
temp_sum = 0;
for k = 1:i-1 % for-5-s
temp_sum = temp_sum + L_matrix(i,k)*U_matrix(k,j); %對應式(1.18)-上部分的求和部分
end % for-5-e
U_matrix(i,j) = A_matrix(i,j) - temp_sum; % 對應式(1.18)-上部分的求差部分
temp_sum_1 = 0;
for p = 1:i-1 % for-6-s
temp_sum_1 = temp_sum_1 + L_matrix(j,p)*U_matrix(p,i); % 對應式(1.18)-下部分的求和部分
end %for-6-e
L_matrix(j,i) = (A_matrix(j,i) - temp_sum_1)/U_matrix(i,i); % 對應式(1.18)-下部分的求差再求商部分
end % for-4-e
end % for-3-e
%% 第三步:回代計算y
x_matrix = zeros(row_a,1);
% 後向回代
% 下三角回代----計算中間矩陣Y
y_matrix(1,1) = B_matirx(1,1);
for i = 2:row_a % for-7-s
temp_sum_2 = 0;
for j = 1:i-1 % for-8-s
temp_sum_2 = temp_sum_2 + L_matrix(i,j)*y_matrix(j,1);
end % for-8-e
y_matrix(i,1) = B_matirx(i) - temp_sum_2;
end % for-7-e
%% 第四步:回代計算x
% 上三角回代----計算結果矩陣X
x_matrix(row_a,1) = y_matrix(row_a,1)/U_matrix(row_a,col_a);
for i=row_a-1:-1:1 % for-9-s
temp_sum_3 = 0;
for j = i+1:row_a % for-10-s
temp_sum_3 = temp_sum_3 + U_matrix(i,j)*x_matrix(j,1);
end % for-10-e
x_matrix(i,1) = (y_matrix(i,1) - temp_sum_3)/U_matrix(i,i);
end % for-9-e
end
四、測試分析
1)演算法的準確性測試
設輸入矩陣A = [2,2,3;4,7,7;-2,4,5], B
= [3,1,-7]
測試程式碼為:
A = [2,2,3;4,7,7;-2,4,5];
B = [3;1;-7];
[L,U,Y,X] = LU_separetion(A,B);
計算結果為:
L = [1,0,0;2,1,0;-1,2,1]
U = [2,2,3;0,3,1;0,0,6]
Y = [3;-5;6]
X = [2;-2;1]
與參考結果完全相等。
2)Gauss消去法、列主元素消去法以及LU分解法效能對比
設引數矩陣A為的元素為隨機數,取值範圍為[1,100],在相同輸入下測試各演算法的時間代價。
測試函式如下:
function [result] = test_function()
% 初始化結果矩陣[3,10]
result = zeros(3,10);
for i = 100:100:1000
% 產生隨機矩陣
A = randint(i,i,[1 100]);
B = randint(i,1,[1,100]);
% 分別呼叫三個函式
[~,time_1] = GaussElimination(A,B);
[~,time_2] = MainElementElimination(A,B);
[~,~,~,~,time_3] = LU_separetion(A,B);
% 將得到的計算時間結果送入結果矩陣
j = i/100;
result(1,j) = time_1;
result(2,j) = time_2;
result(3,j) = time_3;
end
end
注:測試函式只是進行示例說明,可能中間還存在不嚴謹的地方,如果有相關的意見或想法,可以留言一起探討。
測試結果(取四位小數(四捨五入))如下表所示:
維數 | 100 | 200 | 300 | 400 | 500 | 600 |
700 | 800 | 900 | 1000 |
Gauss消去法 | 0.0250 | 0.1910 | 0.6820 | 1.5680 | 3.2010 | 5.4100 | 8.4680 | 14.5070 | 21.4940 | 29.5890 |
列主元素消去法 | 0.0270 | 0.1970 | 0.6730 | 1.5860 | 3.2000 | 4.4970 | 8.7170 | 13.6610 | 21.0680 | 30.0600 |
LU分解法 | 0.0200 | 0.1480 | 0.5000 | 1.1750 | 2.5320 | 4.2160 | 7.0380 | 10.8440 | 14.1410 | 22.4300 |
測試結果圖例:
五、總結
1、本文分析了LU分解法的詳細實現,並對程式設計實現進行了主要步驟的說明,給出了matlab語言的實現程式碼。
2、測試了Gauss消去法、列主元素消去法以及LU分解法的計算效率,從測試結果可以得出:在相同的輸入情況下,LU分解法比前兩者效率更高。