實驗題目：
1、高斯消元法 Computer Problems P100 2.2，（a）（b）
(1)用MATLAB語言提供的方法解方程組（x=A\b）
(2)寫出直接LU分解函式（參考例2.13），給出A的直接LU分解結果
(3)寫前代、回代函式。結合LU分解，前代、回代解（a）（b）方程組。
(4)直接使用軟體環境提供的函式LU分解函式、解方程組方法比較
說明：(2)、(3)，這幾部分不得直接使用軟體環境提供的函式完成，可以參考課本上的演算法流程實現

實驗程式碼：

%實驗二第一問 解方程組（x=A\b）
A = [2,4,-2;4,9,-3;-2,-1,7]; 
b
 
 = [2;8;10];    
c = [4;8;-6];
x = A\b;
y = A\c;
fprintf('結果為：');
x
y

mylu函式：

function [ L,U,p ] = mylu( A )
%UNTITLED2 此處顯示有關此函式的摘要
%   此處顯示詳細說明
% LU 三角分解
%[L,U,p] = mylu(A) 
%生成單位下三角矩陣L，上三角矩陣U，以及排列向量p
%滿足：L*U = A(p,:)
%本函式可以解決主元位置上有0的情況

%獲得A的行數和列數
[n,n] = size(A);

%定義排列向量p
p = (1:n)';

%for迴圈，用 k 記錄消元步數
for k = 1:n-1

    %
 
為了排除第k列全為0以及主元位置為0
    %先找出第k列的對角元及以下元素的絕對值的最大值
    [r,m] = max(abs(A(k:n,k))); % r為最大值，m 為相對行數
    m = m+k-1;  %此時 m 為最大值在A中對應的 行數

    %如果第k列全為0 ，則跳過消元過程
    if (A([m,k]) ~= 0)

    %避免主元為0，將絕對值最大值交換到主元所在行,現在的主元是對角元
    if( m ~= k)
        A([k m],:) = A([m k],:);  %矩陣A兩行進行交換
        p([k m]) = p([m k]) ;      %排列向量p 兩行進行交換
    end

    %
 
計算第k列元素除以主元后得到的數
     i = k+1:n;
     A(i,k) = A(i,k)/A(k,k);

     %更新矩陣剩餘部分
     j = k+1:n;
     A(i,j) = A(i,j) - A(i,k)*A(k,j);
    end
end

%分離結果
L = tril(A,-1) + eye(n,n); %tril()函式提取矩陣的下三角部分，eye()獲得單位矩陣
U = triu(A);               %tril()函式提取矩陣的上三角部分

end

利用mylu函式程式碼如下：

%實驗二的（2）寫出直接LU分解函式，給出A的直接LU分解結果
A = [2,4,-2;4,9,-3;-2,-1,7]; 

[L,U,p] = my_lu(A);
fprintf('矩陣A的 L*U 分解為：');
L
U
fprintf(' L*U 形成矩陣 LU 為：');
LU = L*U 
fprintf(' 排列向量 p 為：');
p'
fprintf('矩陣 A 為：');
A

（3）兩個函式，前代函式以及後代函式：

function x = myforward(L,x)

%my_forward  前向消元，前代
%對於下三角矩陣L， x = my_forward(L,b)給出 L*x = b 的解x

[n,n] = size(L);

%消元過程
for k = 1:n
    j = 1:k-1 ;
    x(k)= x(k) - L(k,j)*x(j);
    x(k)= x(k)/L(k,k);
end
end

function x = myback(U,x)

% my_back  後向消元
%對於上三角矩陣U，  x = my_back(U,y) 給出 U*x = y 的解x

[n,n] = size(U);

%消元過程
for k = n:-1:1
    j = k+1:n ;
    x(k) = (x(k) - U(k,j)*x(j))/U(k,k);

end
end

使用自編函式求解程式碼：

%第三問 寫前代、回代函式。結合LU分解，前代、回代解（a）（b）方程組。
A = [2,4,-2;4,9,-3;-2,-1,7]; 
b = [2;8;10];    
c = [4;8;-6];
[L,U,p] = mylu(A);

%重新排列 b ，向前消元，L*U*x1 = b，其中 L*y1 = b
y1 = myforward(L,b(p));
%反向回代
x1 = myback(U,y1);
fprintf('利用前代、回代函式，結合LU分解，得到的答案 x 為：')
x = x1

%重新排列 c ，向前消元，L*U*x2 = c，其中 L*y2 = c
y2 = myforward(L,c(p));


%反向回代
x2 = myback(U,y2);
fprintf('利用前代、回代函式，結合LU分解，得到的答案 y 為：')
y = x2

%實驗二的第四問 直接使用軟體環境提供的函式LU分解函式
[L,U,P] = lu(A);
fprintf('得到的 L*U 分解為：');
L
U