改进平方根法：解决地下水数值模拟基础问题

如前所述,用有限元法求解地下水流问题一般归结为求解具有对称、正定系数矩阵的线性代数方程组。下面我们来讨论对它的求解。

设求解的线性代数式为

其中系数矩阵[A]=(a i,j)设为正定阵。而

分别为有n个未知量组成的未知矢量和n个已知量组成的右端矢量。

当[A]的行列式不为零时,式(3-36)一定有唯一的解(薛禹群和谢春红,1980)。解的方法通常有两类:一类是迭代法,如超松弛迭代法等,将在第四章中介绍。这种方法往往费时较多。在结点不是太多,计算机存贮量许可的条件下,最好采用另一类直接法求解。

直接法中最常用的是改进平方根法(LDL T)或Cholesky法。这种方法充分利用矩阵对称的特点,压缩存贮量,求解手续简便易行,求解速度比较快,结果良好。

这种方法先把对角线元素均大于零的对称矩阵[A]分解为三个矩阵的乘积。这三个矩阵一个是下三角形矩阵,一个是上三角形矩阵(为下三角形矩阵的转置矩阵),另一个为对角线矩阵。然后将一个线性代数式(3-36)的求解,化为上三角形代数方程组和下三角形代数方程组两步来求解。即若

其中

[L]T为[L]的转置矩阵

则式(3-36)变为

如果令

则有

因而解式(3-36)[A]{H}={B}的问题变为先由式(3-40)[L]{y}={B}计算出中间变量{y},然后由式(3-39)[D][L]T{H}={y}求解出方程组的解{H}的问题。

由式(3-40)展开得(www.xing528.com)

因此,可自上而下地逐个确定y 1、y 2、…、y n。即由第一式得

把它代入第二式,解得

然后把解得的y 1、y 2代入第三式解出y 3

依次类推,最后可以解得y n。这些公式可以合并为一个递推公式

解得{y}后,下一步为解出方程组的解{H}。由于

所以[D][L]T{H}={y}可以写成下列形式

由此可以自下而上地确定H n、H n-1、…、H 1,即

这就是代数式(3-36)的解。

解找到了,但问题还没有全部解决。剩下的一个问题是如何由矩阵[A]推出组成矩阵[L]的各个元素l i,j(i≥j)。由于

所以

显然上式右端应等于[A]。由此得

利用上述关系式,可以依次求出各l i,j值。归结起来,可以得到下列求l i,j的递推公式为

解得l i,j后,就可以把矩阵[A]分解为[L][D][L]T,然后按式(3-41)、式(3-42)来确定{y},接着再按式(3-43)、式(3-44)来确定{H},从而求得线性代数式(3-36)的解。