直接三角分解的演算法

1. LU法求解線性方程組，matlab編程

線性方程組的三角分解求法其實和常用的高斯消去法等效。
如果要直接利用Matlab內置的三角分解演算法，可在命令窗口直接執行以下命令：
A=[1 4 0 1;1 5 1 0;-3 2 0 3;-4 0 1 4];
b=[11; 12; 7; 5];
[L,U]=lu(A); %L為下三角，U為上三角
x=U\(L\b) %求解x

若要自己編程實現以上演算法，可建立以下函數文件：

function x=GaussMethod(A,b)
%高斯消去法求解線性代數方程組Ax=b

n=size(A,1);
m=zeros(n-1,n-1);
x=zeros(n,1);

for k=1:n-1
for i=k+1:n
m(i,k)=-A(i,k)/A(k,k);
A(i,k:n)=A(i,k:n)+A(k,k:n)*m(i,k);
b(i)=b(i)+b(k)*m(i,k);
end
end

x(n)=b(n)/A(n,n);
for i=n-1:-1:1
p=0;
for j=i+1:n
p=p+A(i,j)*x(j);
end
x(i)=(b(i)-p)/A(i,i);
end

編寫函數後保存。在命令窗口輸入：
A=[1 4 0 1;1 5 1 0;-3 2 0 3;-4 0 1 4];
b=[11; 12; 7; 5];
x=GaussMethod(A,b)
運行後可得到
x =

1.0000
2.0000
1.0000
2.0000

2. 矩陣A 的LU 分解

將矩陣A分解為單位下三角陣L和上三角陣U，就是找出L的元素l_ij以及U的元素u_ij與A的元素a_ij的關系式。下面說明矩陣L，U的元素可以通過n步由矩陣A的元素計算出來，其中第k步求出U的第k行和L的第k列元素。

根據式（4－1）中，由矩陣乘法可知：

（1）由L的第1行和U的第j列元素對應相乘後與A的對應元素相等，即

a_1j＝u_1j

同理可知U的第1行元素：

u_1j＝a_1j （j＝1，2，…，n）

由L的第i行和U的第1列元素對應相乘後與A的對應元素相等，即

a_i1＝l_i1·u₁₁

從而得到L的第一列元素

l_i1＝a_i1/u₁₁ （i＝2，3，…，n）

這樣就確定了U的第1行元素及L的第1列元素。類似地，用矩陣乘法再確定U的第2行及L的第2列，如此繼續。

假設已經得出了U的前k－1行及L的前k－1列（1≤k≤n）的全部元素，現在來確定U的第k行元素和L的第k列元素。

（2）在式（4－1）中，由矩陣乘法知

地球物理數據處理基礎

注意，由於L是單位下三角陣，當r>k時，l_kr＝0，而l_kk＝1，所以

地球物理數據處理基礎

從而有

地球物理數據處理基礎

如此即可得到上三角矩陣U的第k行元素。

同理，可確定L的第k列元素。

由 並且當r>k時，u_rk＝0，則有

地球物理數據處理基礎

從而有

地球物理數據處理基礎

即可計算出矩陣L的第k列元素。

上述步驟進行，就可以定出L及U的全部元素，完成矩陣A的LU分解，即對k＝1，2，…，n，計算A＝LU分解的公式為

地球物理數據處理基礎

這里約定

上述這種矩陣A的LU分解計算順序也可按圖4－1所示逐步進行。

由於以上計算公式（4－7）中不含消去法的中間結果a^（k）_ij，可直接逐框計算，所以稱為緊湊格式。

從計算過程中可以看出，對A進行LU分解時，在求出u_1j（j＝1，2，…，n）後不再需要a_1j（j＝1，2，…，n），求出l_i1（i＝2，3，…，n）後不再需要a_i1（i＝2，3，…，n）。同樣求出u_kj和l_ik後不再需要a_kj和a_ik（j＝k，k＋1，…，n；i＝k＋1，k＋2，…，n；k＝1，2，…，n），因此，在計算過程中可以不另設存放L和U的數組，而是將u_kj和l_ik算出後直接存入矩陣A對應元素a_kj和a_ik的存貯單元中。這也稱為動態演算法，它可節省大量內存。

圖4－1 矩陣LU分解順序圖

對A進行LU分解後，可按式（4－5）和式（4－6）求解線性方程組。

同樣，求解線性方程組時，由式（4－5）和式（4－6）式可知，在求yi（i＝1，2，…，n）時，只用到對應的b_i（i＝1，2，…，n），故可不另設存放y的數組，而將y的元素存入列向量b對應元素b_i的存貯單元中。但應注意，這時的b已變成y。同理求x時，也不用另設存放單元，計算結果直接存入b中，因而在求解線性方程組的最後結果中，b存儲的是x。

另外，由於求yi的順代公式（4－5）與求u_kj的公式（4－7）形式完全相同，故可合並到k循環之內進行，這樣可以簡化程序結構。所以，在進行LU分解的同時，Ly＝b的求解也完成了。

［例1］用直接三角分解法求解線性方程組

解：（1）分解A＝LU，

地球物理數據處理基礎

即

地球物理數據處理基礎

（2）求解Ly＝b，

地球物理數據處理基礎

得到y＝（3，－5，6）^T

（3）求解

得到x＝（2，－2，1）^T

3. lu分解的改進

（i）Doolittle分解
對於非奇異矩陣（任n階順序主子式不全為0）的方陣A，都可以進行Doolittle分解，得到A=LU，其中L為單位下三角矩陣，U為上三角矩陣；這里的Doolittle分解實際就是Gauss變換；
（ii）Crout分解
對於非奇異矩陣（任n階順序主子式不全為0）的方陣A，都可以進行Crout分解，得到A=LU，其中L為下三角矩陣，U為單位上三角矩陣；
（iii）列主元三角分解
對於非奇異矩陣的方陣A，採用列主元三角分解，得到PA=LU，其中P為一個置換矩陣，L，U與Doolittle分解的規定相同；
（iv）全主元三角分解
對於非奇異矩陣的方陣A，採用全主元三角分解，得到PAQ=LU，其中P，Q為置換矩陣，L，U與Doolittle分解的規定相同；
（v）直接三角分解
對於非奇異矩陣的方陣A，利用直接三角分解推導得到的公式（Doolittle分解公式或者Crout分解公式），可以進行遞歸操作，以便於計算機編程實現；
（vi）「追趕法」
追趕法是針對帶狀矩陣（尤其是三對角矩陣）這一大稀疏矩陣的特殊結構，得出的一種保帶性分解的公式推導，實質結果也是LU分解；因為大稀疏矩陣在工程領域應用較多，所以這部分內容需要特別掌握。
（vii）Cholesky分解法（平方根法）和改進的平方根法
Cholesky分解法是是針對正定矩陣的分解，其結果是 A=LDLT=LD(1/2)D(1/2)LT=L1L1T。如何得到L1，實際也是給出了遞歸公式。
改進的平方根法是Cholesky分解的一種改進。為避免公式中開平方，得到的結果是A=LDLT=TLT， 同樣給出了求T,L的公式。
小結：
（1） 從（i）~（iv）是用手工計算的基礎方法，（v）~（vi）是用計算機輔助計算的演算法公式指導；
（2） 這些方法產生的目的是為了得到線性方程組的解，本質是高斯Gauss消元法 。

4. 急呀~~~跪求數學論文!題目是-病態方程組微分方程分解法

線性方程組的解法；矩陣特徵值與特徵向量的計算；非線性方程與非線性方程組的迭代解法；插值與逼近；數值積分；常微分方程初值問題的數值解法和偏微分方程的差分解法。內容豐富，系統性強，其深廣度適合工學碩士生的培養要求。本書語言簡練、流暢，數值例子和習題非常豐富。
商品信息

本書是為工學碩士研究生開設數值分析課而編寫的學位課教材。內容包括：線性方程組的解法；矩陣特徵值與特徵向量的計算；非線性方程與非線性方程組的迭代解法；插值與逼近；數值積分；常微分方程初值問題的數值解法和偏微分方程的差分解法。內容豐富，系統性強，其深廣度適合工學碩士生的培養要求。本書語言簡練、流暢，數值例子和習題非常豐富。
【目錄】
第一章 緒 論 1�
1.1 數值分析的研究對象 1�
1.2 誤差知識與演算法知識 1�
1.2.1 誤差 的來源與分類 1�
1.2.2 絕對誤差、相對誤差與有 效數字 3�
1.2.3 函數求值的誤差估計 5�
1.2.4 演算法及其計算復雜性 7�
1.3 向量范數與矩 陣范數 10�
1.3.1 向量范數 10�
1.3.2 矩陣范數 12�
習 題 18�
第二章 線性方程組 的解法 21�
2.1 Gauss消去法 22�
2.1.1 順序Gauss消去法 23�
2.1.2 列主元素Gauss消去法 25�
2.2 直接三角分解法 28�
2.2.1 Doolittle分解法與Crout分解法 28�
2.2.2 選主 元的Doolittle分解法 34�
2.2.3 三角分解法解帶狀 線性方程組 37�
2.2.4 追趕法求解三對角線性方程 組 41�
2.2.5 擬三對角線性方程組的求解方法 43 �
2.3 矩陣的條件數與病態線性方程組 45�
2.3.1 矩陣的條件數與線性方程組的性態 45�
2.3.2 關於病態線性方法組的求解問題 48�
2.4 迭代 法 51�
2.4.1 迭代法的一般形式及其收斂性 51 �
2.4.2 Jacobi迭代法 55�
2.4.3 Gauss�Seidel迭代法 60�
2.4.4 逐次超鬆弛迭 代法 64�
習 題 69�
第三章 矩陣特徵值與特 征向量的計算 74�
3.1 冪法和反冪法 74�
3.1.1 冪 法 74�
3.1.2 反冪法 79�
3.2 Jacobi方法 81�
3.3 QR方法 87�
3.3.1 矩陣的QR分解 87�
3.3.2 矩陣的擬上三角化 92�
3.3.3 帶雙步位移的QR方法 95�
習 題 100�
第四章 非線性方程與非線性方法組的迭代 解法 103�
4.1 非線性方程的迭代解法 103�
4.1.1 對分法 103�
4.1.2 簡單迭代法及其收斂 性 104�
4.1.3 簡單迭代法的收斂速度 109�
4.1.4 Steffensen加速收斂方法 112�
4.1.5 Newton法 115�
4.1.6 求方程m重根的 Newton法 120�
4.1.7 割線法 123�
4.1.8 單點割線法 127�
4.2 非線性方程組的迭代 解法 131�
4.2.1 一般概念 131�
4.2.2 簡單迭代法 134�
4.2.3 Newton法 138�
4.2.4 離散Newton法 140�
習 題 142�
第五章 插值與逼近 144�
5.1 代數插值 144�
5.1.1 一元函數插值 144�
5.1.2 二元函數插值 152�
5.2 Hermite插值 156�
5.3 樣條插 值 160�
5.3.1 樣條函數 160�
5.3.2 三次樣條插值問題 166�
5.3.3 B樣條為基底的三次樣 條插值函數 168�
5.3.4 三彎矩法求三次樣條插值 函數 172�
5.4 三角插值與快速Fourier變換 177�
5.4.1 周期函數的三角插值 177�
5.4.2 快速Fourier變換 180�
5.5 正交多項式 183�
5.5.1 正交多項式概念與性質 183�
5.5.2 幾種常 用的正交多項式 187�
5.6 函數的最佳平方逼近 193 �
5.6.1 最佳平方逼近的概念與解法 193�
5.6.2 正交函數系在最佳平方逼近中的應用 197�
5.6.3 樣條函數在最佳平方逼近中的應用 203�
5.6.4 離散型的最佳平方逼近 205�
5.6.5 曲線擬 合與曲面擬合 207�
習 題 219�
第六章 數值積 分 226�
6.1 求積公式及其代數精度 226�
6.2 插值型求積公式 228�
6.3 Newton�Cotes求積 公式 230�
6.4 Newton�Cotes求積公式的收斂性與數 值穩定性 236�
6.5 復化求積法 237�
6.5.1 復化梯形公式與復化Simpson公式 237�
6.5.2 區 間逐次分半法 242�
6.6 Romberg積分法 244�
6.6.1 Richardson外推技術 244�
6.6.2 Romberg 積分法 247�
6.7 Gauss型求積公式 249�
6.7.1 一般理論 249�
6.7.2 幾種Gauss型求積公式 255�
6.8 二重積分的數值求積法 263�
6.8.1 矩形域上的二重積分 263�
6.8.2 一般區 域上的二重積分 266�
習 題 267

5. 矩陣有哪幾種特殊分解

矩陣分解
(decomposition, factorization)是將矩陣拆解為數個矩陣的乘積，可分為三角分解、滿秩分解、QR分解、Jordan分解和SVD(奇異值)分解等，常見的有三種:1)三角分解法 (Triangular
Factorization)，2)QR
分解法
(QR
Factorization)，3)奇異值分解法
(Singular
Value
Decomposition)。
下面分別簡單介紹上面三個分解演算法：
1、三角分解法
三角分解法是將原正方
(square)
矩陣分解成一個上三角形矩陣或是排列(permuted)
的上三角形矩陣和一個
下三角形矩陣，這樣的分解法又稱為LU分解法。它的用途主要在簡化一個大矩陣的行列式值的計算過程，求逆矩陣，和求解聯立方程組。不過要注意這種分解法所得到的上下三角形矩陣並非唯一，還可找到數個不同
的一對上下三角形矩陣，此兩三角形矩陣相乘也會得到原矩陣。
MATLAB以lu函數來執行lu分解法，
其語法為[L,U]=lu(A)。
2、QR分解法
QR分解法是將矩陣分解成一個正規正交矩陣與上三角形矩陣,所以稱為QR分解法,與此正規正交矩陣的通用符號Q有關。
MATLAB以qr函數來執行QR分解法，
其語法為[Q,R]=qr(A)。
3、奇異值分解法
奇異值分解
(singular
value
decomposition,SVD)
是另一種正交矩陣分解法;SVD是最可靠的分解法，但是它比QR
分解法要花上近十倍的計算時間。[U,S,V]=svd(A)，其中U和V分別代表兩個正交矩陣，而S代表一對角矩陣。
和QR分解法相同，
原矩陣A不必為正方矩陣。使用SVD分解法的用途是解最小平方誤差法和數據壓縮。
MATLAB以svd函數來執行svd分解法，
其語法為[S,V,D]=svd(A)。

6. 實用數學演算法是什麼

常用數學演算法！也是計算方法程序-commonly used mathematical algorithms! Also calculation proceres
計算方法程序：
1.直接三角分解法

程序代碼：

#include <math.h>

#include <stdio.h>

main()

{

int i,j,n,r,k;

float s,t;

float a[100][100],b[100],U[100][100],L[100][100],x[100],y[100];

printf("in put n\n");

scanf("%d",&n);

printf("input a[*][*]\n");

for(i=0;i<n;i++)

for(j=0;j<n;j++)

scanf("%f",&a[i][j]);

printf("in put b[*]\n");

for(i=0;i<n;i++)

scanf("%f",&b[i]);

for(i=1;i<n;i++)

a[i][0]=a[i][0]/a[0][0];

­

for(k=1;k<n;k++)

{

for(j=k;j<n;j++)

{

s=0;

for (i=0;i<k;i++)

s=s+a[k][i]*a[i][j];

a[k][j]=a[k][j]-s;

}

for(i=k+1;i<n;i++)

{

t=0;

for(j=0;j<k;j++)

t=t+a[i][j]*a[j][k];

a[i][k]=(a[i][k]-t)/a[k][k];

}

}

­

for(i=0;i<n;i++)

for(j=0;j<n;j++)

{ if(i>j)

{ L[i][j]=a[i][j]; U[i][j]=0;}

else

{ U[i][j]=a[i][j];

if(i==j) L[i][j]=1;

else L[i][j]=0;

}

}

printf("\nL[*][*]=");

for(i=0;i<n;i++)

{ printf("\n");

for(j=0;j<n;j++)

printf(" %0.2f",L[i][j]);

}

­

printf("\nU[*][*]=");

for(i=0;i<n;i++)

{ printf("\n");

for(j=0;j<n;j++)

printf(" %0.2f",U[i][j]);

}

printf("\n");

y[0]=b[0];

for(k=1;k<n;k++)

{

s=0;

for (i=0;i<k;i++)

s=s+L[k][i]*y[i];

y[k]=b[k]-s;

}

x[n-1]=(y[n-1]/U[n-1][n-1]);

for(k=(n-2);k>=0;k--)

{

t=0;

for(i=k+1;i<n;i++)

t=t+U[k][i]*x[i];

x[k]=(y[k]-t)/U[k][k];

}

for(i=0;i<n;i++)

printf(" %0.2f",y[i]);

printf("\n");

printf("The answer is \n");

for(i=0;i<n;i++)

printf(" %0.2f",x[i]);

printf("\n");

}

/*不是完全接照實驗要求做的，沒有檢測這一步。只是照書上做的，如有需要，再加個IF語句判定*/

2.迭代法

程序代碼

#include <stdio.h>

#include <math.h>

double f(double x)

{

return 2*sin(x);

}

main()

{

int n=0;

double x,y,e;

printf("請輸入X的初值 X0:\n");

scanf("%lf",&x);

printf("請輸入精度e的值 e:\n");

scanf("%lf",&e);

do{

printf("x%d=%lf\n",n,x);

n++;

y=f(x);

y=fabs(x-y);

x=f(x);

}while(y>e);

}

­

/*本代碼只針對實驗二中的函數，實際中可將其中的定義的double f(double x)改為要求函數即可*/

3.牛頓法（還是對實驗二中的函數）

程序代碼：

#include <stdio.h>

#include <math.h>

double f(double x)

{

return (2*sin(x)-x);

}

double df(double x)

{

return (2*cos(x)-1);

}

main()

{

int n=0;

double x,y,e;

printf("請輸入X的初值 X0:\n");

scanf("%lf",&x);

printf("請輸入精度e的值 e:\n");

scanf("%lf",&e);

do{

printf("x%d=%lf\n",n,x);

n++;

y=(x-f(x)/df(x));

y=fabs(x-y);

x=(x-f(x)/df(x));

}while(y>e);

}

/*本代碼還只是針對實驗二中的函數，實際中可將其中的定義的double f(double x)、double df(double x)改為要求函數，double df(double x)為double f(double x)的導數。*/

4.雅可比法（針對實驗三）

程序代碼

#include <math.h>

#include <stdio.h>

void main(void)

{

int i,j,n,t,r,f,k=0,h=0;

double e;

float a[100][100],b[100],x[100],y[100];

printf("請輸入元數 n：\n");

scanf("%d",&n);

printf("請輸入系數矩陣 a[*][*]:\n");

for(i=0;i<n;i++)

for(j=0;j<n;j++)

scanf("%f",&a[i][j]);

printf("請輸入 b[*]:\n");

for(i=0;i<n;i++)

scanf("%f",&b[i]);

printf("請輸入精度要求 e:\n");

scanf("%lf",&e);

printf("請輸入最大允許迭代次數 t:\n");

scanf("%d",&t);

printf("請輸入x[*]的初值 x[0]:\n");

for(i=0;i<n;i++)

scanf("%f",&x[i]);

do

{

k++;

for(i=0;i<n;i++)

{

r=0;

for(f=0;f<n;f++)

{ if(f!=i)

r=r+a[i][f]*x[f];

}

y[i]=(1/a[i][i])*(b[i]-r);

if(k==t)

h=0;

else

{

if(fabs(x[i]-y[i])>e)

h=1;

else

h=0;

}

x[i]=y[i];

}

}while(h==1);

printf("答案是：\n");

for(i=0;i<n;i++)

printf("x%d=%f ",i,x[i]);

printf("\n");

printf("迭代次數是：%d\n",k);

}

7. 求直接三角分解法的matlab程序，每一步都需要解釋，謝謝！

你的代碼我幫你解釋了，果然是個復雜的活，如下function hl=zhjLU(A) %函數名為LU，輸入矩陣A,來進行LU分解，返回值hl為A的各階主子式的行列式[n n] =size(A);%返回矩陣A的維數
RA=rank(A); %返回矩陣A的秩if RA~=n %判斷矩陣A的秩RA是否不等於A的維數n，當不等於n時，即小於n時執行if中的語句disp('請注意：因為A的n階行列式hl等於零，所以A不能進行LU分解.A的秩RA如下:'), RA,hl=det(A);%輸出disp中的字元串，RA的值,hl的值(此時它的值為矩陣A的行列式)return %並且退出程序end % if end 語句塊，end標志if語句結束if RA==n %判別A的秩是否等於n，當等於n時，即滿秩時執行下面的語句for p=1:n%從1到n循環，主要是想獲得1至n階主子式的行列式h(p)=det(A(1:p, 1:p));%A(1:p, 1:p)為A的前p行前p列，即A的p階主子式，從而h(p)為A的p階主子式的值end %for語句結束標志hl=h(1:n);%把A的各階主子式的行列式值賦給hl存儲.for i=1:n %從1到n循環，主要是判斷各階主子式是否為0if h(1,i)==0 %判段第i階主子式的行列式是否為0，若為0輸出下面的語句disp('請注意：因為A的r階主子式等於零，所以A不能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：'), hl;RA%輸出disp中的字元串，各階主子式的行列式的值，以及矩陣A的秩return %當有一個為0時退出程序end %if語句結束標志end %for語句結束標志if h(1,i)~=0 %如果執行到這一句，說明上邊的for循環都執行了且沒有return出去，則此時i的值為n，判斷第n階行列式是否為0，即還是判斷A的行列式是否不為0，若不為0則輸出下面語句disp('請注意：因為A的各階主子式都不等於零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：')%輸出disp中的字元串for j=1:n%對1到n循環U(1,j)=A(1,j);%把A的第一行中的各個值賦給U中第一行的各個變數end%for循環結束標志for k=2:n%從第2列到第n列進行循環，即依次給LU的第k列賦值for i=2:n%從2到n循環 ,與j的循環是相互聯系的，即當i>j時對下三角矩陣L賦值，當j>=i時對上三角U賦值for j=2:n%從2到n循環 ,與i的循環是相互聯系的，即當i>j時對下三角矩陣L賦值，當j>=i時對上三角U賦值L(1,1)=1;L(i,i)=1; %始終保持L對角線上的元素為1if i>j %當i>j時，此時就是要對L進行賦值了，因為當i<j時,L(i,j)=0,不用管L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);%根據LU分解的關系知,L(i,1)處的值=A(i,1)/U(1,1)，從而對L(i,1)進行賦值，即對L的第一列賦值L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);%還是根據LU分解的關系L(i,k)處的值等於(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k)，對L(i,k)賦值，即對L的第k列賦值else %否則,即i<j時，此時就是要對U進行賦值了，因為當i>j時,U(i,j)=0,不用管U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);%根據LU分解的關系知U(k,j)的值等於A(k,j)-L(k,1:k-1)*U(1:k-1,j)進而對U的第k列賦值end %if結束標志end %for j=2:n結束標志end %for i=2:n結束標志end %for k=2:n結束標志hl;RA,U,L %輸出矩陣的秩RA，上三角函數U,下三角函數Lend %if h(1,i)~=0結束標志end %RA==n結束標志 %%%%%%%%%%%%%%%%%%%%%%%%以上代碼中有很多多餘的，當判斷A的秩為n之後其他的主子式的行列式都不為0了，判斷是多餘的，故我進行了改進,如下%%%%%%%%%%%%%%%%%%%%%%% function hl=LUfenJie(A) %函數名為LU，輸入矩陣A,來進行LU分解，返回值hl為A的各階主子式的行列式[n n] =size(A);%返回矩陣A的維數
RA=rank(A); %返回矩陣A的秩if RA~=n %判斷矩陣A的秩RA是否不等於A的維數n，當不等於n時，即小於n時執行if中的語句disp('請注意：因為A的n階行列式hl等於零，所以A不能進行LU分解.A的秩RA如下:'), RA,hl=det(A);%輸出disp中的字元串，RA的值,hl的值(此時它的值為矩陣A的行列式)return %並且退出程序elsefor p=1:n%從1到n循環，主要是想獲得1至n階主子式的行列式h(p)=det(A(1:p, 1:p));%A(1:p, 1:p)為A的前p行前p列，即A的p階主子式，從而h(p)為A的p階主子式的值end %for語句結束標志hl=h(1:n);%把A的各階主子式的行列式值賦給hl存儲.
disp('請注意：因為A的各階主子式都不等於零，所以A能進行LU分解.A的秩RA和各階順序主子式值hl依次如下：')%輸出disp中的字元串for j=1:n%對1到n循環U(1,j)=A(1,j);%把A的第一行中的各個值賦給U中第一行的各個變數end%for循環結束標志for k=2:n%從第2列到第n列進行循環，即依次給LU的第k列賦值for i=2:n%從2到n循環 ,與j的循環是相互聯系的，即當i>j時對下三角矩陣L賦值，當j>=i時對上三角U賦值for j=2:n%從2到n循環 ,與i的循環是相互聯系的，即當i>j時對下三角矩陣L賦值，當j>=i時對上三角U賦值L(1,1)=1;L(i,i)=1; %始終保持L對角線上的元素為1if i>j %當i>j時，此時就是要對L進行賦值了，因為當i<j時,L(i,j)=0,不用管L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1);%根據LU分解的關系知,L(i,1)處的值=A(i,1)/U(1,1)，從而對L(i,1)進行賦值，即對L的第一列賦值L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k);%還是根據LU分解的關系L(i,k)處的值等於(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k)，對L(i,k)賦值，即對L的第k列賦值else %否則,即i<j時，此時就是要對U進行賦值了，因為當i>j時,U(i,j)=0,不用管U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);%根據LU分解的關系知U(k,j)的值等於A(k,j)-L(k,1:k-1)*U(1:k-1,j)進而對U的第k列賦值end %if結束標志end %for j=2:n結束標志end %for i=2:n結束標志end %for k=2:n結束標志hl;RA,U,L %輸出矩陣的秩RA，上三角函數U,下三角函數Lend %RA~=n結束標志