三维有限元法计算过程

三维有限元法计算过程
三维有限元法的计算过程：
1）网格单元剖分；
2）线性插值；
3）单元分析；
4）总体刚度矩阵合成；
5）求解线性方程组等部分组成。

一、偏微分方程对应泛函的极值问题
矿井稳恒电流场分布示意图
主要任务是分析在给定边界条件下，求解稳定电流场的Laplace 方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题：
)
()((0)(0)()()(000z z y y x x I F u n u
n u F z u z y u y x u x Lu w D ---=⎪⎪⎪⎩
⎪⎪⎪⎨⎧=+∂∂=∂∂=∂∂∂∂+∂∂∂∂+∂∂∂∂≡ΓΓ+Γδδδγσσσ 上述微分方程边值问题等价于下面泛函的极小值问题：
dS U dxdydz fU z U y U x U U J w D ⎰⎰⎰⎰⎰Γ+Γ+ΓΩ
+-∂∂+∂∂+∂∂=222221
}])()()[(2{][γσσ
二、网格剖分
∞
1
ρi
i
h ρ...
...
1、网格单元的类型
图2-5 网格单元类型
2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则
1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而
不能成为其它单元内点；
2）、如果求解区域对称，那么单元剖分也应该对称；
3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓
的区域单元密度应小。

4）、网格单元体的大小变化应逐步过渡。

根据上述剖分原则，以x 、y 、z 坐标轴原点o 为中心，分别向x 、y 、z 方向的两侧作对称变步长剖分，距o 越远，步长应越大。

常用的变步长方法有：
c i x x i i )1(1+=∆-∆+ c x x i i =∆∆+/1(i ≠0)
c x x i i =∆-∆+1
1
1(i ≠0) 以上各式中c 为常数，1+∆i x 、i x ∆为同一坐标轴上相邻步长值。

以x 方向为例，可知，x 正方向与负方向对称，只相差一负号。

若令00=∆x ，只要给出距原点最近节点的坐标1x ∆，由上式即可求出其它相应的步长i x ∆。

同理可求得y 、
z 方向上的变步长i y ∆、i z ∆。

3、网格剖分方法
图2-6 平面内节点编号示意图
图2-7 长方体单元编号示意图
在程序设计的过程中，将采用四面体单元对研究区域进行分析。

因此对上面的长方体网格单元还需作进一步的处理。

任取一编号为m 的长方体单元，显然根据m 的值可求出其八个节点的节点编号。

为讨论方便，将长方体单元的各顶点的节点号，将其重新编号为0、1、2、3、4、5、6、7，如图2-8所示。

图2-8 长方体单元节点编号示意图
可知长方体单元可划分为六个四面体单元，其方法为用1、2、5、6顶点组成的平面将长方体分成二个三棱柱体，每个三棱柱体又可分成三个四面体单元，其结果如下：（0，1，2，4）、（1，2，4，5）、（2，4，5，6）、（1，2，3，5）、（2，3，5，6）、（3，5，6，7）。

其组合规律为：（j i +，2/)1(1+++j j i ，
2/)5(2j j i -++，j i ++4），1,0=i ，2,1,0=j 。

在具体的计算过程中，按长方体网格编号逐一计算，在每个长方体单元内再按以上的组合规律依次对六个四面体计算。

这样处理后四面体的总数为
1)(L *L_n *6Y -。

三、线性插值分析
任取一四面体单元如图2-9所示，设其顶点相应的节点序号为),,,(m l j i ，并设各节点电位为k u ，m l j i k ,,,=。

对应坐标为),,(i i i z y x ，),,(j j j z y x ，),,(l l l z y x ，
),,(m m m z y x ，当单元足够小时，设四面体内部电导率为常数、其电位呈
现线性变化。

即：
z y x U 4321αααα+++= （2-29）
图2-9 四面体单元示意图
对m l j i ,,,四个节点有：
⎪⎪⎩⎪⎪⎨
⎧+++=+++=+++=+++=m
m m
m l l l l j
j j j i i i
i z y x u z y x u z y x u z y x u 4321432143214321αααααααααααααααα （2-30）
整理式(2-29)、（2-30）得：
m m l l j j i i u N u N u N u N U +++= （2-31）
式中u 和k N 都是的),,(z y x 坐标函数，其中k N 的值为：
)(1
z d y c x b a D
N k k k k k +++=
，m l j i k ,,,= （2-32） 可知：k k k k d c b a ,,,（m l j i k ,,,=）、D 只与节点坐标有关。

为了分析式（2-32）所对应k N 的几何意义，现引入面积坐标和体
积坐标的概念。

如图2-10所示为ijl ∆，),(i i y x 、),(j j y x 、),(l l y x 为其顶点坐标，则三角形ijl ∆的面积为：
图2-10 三角形平面图
l
l j j i
i
ijl y x y x y x 1112
1
=
∆ 设),(y x p 点为三角形中任一点，则 l
j
l j l l
j j
i i i x x z
y y x y x y x
z c x b a 1111+-=++l
l
j j y x y x y
x
111= （2-33）
上式值为图2-10中三角形pjl ∆面积的两倍，所以：
ijl
pjl
ijl i i i i y c x b a y x N ∆∆=
∆++=
2)(),( （2-34） 同理可得：ijl
pil j y x N ∆∆=),(，ijl
pij l y x N ∆∆=
),(
定义),(y x N i ，),(y x N j ，),(y x N l 为),(y x p 在三角形ijl ∆内的面积坐标，具有以下性质：
1） 1),(),(),(=++y x N y x N y x N l j i （2-35） 2） 1),(≤y x N i 3） ⎩⎨
⎧≠==n
k n k y x N n n k 0
1
),(
4） 与节点i 对边平行的直线上各点的),(y x N i 值不变。

利用面积坐标的定义，可以方便地确定三角单元中任一点的面积坐标，且可建立一个面积坐标网，如图2-11所示。

因为面积坐标与通常使用的直角坐标系的选择无关，可利用该坐标网构造任意性线、非线性的插值函数。

图2-11 三角形面积坐标网
采用与上述相似的分析方法，可知),,(z y x N k (m l j i k ,,,=)为任一点
),,(z y x p 在四面体ijlm 内的体积坐标，具有与三角单元面积坐标相似的性质。

ijlm pjlm i D D z y x N =
),,( （2-36）
ijlm pilm j D D z y x N =
),,( （2-37）
ijlm pijm l D D z y x N =
),,( （2-38）
ijlm
pijl m D D z y x N =
),,( （2-39）
上式中ijlm D 为其对应节点所构成四面体的体积。

在式（2-31）假设四面体内电位为线性变化，P 点的电位值是通过该点的体积坐标将各节点的电位值线性化计算得出。

如果四面体内的电位值为非线性变化，则只要改变相应点体积坐标的系数即可。

形如：
m m l l j j i i u QN u PN u ON u SN U +++= （2-40）
式中S 、O 、P 、Q 为P 点体积坐标对应的系数，由于四面体内电位的非线性变化规律比较复杂，系数值的确定比较困难，因此在以下各章节的讨论中只考虑电位值的线性变化，对于非线性的变化规律将以后讨论。

四、网格单元分析
整个求解区域Ω被剖分为一系列小四面体单元后，泛函J （u ）也
相应地被离散为各单元上的泛函)(u J e ，故有∑=m
e u J u J )()(。

m 为四面体
单元总个数。

在任一四面体单元上泛函形式为：
⎰⎰⎰⎰⎰ΓΩ+-∂∂+∂∂+∂∂=e
e
dS u dxdydz fu z u y u x u
u J e 2222
21
}])()()[(2{)(σγσ
对任意节点的微分得：
⎰⎰⎰⎰⎰ΓΩ∂∂+∂∂-∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂=∂∂e
e
ds
u u
u dxdydz u u f z u u z u y u u y u x u u x u u u J k
k k k k k e σγσ})]()()([{)( （2-42）
其中m l j i k ,,,= 令 dxdydz z
u
u z u y u u y u x u u x u k k k e
)]()()([
∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂=Φ⎰⎰⎰Ωσ 则：
⎰⎰⎰⎰⎰ΓΩ∂∂+∂∂-Φ=∂∂e
e
ds u u
u dxdydz u u f u u J k k k e σγ)( （2-43）
现对各项分别进行讨论： 1、对于Φ项
由前面的讨论可知：∑=k
k k U N U ，)(1
z d y c x b a D
N k k k k k +++=
所以
m m l l j j i i
U x N U x N U x N U x N x U ∂∂+∂∂+∂∂+∂∂=∂∂ （2-44） x
N x U
U k k ∂∂=∂∂∂∂)( （2-45） 又因为：
D
b x N k k =∂∂，k k N U U
=∂∂ 所以：
)(1)(2m m k l l k j j k i i k k U b b U b b U b b U b b D
x U U x U +++=∂∂∂∂⋅∂∂ （2-46） 同理可得：
)(1)(2m m k l l k j j k i i k k U c c U c c U c c U c c D y U U y U +++=∂∂∂∂⋅∂∂ （2-47） )(1)(2m m k l l k j j k i i k k U d d U d d U d d U d d D
z U U z U +++=∂∂∂∂⋅∂∂ （2-48） 则：
)()()(z
U
U z U y U U y U x U U x U k k k ∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂+∂∂∂∂⋅∂∂ )]
()()[(1
2
m k m l k l j k j i k i m k m l k l j k j i k i m k m l k l j k j i k i U d d U d d U d d U d d U c c U c c U c c U c c U b b U b b U b b U b b D
+++++++++++=
]
)()()()[(1
2
m k m k m k m l k l k l k l j k j k j k j i k i k i k i U d d c c b b U d d c c b b U d d c c b b U d d c c b b D
+++++++++++=
（2-49） 记：
m
k m k m k m l k l k l k l j
k j k j k j i k i k i k i k U d d c c b b U d d c c b b U d d c c b b U d d c c b b F )()()()(+++++++++++=
所以：
k
e
k F V dv F D
e
362
σ
σ
=
=Φ⎰⎰⎰
Ω （2-50）
2、对于
dxdydz u u
f
e
k
⎰⎰⎰
Ω∂∂项
由前面分析可知：)()()(000z z y y x x I f ---=δδδ，
k k
N U U
=∂∂所以有： dxdydz z z y y x x z y x N I dxdydz u u
f
Ve
k Ve
k
)()()(),,(000---=∂∂⎰⎰⎰⎰⎰⎰
δδδ （2-51） 由δ函数的性质可知：
I z y x N dxdydz u u
f
k Ve
k
⋅=∂∂⎰⎰⎰
),,(000 ⎩⎨⎧=)
,,(0
),,(k k k k k k z y x z y x I 供电点不在点供电点在点 （2-52）
3、对于
⎰⎰
Γ∂∂e
ds u u
u
k
σγ项
由式（2-25）分析可知，该项只与研究区域Ω的第三类边界条件有关，所以当四面体完全位于求解区域内部时（即非边界单元），该项值为零。

在边界处则须对该项进行计算。

1） 关于γ值的计算
由电场理论知识可知，对于点源和均匀介质的情况，其电位函有如下的形式：
r
C
z y x U =
),,( （C 为常数） （2-53） 则：
θcos 3r
U
n r r C n U -=⋅=∂∂ （2-54） 式中θ为由点源到计算点的径向矢量r 和外法线矢量n 之间的夹角。

式（2-54）可改写成如下形式：
0cos =+∂∂θr
U
n U （2-55） 对比混合边界条件（2-7）的公式可知：
θγcos 1r
= （2-56）
设边界四面体单元e Ω如图2-9所示，由节点l j i ,,组成的面为边界面，另一节点m 位于研究区域内。

在三角形ijl ∆内，γ值是坐标),,(z y x 及θcos 的函数，为简化计算，当三角形的面积较小时，用其重心处的γ值来代替整三角形的γ值。

设三角形重心P 的坐标为),,(t t t z y x ，则：
⎪⎪⎪
⎩
⎪⎪
⎪⎨⎧++=++=++=)(31)(31
)(31t j i t t j i t t j i t z z z z y y y y x x x x （2-57）
可知：
202020)()()(z z y y x x r t t t -+-+-= （2-58）
由空间解析几何与向量代数理论可知，op
op
n n n n ⋅⋅=θcos
式中，op n
为由供电点O 到三角形重心的径向矢量，
),,(000z z y y x x n t t t op ---=
（2-59）
lj li n n n ⨯=
其中li n 、lj n
分别为节点i l ,和节点j l ,所在直线的径向矢量。

),,(l i l i l i li z z y y x x n ---=
（2-60） ),,(l j l j l j lj z z y y x x n ---=
（2-61）
可得：lj li n n
⨯
k
x x y y y y x x j z z x x x x z z i y y z z z z y y l j l i l j l i l j l i l j l i l j l i l j l i
)])(())([()])(()
)([()])(())([(-----+-----+-----= （2-62） 令 ))(())((l j l i l j l i y y z z z z y y A -----=
))(())((l j l i l j l i z z x x x x z z B -----=
))(())((l j l i l j l i x x y y y y x x C -----=
则： ),,(C B A n =
所以：
222202020000])()()[()()()(cos C
B A z z y y x x z z
C y y B x x A r t t t t t t ++⋅-+-+--+-+-=
=θ
γ 2） 关于
⎰⎰
Γ∂∂e
ds u u u
k
值的计算
由式（2-31）计算可知
k k
N U U
=∂∂，所以： dS U N N U N N U N N U N N dS U U
U
e
e
m k m l k l j k j i k i k
⎰⎰⎰⎰ΓΓ+++=∂∂)( （2-64） 为方便计算，设四面体单元e Ω的边界面ijl ∆所在的平面为XOY 平面，l 为原点，li 为X 轴，并设节点m l j i ,,,的坐标分别为)0,0,(i x ，)0,,(j j y x ，)0,0,0(和
),,(m m m z y x 。

再设三角形内的一动点P 的坐标为)0,,(y x 。

可知：m j i m
m
j m
j
i z y x z y y x x x D ==
0001
111，2j i y x e =∆ 由式（2-36）~（2-39）所示的体积坐标可知：
j i j i m m j m j i y x y x x x z y y y x x x
D N -
==
000001
1
11
1,j m
m m i
j y y
z y y x x x D N =
=0
000001111
1 i
j i i j m
m j m j i
l x x y y x x x z y y y x x x x D N -+-==
1000
01
111
1,00
00
00
001111
1==
y y x x x D N j j i m
因为1=+++m l j i N N N N 。

所以在Δ上1=++l j i N N N 。

故i N ，j N ，l N 中
只有两个自由量，设为i N ，j N ，则
dxdy N N j
i
⎰⎰∆
dx y y y x y x x x dy j
j i j y y
y x x x y x i j j
i j i j
j
)(
-=⎰⎰-+
dy y x x x y x y x y y x y y x x x y x y
i
i i j i j j j j i j i y j i j
)}(])()[(2{22220
----+=⎰
12
24
e
y x j i ∆=
=
（2-65） dx y y dy dxdy N N
j
j
i j i j
j
y y
y x x x y x j
j j
⎰⎰⎰⎰-+
∆
=0
22
6
12
e
y x j i ∆=
=
（2-66） 因为i N ，j N ，l N 中只有两个自由量，而且三角形Δ三个节点i,j,l 的选取排列是任意的。

因此：
l j i t t
k e
t k e dxdy N N t k ,,6
12
=⎪⎪⎩⎪⎪⎨⎧=∆≠∆=⎰⎰∆
（2-67）
所以：
ds u N N u N N u N N u N
N ds u u
u
m k m l k l j k j i k
i
k
)(+++=∂∂⎰⎰Γ
Γ
σγσγ
dxdy u N N u N N u N N u N N m k m l k l j k j i k i )(+++=⎰∆
σγ
m
k l k j k i k u u u u e e e e e e e e e
m
l j i ====⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎣
⎡⎥⎥⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎣
⎡∆∆∆∆∆∆∆∆∆=0000061212012612012126σγσγσγσγσγσγσγσγσγ （2-68）
记上式为k ψ
五、系数矩阵合成
由单元分析可知：
k k k k e N I F Ve
u u J ψ+⋅-=∂∂36)(σ
（2-69） 其中k ψ项只对参于对边界单元的计算。

由泛函极值的必要条件
0)
(=∂∂k
e u u J 可知： 036=ψ+⋅-k k k N I F Ve
σ
即：
k k k N I F Ve
⋅=ψ+36σ
（2-70）
式中 ⎩⎨
⎧≠==)
,,(),,(0
),,(),,(1
000000z y x z y x z y x z y x N k k k k k k k 供电点不在节点供电点在节点
对m l j i k ,,,=，上式可写成如下的矩阵形式[4]：
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢
⎢⎢⎢⎣
⎡m l j i m l j i e mm e lm e ll e jm e
jl e jj e
im
e
il e ij e ii IN IN IN IN u u u u k k k 对k k k k k k k 称 （2-71） 系数矩阵即为四面体单元Ve 的单元刚度矩阵e K ，显然单元刚度矩阵为对称阵，其阶数等于单元节点数。

在具体求解计算过程中，应将单元系数矩阵合成总的系数矩阵K 。

即
∑=m
e K K ,其中 C K C K e T e =为单元系数矩阵在总系数矩阵中的形式；C
为n ⨯4阶的转换矩阵，其作用将单元体系数矩阵变成总系数矩阵，其元素值为11=i C ，12=j C ，13=l C ，14=m C 其余值全为零。

整个求解区域离散为四面体单元，每一个单元体都可求得其对应的系数矩阵，将各单元刚度矩阵叠加即可得总系数矩阵。

设整个求解区域有n 个节
点，则K 是一个n n ⨯阶的系数矩阵。

合成总矩阵后，可得由n 个方程组成的大型线性方程组，其矩阵表达形式为：
⎥⎥⎥
⎥⎥
⎥
⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------1n i 101n i 101n 1,n j 1,n 1,1n 1,0
n 1n i,ij i1i01n 1,1j 11101n 0,0j 0100IN IN IN IN u u u u k k k k k k k k k k k k k k k k
（2-72） 可记作：
F
KU = （2-73）
系数矩阵K 的特点
总系数矩阵K 的阶数为n n ⨯，其数据量较大，在计算时对计算机的内存要求较高、计算时间较长。

为解决这一问题，可以利用总系数矩阵的某些特性，得到一定程度的解决。

有限元总系数矩阵具有以下特点：
1)对称性 总系数矩阵是一个对称矩阵,运行时系数矩阵可以只存贮矩阵的上三角部分或下三角部分,可节省一半的存贮量。

2）稀疏性 由于有限元方法采用四面体单元结点的相互关联性，这决定了总系数矩阵具有大型、稀疏的特征。

实际上，只有当结点i 是结点j 的相邻点时0k ≠ij 即总系数矩中的某一行中只有几个非零元素，对于一个四面体单元的四个结点，与任一结点对应的整体刚度矩阵中最多只有15个非零元素。

当空间结点数较多，网格剖分较密，总系数矩阵很大时其矩阵的稀疏性表现得越突出。

另一方面，总系数矩阵的结构（即非零元素的排列方式）依赖于结点的编号规律。

相邻结点的编号相差较大，其系数矩阵的稀疏性越明显。

3)带状分布特性 对于有限元数值分析中的总系数矩阵,非零元素分布在
以主对角线为中心的斜带形区域中,由于对称性,存贮系数时,可只取上三角(或下三角)部分,因而每行元素只存贮半个带宽的元素就可以了。

4)正定性有限元数值分析中,方程离散后形成的系数矩阵具有正定性特点,即主对角线元素在所处行与列为最大且为正值。

六、求解线性方程组
有限元法是求解复杂电磁场问题的有效工具之一，大型工程电磁场有限元计算最后都归结为大型稀疏有限元方程的求解。

在求解区域上对泊松方程进行离散并设置各节点参数及其边界条件之后，便得到一个大型的线性方程组。

求解此方程组，就可以得到各节点的电位值，进而计算出不同装置形式下的视电阻率值。

目前求解线性方程组F
KU 的方法很多：
1）、Gauss-Sidel迭代法；
2）、超松驰迭代法；
3）、高斯消去法（Gauss）；
4）、乔列斯基分解法(LDL T)等。

这些方法在求解线性方程组方面已经取得了较为广泛的应用，并且具有计算精度较高、稳定性强等优点。

本文在程序设计的过程中，主要采用了常规的Gauss消去法和适合大型稀疏矩阵方程的LDL T分解法，并对其计算精度进行对比。

1、Gauss消去法
Gauss消去法是计算机上常用的求解线性方程组的有效算法，此方法分为消元过程和回代过程。

消元过程是将原方程组化为上三角形方程组的过程，而
回代过程是求解上三角形方程的过程
对于线性方程组：
F KU = （3-1）
其中：n n )k (K ⨯=ij 为对称正定型稀疏矩阵
T n 21)u ,......,u ,(u U = T n 21)f ,......,f ,(f F =
为方便起见，计（3-1）式为(1)(1)F U K =，即：
⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡(1)n (1)2(1)1n 21(1)nn (1)n2
(1)n1
(1)2n (1)22(1)21(1)
1n
(1)12(1)11f f f u u u k k k k k k k k k （3-2）
第一次消元：
设0k (1)
11≠，记(1)11(1)11/k k l i i =，
（n ,,3,2 =i ），以1l i -乘式（3-2）中的第一个方程加到第i （n ,,3,2 =i ）个方程上，得(2)(2)F U K =，即：
⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡(2)n (2)2(1)1n 21(2)nn (2)n2(2)2n (2)22(1)
1n
(1)12(1)11f f f u u u k k k k k k k （3-3）
其中：)
1(11)1()2(k l k k j i ij ij -= （n ,,3,2, =j i ）
)1(11)1()2(f l f f i i i -= （n ,,3,2 =i ） 第二次消元：
设0k )
2(22≠，记)2(22)2(22k /k l i i =，（n ,,4,3 =i ），以2l i -乘式（3-3）中的第二个方程加
到第i （n i ,,4,3 =）个方程上，得(3)(3)F U K =，即：
⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣
⎡(3)n (3)3(2)2(1)1n 21(3)nn (3)n3(3)3n
(3)33(2)2n (2)23(2)22
(1)
1n
(1)13(1)12(1)11f f f f u u u k k k k k k k k k k k （3-4） 第1-k 次消元完成后，得与式（3-1）等价的方程组，)()(F U K k k =，即：
⎥⎥⎥
⎥⎥⎥
⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥
⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡)(n )((2)2(1)1n 21)
(nn )(n )(n (2)2n (2)22
(1)
1n
(1)
12(1)11f f f f u u u u k k k k k k k k k k k k k k k k
k k (k)kk
（3-5） 其中：)
1(11)1()(k l k k -----=k j k ik k ij k ij （n ,,1,, +=k k j i ） )1(11)1()(f l f f -----=k k ik k i k i （n ,,1, +=k k i ）
如此进行下去，共经过1-n 次消元计算后得到与式（3-1）等价的上三角形方程组，即：
⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣
⎡(n)n (2)2(1)1n 21(n)nn (2)2n (2)23(2)22
(1)
1n
(1)13(1)12(1)11f f f u u u k k k k k k k k （3-6） 然后利用回代公式求解线性方程组的解为：
⎪⎩
⎪⎨⎧-==∑+=)(n 1)()()
()(k /)u k f (u k /f u i ii
i s s i is i i i n nn n n n （1,2,,2n ,1n --=i ） （3-7） 2、LDL T 分解法
对有限元方程组F KU =的系数矩阵K 采用半带宽压缩存储技术，用二维数
组)M KN(N,存放系数矩阵K ，其中N 为网格节点总数，M 为半带宽。

将系数矩阵K 矩阵分解为：
T LDL K = （3-8） 式中，L 为带状的下三角矩阵，0l =ij )M (>-j i ；D 为对角阵ii
ii l 1d =
则： ∑
-=-=1
l l l k l j t t tt
jt it ij ij m
)t ;N 3,2,1(i j i m == （3-9）
式中：⎩⎨
⎧+>-+≤=1
M M
1M 0
t i i i m （3-10）
于是，解有限元方程F KU =等价于解方程F U LDL T =，U L Y T =，
则得方程组： ⎩⎨⎧==Y
U L F
LDY T （3-11）
)l f l (
f y 1
∑-=-=i t t tt
t
it i i m )N ,...,2,1(-i （3-12） 式中m t 值与式（3-10）相同，最后，矩阵U 各元素的计算公式为：
ii i t t ti
i i m
l /)u l
y (u t 1
∑+=-
= 1,,1N N, -=i （3-13）
⎩⎨
⎧+-≤+++->-=1
M N 1
M 1M N 1N t i i i m （3-14） 为比较以上两种方程组求解方法的计算精度及其稳定性，令式（2-73）中的U 向量值为1，求解右端项F 。

在程序计算过程中，利用总系数矩阵K 和所得的右端项F ，计算向量U 的值，比较U 与1的大小关系，即可判断比较方程组求解方法的计算精度及其稳定性。

具体计算结果如表3-2所示。

表3-2 计算结果比较
经过对上述两种方程组求解算法的计算结果比较可知，在相同条件下LDL T 法计算精度较高，稳定性好，但速度较慢；Gauss 消去法计算速度快，但其计算精度较低，稳定性差。

因此在计算过程中将采用LDL T 分解算法。

3、LU 分解法
高斯消元的过程相当于下述矩阵乘法运算
[][](2)(1)||M M A b U y =
因此由分块矩阵乘法可得
(2)(1)M M A U =(2)(1)M M b y = （3-2-15）
令 (2)(1)1(1)1(2)1()()()L M M M M ---==
可见只要消元过程能进行到底，就有以下等价关系
⎩
⎨⎧==−−→←=−−→←===y Ux b
Ly b LUx b Ax Ux
y LU
A （3-2-16） 消元过程相当于分解A 为单位下三角阵L 与上三角阵U 的乘积，解方程组Ly=b ，回代过程就是解方程组Ux=y 。

以上分析与结论完全可以推广到n 阶线性方程组，相应的式子（3-2-16）中的L 为n 阶单位下三角阵,U 为上三角阵
21
12
1
11n n l L l l ⎡⎤⎢⎥⎢⎥=⎢⎥⎢
⎥⎣⎦ 11121(1)
(1)222(1)n n n nn a a a a a U a -⎡⎤
⎢⎥⎢
⎥=⎢⎥⎢⎥⎣
⎦
解线性方程组Ax=b 可利用矩阵A 的LU 分解转化为解两个三角形方程，这种解方程组的方法称为LU 分解法或三角状分解法。

关于矩阵A 的LU 分解的存在唯一条件有如下两个定理。

由前面的分析可知，可以用顺序消去法求矩阵A 的LU 分解，为了简化分解，下面通过比较元素法导出一个直接求L 与U 的元素的公式。

令
A=LU （3-2-17） 其中
２６
213132
12
3
1111n n n l
L l l l l l ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦ 1112
13
12223233
3n n n nn u u u u u u u U u u u ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢
⎥⎢⎥⎣⎦
111111
(1,2,3
,)(2,3,4
,)
j j i i U a j n a l j n u ==⎧⎪
⎨==⎪⎩
1
1
1
1
(,1,,)
(1,2,
,)
k jk ks sj kj
s k ik is sk ik kk
s a l u u j k k n a l u l u j k k n -=-==+=+=+=++∑∑ （3-2-18）
1
1
1
1(,1,,)
(1,2,
,)
k kj kj ks sj
s k ik is sk s ik kk u a l u j k k n a l u l j k k n u -=-=⎧
=-=+⎪⎪⎪⎡⎤⎨
-⎢⎥⎪
⎣
⎦⎪==++⎪⎩
∑∑ （3-2-19）
（书P135：电测深曲线三维正演数值模拟）。