19.5.3 有限元方法 (FEM)

在现代计算机出现后, 有限元方法成为求解偏微分方程的最重要的技术. 也容易说明这种强有力的方法产生的结果.

对于不同类型的应用, 有限元方法有很不相同的实施方式, 故这里仅介绍其基本思想. 它类似于数值求解常微分方程边值问题的里茨法 (参见第 1266 页 19.4.2.2), 也涉及样条插值 (参见第 1293 页 19.7).

有限元方法有如下步骤.

1. 定义变分问题

由已知边值问题形成变分问题. 以如下边值问题为例说明其过程:

在 $G$ 的内部 $\mathrm{\Delta }u=u_{xx}+u_{yy}=f$ ,在 $G$ 的边界上 $u=0$ .(19.144)

在微分方程 (19.144) 两边乘以在 $G$ 的边界为零的适当的光滑函数,并在整个区域 $G$ 上积分,得到

$$\begin{array}{}\text{(19.145)}& {\iint }_{\left(G\right)}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)v\mathrm{d}x\mathrm{d}y={\iint }_{\left(G\right)}fv\mathrm{d}x\mathrm{d}y.\end{array}$$

应用高斯求积公式 (参见第 945 页 13.2.2.1,(2)),将 $P\left(x,y\right)=-v{u}_y,Q\left(x,y\right)=v{u}_x$ 代入 (13.121), 从 (19.145) 得到变分方程

$$\begin{array}{}\text{(19.146a)}& a\left(u,v\right)=b\left(v\right)\end{array}$$

其中

$$\begin{array}{}\text{(19.146b)}& a\left(u,v\right)=-{\iint }_{\left(G\right)}\left(\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y}\right)\mathrm{d}x\mathrm{d}y,b\left(v\right)={\iint }_{\left(G\right)}fv\mathrm{d}x\mathrm{d}y.\end{array}$$

2. 三角形剖分

将积分区域 $G$ 分解为简单的子区域,一般用三角形剖分,即区域 $G$ 被三角形覆盖, 且其中相邻的三角形或共有整条边, 或仅有一个公共顶点. 每个带曲线边界的区域可由三角形的联合很好地逼近 (图 19.7).

注 为避免数值上的困难, 三角形剖分中应该不包含钝角三角形.

单位正方形的三角形剖分见图 19.8. 从坐标为 $x_{\mu }=\mu h,y_{\nu }=\nu h(\mu ,\nu =0,1$ , $2,\cdots ,N;h=1/N)$ 的网格剖分节点出发,有 ${(N-1)}^2$ 个内节点. 考虑到选取解函数,常用以 $\left(x_{\mu },y_{\nu }\right)$ 为公共顶点的六个三角形组合成的面元 $G_{\mu \nu }$ (在其他的情况下, 三角形可能不是 6 个. 这些面元显然互不排斥).

3. 求解

对所求的函数 $u\left(x,y\right)$ 在每一个三角形上定义假设的近似解. 带相应假设解的三角形称为有限元. $x$ 和 $y$ 的多项式是最合适的选择. 在许多情况下,线性近似

$$\begin{array}{}\text{(19.147)}& \tilde{u}\left(x,y\right)=a_1+a_{2}x+a_{3}y\end{array}$$

已经足够. 假设近似函数在相邻三角形间必须连续, 故最终的解也是连续的.

(19.147) 中的系数 $a_1,a_2,a_3$ 由在三角形顶点上的函数值 $u_1,u_2,u_3$ 唯一确定. 这同时保证了相邻三角形间的连续性. 假设解 (19.147) 包含了作为未知参数的要求函数的近似值 $u_i$ . 在整个区域 $G$ 上用来逼近要求的解 $u\left(x,y\right)$ 的假设解选为

$$\begin{array}{}\text{(19.148)}& \tilde{u}\left(x,y\right)=\sum _{\mu =1}^{N-1}\sum _{\nu =1}^{N-1}{\alpha }_{\mu \nu }{u}_{\mu \nu }\left(x,y\right).\end{array}$$

确定适当的系数 ${\alpha }_{\mu \nu }$ . 对函数 $u_{\mu \nu }\left(x,y\right)$ 必须成立: 根据 (19.147),它们在 $G_{\mu \nu }$ 的每一个三角形上是满足下面条件的线性函数:

(1)

$$\begin{array}{}\text{(19.149a)}& u_{\mu \nu }\left(x_k,y_l\right)=\{\begin{array}{ll}1,& k=\mu ,l=\nu ,\\ 0,& \text{ 其他. }\end{array}\end{array}$$

(2)

$$\begin{array}{}\text{(19.149b)}& u_{\mu \nu }\left(x,y\right)\equiv 0\text{,对任意}\left(x,y\right)\notin G_{\mu \nu }\text{.}\end{array}$$

$u_{\mu \nu }\left(x,y\right)$ 在 $G_{\mu \nu }$ 上的图形表示见图 19.9. 在 $G_{\mu \nu }$ 上,即在图 19.8 中从 1 到 6 所有的三角形上,计算 $u_{\mu \nu }\left(x,y\right)$ . 这里仅说明三角形 1 上的计算:

$$\begin{array}{}\text{(19.150)}& u_{\mu \nu }\left(x,y\right)=a_1+a_{2}x+a_3,\end{array}$$

满足

$$\begin{array}{}\text{(19.151)}& u_{\mu \nu }\left(x,y\right)=\{\begin{array}{ll}1,& x=x_{\mu },y=y_v,\\ 0,& x=x_{\mu -1},y=y_{v-1},\\ 0,& x=x_{\mu },y=y_{v-1}.\end{array}\end{array}$$

由 (19.151) 得 $a_1=1-v,a_2=0,a_3=1/h$ ,随之对三角形 1 有

$$\begin{array}{}\text{(19.152)}& u_{\mu \nu }\left(x,y\right)=1+\left(\frac{y}{h}-v\right).\end{array}$$

类似有

$$\begin{array}{}\text{(19.153)}& u_{\mu \nu }\left(x,y\right)=\{\begin{array}{ll}1-\left(\frac{x}{h}-\mu \right)+\left(\frac{y}{h}-v\right),& \text{ 三角形 }2,\\ 1-\left(\frac{x}{h}-\mu \right),& \text{ 三角形 }3,\\ 1-\left(\frac{y}{h}-v\right),& \text{ 三角形 }4,\\ 1+\left(\frac{x}{h}-\mu \right)+\left(\frac{y}{h}-v\right),& \text{ 三角形 }5,\\ 1+\left(\frac{x}{h}-\mu \right),& \text{ 三角形 }6.\end{array}\end{array}$$

4. 计算解的系数

解的系数 ${\alpha }_{\mu \nu }$ 由解 (19.148) 对每个解函数 $u_{\mu \nu }$ 都满足变分问题 (19.146a) 来确定,即在 (19.146a) 中用 $\tilde{u}\left(x,y\right)$ 代替 $u\left(x,y\right)$ ,用 $u_{\mu \nu }\left(x,y\right)$ 代替 $v\left(x,y\right)$ . 于是得到关于未知系数的线性方程组

$$\begin{array}{}\text{(19.154)}& \sum _{\mu =1}^{N-1}\sum _{\nu =1}^{N-1}{\alpha }_{\mu \nu }a\left(u_{\mu \nu },u_{kl}\right)=b\left(u_{kl}\right)\left(k,l=1,2,\cdots ,N-1\right),\end{array}$$

其中

$$a\left(u_{\mu \nu },u_{kl}\right)={\iint }_{G_{kl}}\left(\frac{\partial u_{\mu \nu }}{\partial x}\frac{\partial u_{kl}}{\partial x}+\frac{\partial u_{\mu \nu }}{\partial y}\frac{\partial u_{kl}}{\partial y}\right)\mathrm{d}x\mathrm{d}y,b\left(u_{kl}\right)={\iint }_{G_{kl}}f{u}_{kl}\mathrm{d}x\mathrm{d}y.$$

(19.155)

在计算 $a\left(u_{\mu \nu },u_{kl}\right)$ 时,必须记住仅需计算区域 $G_{\mu \nu }$ 和 $G_{kl}$ 具有非空交集情况下的积分. 这些区域在表 19.1 中用阴影表示.

因总是在面积为 $h^2/2$ 的三角形区域上积分,故关于变量 $x$ 的偏导数为

$$\begin{array}{}\text{(19.156a)}& \frac{1}{h^2}\left(4{\alpha }_{kl}-2{\alpha }_{k+1,l}-2{\alpha }_{k-1,l}\right)\frac{h^2}{2}.\end{array}$$

类似得到关于变量 $y$ 的偏导数

$$\begin{array}{}\text{(19.156b)}& \frac{1}{h^2}\left(4{\alpha }_{kl}-2{\alpha }_{k,l+1}-2{\alpha }_{k,l-1}\right)\frac{h^2}{2}.\end{array}$$

计算 (19.154) 的右端项 $b\left(u_{kl}\right)$ 得

$$\begin{array}{}\text{(19.157a)}& b\left(u_{kl}\right)={\iint }_{G_{kl}}f\left(x,y\right)u_{kl}\left(x,y\right)\mathrm{d}x\mathrm{d}y\approx f_{kl}{V}_P,\end{array}$$

其中 $V_P$ 为由 $u_{kl}\left(x,y\right)$ 确定的 $G_{kl}$ 上高度为 1 的棱锥的体积 (图 19.9). 因为

$$\begin{array}{}\text{(19.157b)}& V_P=\frac{1}{3}\cdot 6\cdot \frac{1}{2}{h}^2\text{,故其近似值为}b\left(u_{kl}\right)\approx f_{kl}{h}^2\text{.}\end{array}$$

于是变分方程 (19.154) 导致确定解系数的线性方程组

表 19.1 有限元法附表

面域	图示	三角形	$\frac{\partial u_{kl}}{\partial x}$	$\frac{\partial u_{\mu \nu }}{\partial x}$	$\sum \frac{\partial u_{kl}}{\partial x}\frac{\partial u_{\mu \nu }}{\partial x}$
		$G_{kl}{G}_{\mu \nu }$
(1) $\begin{array}{l}\mu =k\\ \nu =l\end{array}$		$\begin{array}{ll}1& 1\end{array}$	0	0	$\frac{4}{h^2}$
		22	$-1/h$	$-1/h$
		33	$-1/h$	$-1/h$
		$\begin{array}{ll}4& 4\end{array}$	0	0
		55	$1/h$	$1/h$
		6	$1/h$	$1/h$
					0
(2) $\begin{array}{l}\mu =k\\ \nu =l-1\end{array}$		$\begin{array}{ll}1& 5\\ 2& 4\end{array}$	0	$1/h$
			$-1/h$	0
(3) $\begin{array}{l}\mu =k+1\\ \nu =l\end{array}$		$\begin{array}{ll}2& 6\\ 3& 5\end{array}$	$-1/h$	$1/h$	$-\frac{2}{h^2}$
			$-1/h$	$1/h$
					0
(4) $\begin{array}{l}\mu =k+1\\ \nu =l+1\end{array}$			$-1/h$	0
		$\begin{array}{ll}3& 1\\ 4& 6\end{array}$	0	$1/h$
(5) $\begin{array}{l}\mu =k\\ \nu =l+1\end{array}$		$\begin{array}{ll}4& 2\\ 5& 1\end{array}$	0	$1/h$	0
			$-1/h$	0
(6) $\begin{array}{l}\mu =k-1\\ \nu =l\end{array}$		$\begin{array}{ll}5& 3\\ 6& 2\end{array}$	$1/h$	$-1/h$	$-\frac{2}{h^2}$
			$1/h$	$-1/h$
					0
			$1/h$	0
(7) $\begin{array}{l}\mu =k-1\\ \nu =l-1\end{array}$		$\begin{array}{ll}6& 4\\ 1& 3\end{array}$	0	$-1/h$

---

$$4{\alpha }_{kl}-{\alpha }_{k+1,l}-{\alpha }_{k-1,l}-{\alpha }_{k,l+1}-{\alpha }_{k,l-1}=h^2{f}_{kl}\left(k,l=1,2,\cdots ,N-1\right).\left(19.158\right)$$

---

(19.158)

注 (1) 若解的系数可由 (19.158) 确定,则由 (19.148) 得到的 $\tilde{u}\left(x,y\right)$ 表示其显式近似解,可对 $G$ 上任意点计算其值.

(2)若积分区域必须被不规则的三角形网格覆盖, 则要引入三角形坐标 (也称重心坐标). 这样容易确定点关于三角形网格的位置, 又因为每个三角形容易变换为以 $\left(0,0\right),\left(0,1\right),\left(1,0\right)$ 为顶点的单位三角形,从而使得 (19.155) 中多重积分的计算更容易.

(3) 若需要改善解的精度或可微性, 则必须应用分片二次或三次函数以得到假设的近似 (见 [19.28]).

(4) 在实际应用中, 通常得到高维的线性方程组. 这正是发展许多特殊方法的原因, 例如, 方程组的结构依赖于三角形网格的自动剖分和单元的实际列举. 关于有限元方法的详细介绍见 [19.21], [19.13], [19.28].