直接三角分解的算法

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结束标志