Doolittle分解法（LU分解）詳細分析以及matlab的實現
一、基本介紹
前面介紹的Gauss消去法實際上做的事情是將係數矩陣A做了一個三角分解，即：
A=LU
式（1）
其中，L為單位下三角陣，U為上三角陣，該分解唯一。若A為非奇異，則U也非奇異。
實際消元過程如下所示：
第1步對應將係數矩陣
其中，有
消元的第2步對應為
其中，有
以此類推，第n-1步消元后，有：
則消去過程對應的矩陣變換為
由式（1.8）可匯出LU分解式（1），此處L=L1L2...Ln-1仍然是單位下三角陣，U為上三角陣，具體形式為：
在實際程式設計計算的時候，只需要最終的兩個矩陣L和U，此時，原方程Ax=b的求解就轉化為了兩個三角形方程組的求解
採用回代的方式即可求出中間y和結果x。
Doolittle分解可直接通過矩陣乘法匯出計算過程。設A=LU，即
由矩陣乘法及兩矩陣相等，可得
第一步，應l11=1,故u1j
= a1j(j=1,2,...,n)
一般地，若U的前i-1行及L的前i-1列元素已計算出來，則
第i步，由
得
又由
得
綜上可得，A的LU分解公式如下
對i = 1,2,...,n
二、演算法分析
在程式設計實現的時候，真正需要主要的公式並不多，很多看似複雜的推到，實際上只是舉例幾步，使得我們能夠舉一反三，更加清楚地理解推導式的由來。那麼在這個程式碼中，我們需要使用到的公式包括：（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分解法比前兩者效率更高。