在有限元程序中，刚度方程[K]建立完毕，节点力向量F经过了非齐次边界条件处理、等效节点力处理后，都搞成了已知量。此时，就可以解F=KD方程组，来求节点位移向量D了。
求解F=KD方程组的方法有很多，主要可以分为精确解法和迭代解法两种。顾名思义，精确解法就是直接解出D向量的精确值。迭代法则用于复杂的方程组，在精确解难以求出或求解成本太高时，通过一次次迭代近似求解，逼近精确解，当达到可接受的精度时，迭代停止。
对于商用有限元软件，由于模型太大太复杂，一般都使用迭代求解。
本篇介绍适用于小型方程组且能够得到精确解答的情况，使用LDLT方法解刚度方程F=KD，并用C++实现。
（1）LDLT方法的原理
LDLT分解法是一种基本的线性方程组直接解法。其优势是可以节约内存空间和减少计算量，其公式推导的资料在网上有一大堆，有兴趣的同学可以自行查找。此处仅简要说明。
对于对称正定矩阵A（例如，刚度矩阵[K]），可以将A分解为A =LDLT.并且该分解是唯一的。其中L是单位下三角矩阵，D为对角线矩阵，且对角线元素不为零。其中各矩阵的形式为：
A=
各矩阵元素之间的关系为：
为减少乘法次数，引入辅助量，则公式可写成
当设计程序时，为了减少存储空间，又由于矩阵A在经过上述计算后就不需要再使用了，因此可以采用原位存储的方法。也就是可以将dk存于矩阵A的对角线元素即akk中，uik存于矩阵A的上三角部分即aki中，而lik则可以存于矩阵A的下三角部分即aik中。
当需要求解对称线性方程组Ax=b（例如，刚度方程kd=f）时，由A=LDLT分解，有(LDLT)x=b，这可以分解成三个简单的方程组：Ly=b，Dz=y和LTx=z。由此易得求解公式为：
同样，向量b在经过计算后也不再需要，因此也可以采用相同的方法处理。可以在计算yi时将yi存于bi中，然后计算zi时也将zi存于bi中，最后计算xi时还是存于bi中，由bi输出计算结果即可。
（2）LDLT方法解刚度方程[F]=[K][D]的程序流程
1）前期准备：
采用一维变带宽存储的总刚矩阵pGK和辅助数组pDiag（参见我前面写的文章有限元刚度矩阵的一维变带宽存储用C++实现）；
经过非齐次边界条件和非节点力处理后的节点力数组pLoadVect（参见我前面写的文章有限元非齐次边界条件的刚度方程的化简，用C++实现和有限元模型中非节点外力的处理用C++实现）
用于存储节点位移的数组pDisp[n]，n为总刚矩阵[K]的阶数，也就是节点自由度总数。
2）程序流程图
（3）LDLT方法解刚度方程[K]的C++程序实现
bool LDLTDCMP (int n, double * *a) //将总刚矩阵K做LDLT分解，把L矩阵和u矩阵存储到完整总刚矩阵a[nTotalDOF][nTotalDOF]中
//n为外界传入的总刚矩阵阶数，即nTotalDOF；
//a是指向指针的指针（即，指向二维数组的指针），即为nTotalDOF行，nTotalDOF列的二维数组，存储完整的总刚矩阵，由外界传入
{
int k;
int m;
int i;
for (k = 0; k < n; k++)
{
/* Calculate d[k], and saved in a[k][k] */
for (m = 0; m < k; m++)
a[k][k] = a[k][k] - a[m][k] * a[k][m]; //从下一个for循环可知，a[m][k]=l[m][k],a[k][m]=u[m][k]
if (a[k][k] == 0)
{
cout<<"总刚系数有错！"<<endl;
return false;
}
for (i = k + 1; i < n; i++) //在此循环中i!=k，故a[k][k]=d[k]不会被覆盖！
{
/* temporary varible u[i][k] is saved in a[k][i]*/
for (m = 0; m < k; m++)
a[k][i] = a[k][i] - a[m][i] * a[k][m]; //a[k][i]=u[i][k]存储在矩阵a[n][n]的上三角区
/* Calculate l[i][k], and saved in a[i][k]*/
a[i][k] = a[k][i] / a[k][k]; //将a[i][k]=l[i][k]存储在a[n][n]矩阵的下三角区
}
}
return true;
}
void LDLTBKSB (int n, double * *a, double *b) //已知k=lu，计算kd=f中的d
//n为外界传入的总刚矩阵阶数，即nTotalDOF；
//a为LDLTDCMP（）函数计算修改了的a[n][n]矩阵，存储了L和U
//b为长度为nTotalDOF的一维数组，用于存储Ly=b，Dz=y，LTx=z中的y,z,x。最终，x即为kd=f中的位移向量d
{
int i;
int k;
for (i = 0; i < n; i++)
{
/* Calculate y[i], and saved in b[i] */
for (k = 0; k < i; k++)
b[i] = b[i] - a[i][k] * b[k]; //其中，a[i][k]=l[i][k]
}
for (i = n - 1; i >= 0; i--)
{
/* Calculate z[i], and saved in b[i] */
b[i] = b[i] / a[i][i];
/* Calculate x[i], and saved in b[i] */
for (k = i + 1; k < n; k++)
b[i] = b[i] - a[k][i] * b[k]; //其中，a[k][i]=u[i][k]
}
}
bool MyLDLTSolve(int nFreeDOF,double * *RpGK,double * pB) //自己写的解平衡方程求位移 // 用LDLT方法求解KD=F中，nFreeDOF未知节点位移部分的位移向量d的值，放入pD数组的对应位置中，即求解pDisp[0~nFreeDOF]是存在pD中
////nFreeFOF为总刚矩阵未知位移部分的阶数；RpGK为指向二维数组的指针，即指向总刚矩阵中未知位移部分组成的新矩阵(即，RpGK=pGK[0~nFreeDOF][0~nFreeDOF])，将要被修改为L和U矩阵；
//pB为存储最终计算结果，即kd=f中的d向量的一维数组，是外界传入的长度为nFreeDOF的一维数组，里面存入了总刚矩阵nFreeDOF部分的包含齐次化项的力向量，即pB=pLoadVect[0~nFreeDOF]
{
bool ind;
ind=LDLTDCMP (nFreeDOF, RpGK); //将总刚矩阵RpGK转存为了L和U矩阵
if(ind) //判断LDLT分解是否成功
{
LDLTBKSB (nFreeDOF, RpGK, pB);
return true;
}
else
{
cout<<"LDLT解方程出错"<<endl;
return false;
}
}
(4)缺点：这样普通的LDLT解方程，需要总刚矩阵[K]完全存储的临时二维数组RpGK，则一维变带宽的优势没有任何体现。下一篇中将会讨论LDLT解方程不再增加临时二维数组RpGK，直接在pGK一维数组中搜索的方法。