提示：本篇博客仅提供函数代码，不进行数学推导！需要了解数学推导的uu请移步其他博客！
笔者这学期选修了数值计算这门课程。一次实验报告中需要自行编写LDLT分解的函数，并在脚本文件中调用此函数，返回求得的x、L、D。
可能由于这个方法较为冷门？（瞎猜的，别cue我）网上相关的资料很少，能找到的、不用开vip就能看的文章只有个位数。。。因此能参考的内容实在太少。
求人不如求己！ 自！己！来！
所以笔者的灵感来自于老师已经给出的Cholesky分解、LU分解的函数文件上进行修改，最终完成了LDLT方法的求解函数。（事实上，老师布置这道题的初衷就是希望我们能在他给出的函数文件上进行修改，进而得到LDLT分解的函数）
秉着乐于助人的原则（always），遂把此代码分享在各个平台。愿在世界各地因同样问题正在挠头的你们，能看到这篇博客吧。
ps：笔者是刚入门matlab和数值计算的菜鸟，若以下代码有更加简略的写法，欢迎评论区留言或私信我~
提示：以下是本篇文章正文内容
1.假设我们要使用LDLT分解解如下一个方程组：
{ 4 x 1 − 2 x 2 − 4 x 3 = 10 − 2 x 1 + 17 x 2 + 10 x 3 = 3 − 4 x 1 + 10 x 2 + 9 x 3 = − 7 \begin{cases} \ 4x_1- 2x_2-4x_3=10\\ \ -2x_1+17x_2+10x_3=3\\ \ -4x_1+10x_2+9x_3=-7 \end{cases} ⎩ ⎨ ⎧ 4x1−2x2−4x3=10 −2x1+17x2+10x3=3 −4x1+10x2+9x3=−7
2.那么系数矩阵为A，A是如下的一个三阶对称矩阵：
[ 4 − 2 − 4 − 2 17 10 − 4 10 9 ] \left[ \begin{matrix} 4 & -2 & -4 \\ -2 & 17 & 10 \\ -4 & 10 & 9 \end{matrix} \right] 4−2−4−21710−4109经过初等变换后，A变为上三角矩阵DU。我们仅取DU的对角线元素，并赋值给全零阵D。
则DU如下：
[ 4 − 2 − 4 0 16 8 0 0 1 ] \left[ \begin{matrix} 4 & -2 & -4 \\ 0 & 16 & 8 \\ 0 & 0 & 1 \end{matrix} \right] 400−2160−481
所以根据上面的规则，矩阵D如下：
[ 4 0 0 0 16 0 0 0 1 ] \left[ \begin{matrix} 4 & 0 & 0 \\ 0 & 16 & 0 \\ 0 & 0 & 1 \end{matrix} \right] 4000160001
至于D为什么长这个样子，大家去看LDLT的数学推导就好，此处不再展开。
求得LDLT的D之后，接下来求L。这里有着繁琐的数学推导，此处不再展开*2。
废话少说，直接上LDLT分解的函数代码：
%LDL分解MATLAB程序function[x,L,D]=mldl(A,b)%用途：用平方根法解线性方程组Ax=b,A=LL'%格式：[x,L]=mchol(A,b),A为系数矩阵,b为右端向量%返回：x-解向量,L-下三角阵,D-对角阵n=length(b);
L=zeros(n,n);
L(logical(eye(size(L)))) = repmat(1,1); %L主对角线全是1DU = diag(basis_change(A)); %DU是把A经过初等行变换之后的上三角矩阵D = zeros(3);
D(logical(eye(size(D))))=DU; %把D对角线上的元素更改为DU的对角线上的元素for k=2:n
L(k:n,1) = A(k:n,1)/D(1,1);
L(k+1:n,k)=[A(k+1:n,k)-L(k+1:n,1:k-1)*D(k-1,k-1)*(L(k,1:k-1)')]/D(2,2);
end%解下三角方程组Ly=by=zeros(n,1);
y(1)=b(1);for k=2:n
y(k)=b(k)-L(k,1:k-1)*y(1:k-1);end%解上三角方程组DL'x=yx=zeros(n,1);
x(n)=y(n)/(D(n,n)*L(n,n)');
DLT = D*L';for k=n-1:-1:1
x(k)=(y(k)-DLT(k,k+1:n)*x(k+1:n))/DLT(k,k);end
最后求解x的时候，实际上还是偷懒用的LU分解的方法，因为我实在写不出x的推导式。。。
代码中的basis_change函数是我搜到的能够使用初等变换将矩阵变为上三角矩阵的代码，文章链接放在这里，把函数名修改成basis_change即可：
matlab中将任意矩阵转换成上三角矩阵的源码
最后在命令行或脚本文件中调用函数求解即可：
clc,clear
% 定义系数矩阵A和右端项b
A = [4 -2 -4;-2 17 10;-4 10 9]
b = [10 3 -7]';
[x,L,D] = mldl(A,b)
即可返回求得的x、L和D。
