% exp3_2.m --- 矩陣分解命令的學習

function study_lu_and_chol

% ---------- LU分解 ----------
% [簡介] 設A是非奇異矩陣,則有如下列主元的LU分解
%             PA = LU
%        其中P是置換矩陣(又稱排列矩陣),L是單位下三角矩陣且|l(ij)|<=1,U是上三角矩陣
%        矩陣A如果有了上述分解則求解線性方程組 Ax = b 就等價于分別求解兩個三角方程組
%            Ly = Pb  和  Ux=y
%        而求解三角方程組是非常容易的.用這種方法求解方程組與列主元的Gauss消去法是等價的
%        但 LU 分解還有其它優點. 參見 P65 例9 和 P70 實驗課題(一)

A = [10  -7  0
     -3   2  6
      5  -1  5];
[L U P] = lu(A)

% 如果再給b,則下面是解兩個三角方程組的程序. 參見P49(3-15)式
b = [7 4 6]';
[n n]=size(A);
x = zeros(n,1);

% 前代求解:Ly = Pb(解用x儲存)
b = P*b;
x(1) = b(1);
for k = 2:n
    x(k) = b(k)-L(k,1:k-1)*x(1:k-1);
end

% 回代求解:Ux = y (這里先y=x,解仍用x儲存)
x(n) = x(n)/U(n,n);
for k = n-1:-1:1
    x(k) = ( x(k)-U(k,k+1:n)*x(k+1:n) ) / U(k,k);
end
x

% ---------- Cholesky分解 ----------
% [簡介] 設A是(實對稱)正定矩陣,則有如下唯一的分解
%             A = R'*R
%        其中 R 是上三角矩陣且主對角元全為正

% 例
A = [4    -1     1
    -1  4.25  2.75
     1  2.75   3.5];
R = chol(A)

% [注1]  有了Chol分解,求解 Ax = b 又可轉化為分別求解兩個三角方程組
%        這是容易的,請同學們自己完成(常稱為平方根法)
% [注2]  P51 的 chol 分解 A = LDL'
%        避免了前面分解中的開方運算(開方運算比四則運算速度是慢的)



% ******* 你的實驗  *******

% ★【實驗一】
% 首先完成下面Gauss列選主元的消去法程序,單獨存為 gauss.m 文件(注意一定要與函數名相同)
% 可參考P43 圖3-2
% 然后找一個例子調用此程序驗證是否正確(調用方法同Matlab內部函數調用完全一樣)

% function x = gauss(A,b)
% [n,n] = size(A);
% x = zeros(n,1);
% 
% Aug = [A,b];                         % 增廣矩陣
% 
% for k = 1:n-1
%     [piv,r] = max(abs(Aug(k:n,k)));  % 找列主元所在子矩陣的行r
%     r = r + k - 1;                   % 列主元所在大矩陣的行
%     
%     if r>k
%         ?                            % 對Aug實施行交換(一行命令就可以了,怎么寫?)
%     end
%     
%     if Aug(k,k)==0, error('?'), end  % 程序遇到error會中斷執行并顯示其中的提示內容
%     
%     % 把增廣矩陣消元成為上三角
%     for p = k+1:n
%         mult = Aug(p,k)/Aug(k,k);    % 消元乘子
%         Aug(p,k:n+1) = ?;
%     end
% end 
% 
% % 解上三角方程組
% A = Aug(:,1:n); b = Aug(:,n+1);
% x(n) = ?;
% for k = n-1:-1:1
%     x(k) = ?;
% end



% 【實驗二】編下面程序
%  (1) 追趕法(P45) 
%  (2) Cholesky分解法(P51)

%  下面是追趕法的程序你可以參考( 參見 P45 ),注意向量 a 的下標與書上不同

% function d = tridiag(a,b,c,d)
% % d = tridiag(a,b,c,d) --- 求解三對角線性方程組的追趕法( 參見 P45 )
% % 方程組的形式
% %         |b1  c1            | |x1|     |d1|
% %         |a1  b2  c2        | |x2|     |d2|
% %         |   a2   b3  c3    | |x3|  =  |d3|
% %         |        a3  b4  c4| |x4|     |d4|
% %         |            a4  b5| |x5|     |d5|
% % 輸入: a,b,c,d 是四個向量 
% % 輸出: d 方程組的解

% n = length(d);
% k = 1;
% d(k) = d(k)/b(k);
% c(k) = c(k)/b(k);
% for k = 2:n-1
%     b(k) = b(k) - a(k-1)*c(k-1);
% 	  c(k) = c(k)/b(k);
%     d(k) = ( d(k) - a(k-1)*d(k-1) )/b(k);
% end
% k = n;
% b(k) = b(k)-a(k-1)*c(k-1);
% d(k) = ( d(k) - a(k-1)*d(k-1) )/b(k);
% 
% % d(n) = d(n);
% for k = n-1:-1:1 
% 	d(k) = d(k) - c(k)*d(k+1);
% end