在有限元程式中，剛度方程[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（參見我前面寫的文章有限元非齊次邊界條件的剛度方程的化簡，用

用於儲存節點位移的陣列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一維陣列中搜索的方法。