python 有限元分析_用python实现简单的有限元方法（一）

2022年 11月 5日

python 有限元分析_用python实现简单的有限元方法（一）

华中师范大学 hahakity

有限元算法(Finite Element Method，简称 FEM)是一种非常流行的求解偏微分方程的数值算法。有限元被广泛应用于结构受力分析、复杂边界的麦克斯韦方程求解以及热传导等问题。这一节介绍有限元方法的基本原理，以及如何用 Python 从头实现一个有限元算法，数值求解麦克斯韦方程。

学习内容筑基：加权残差法 (Weighted Residual Method)

心法：有限元与有限差分算法的区别

加权残差法(Weighted Residual Method )

理解了加权残差法，有限元算法就理解了一半。

这里参考文献【1】中例子，介绍加权残差法，为有限元算法筑基。

举例：求解

$equation?tex=x%28t%29$ , 使其满足如下微分方程，

$equation?tex=%7Bdx+%5Cover+dt+%7D+%2B+x+%3D+0$

初始条件

$equation?tex=x%280%29%3D1$ ,求解区域

$equation?tex=t+%5Cin+%5B0%2C+1%5D$ 。

这个方程的解析解是

$equation?tex=e%5E%7B-t%7D$ ,

$equation?tex=x%28t%29+%3D+e%5E%7B-t%7D+%3D+1+-+t+%2B+%7B1+%5Cover+2%21%7Dt%5E2++-+%7B1+%5Cover+3%21%7Dt%5E3++%2B+%7B1+%5Cover+4%21%7Dt%5E4+%2B+%5Ccdots$

在计算机上展开到无穷阶，就能给出方程的精确解。现实中显然只能考虑有限项，后面忽略的项称作截断误差。假设截断到二阶，

$equation?tex=x%28t%29+%5Capprox+1+-+t+%2B+%7B1+%5Cover+2%21%7Dt%5E2++$

代入原方程产生残差

$equation?tex=R$ ，

$equation?tex=%7Bdx+%5Cover+dt+%7D+%2B+x+%3D+R%2C+%5C%3B+R%5Cneq+0$

$equation?tex=R$ 越接近0，表示近似解越接近精确解。

明知道精确解

$equation?tex=e%5E%7B-t%7D$ 不在

$equation?tex=%281%2C+t%2C+t%5E2%29$ 这三个函数支撑起的函数空间中，我们仍可做出改进。

使用待定系数展开，

$equation?tex=x%28t%29+%5Capprox+1+%2B+c_1+t+%2B+c_2+t%5E2++$

此时残差函数

$equation?tex=R%28t%29+%3D+%7Bd+x%5Cover+dt%7D+%2B+x%28t%29+%3D+1+%2B+%281+%2Bt%29c_1+%2B+%282+t+%2B+t%5E2%29+c_2$

通过调节

$equation?tex=c_1%2Cc_2$ , 最小化

$equation?tex=R%28t%29$ 的绝对值，就能得到比

$equation?tex=1+-+t+%2B+%7B1+%5Cover+2%21%7Dt%5E2++$ 更好的解。

最优情况显然是

$equation?tex=R%28t%29+$ 对任意

$equation?tex=t+$ 都等于零，但这个做不到。所以存在不同的方法，减弱约束。下面介绍四种方法，点匹配法：在有限个点上，要求

$equation?tex=R%28t%29$ 等于0

子区域法：在有限子域上，要求

$equation?tex=R%28t%29$ 均值等于0

伽辽金法：对有限检验函数，要求

$equation?tex=R%28t%29$ 的加权求和等于0

最小二乘法：通过求

$equation?tex=%5Cint_0%5E1R%5E2%28t%29dt$ 的最小值，解待定系数

(1)点匹配法 (Point Matching)

最简单的考虑是，既然有

$equation?tex=c_1%2Cc_2$ 两个待定系数，那我们只需要任选两个点，让

$equation?tex=R%28t%29$ 在那两个点上等于 0，就可以得到两个方程，刚好可求解出

$equation?tex=c_1%2Cc_2$ ，

比如，选

$equation?tex=t+%3D1%2F3%EF%BC%8C2%2F3$ 两个点，得到两个方程。对第一个方程，

$equation?tex=R%5Cleft%28%7B1+%5Cover+3%7D+%5Cright%29+%3D+1+%2B+c_1%5Cleft%281+%2B+%7B1+%5Cover+3%7D%5Cright%29+%2B+c_2%5Cleft%282+%5Ctimes+%7B1%5Cover+3%7D%2B+%28%7B1%5Cover+3%7D%29%5E2+%5Cright%29+%3D+0$

化简得到，

$equation?tex=%7B4+%5Cover+3%7Dc_1+%2B+%7B7+%5Cover+9%7D+c_2+%3D+-1$

同理，对

$equation?tex=t+++%3D+2%2F3$ , 得到方程，

$equation?tex=%7B5+%5Cover+3%7Dc_1+%2B+%7B16+%5Cover+9%7D+c_2+%3D+-1$

解出

$equation?tex=c_1%3D-0.9310%2C%5C%3B+c_2%3D0.3103$ ，所以近似解为

$equation?tex=x%28t%29+%5Capprox+1+-+0.9310%5C+t+%2B+0.3103%5C+t%5E2++$

(2)子域法(Subdomain method)

根据待定系数的个数 n，将求解区域分成 n 份，要求

$equation?tex=R%28t%29$ 在每份上的积分等于0。对于这个例子，分

$equation?tex=t+%5Cin+%5B0%2C+1%2F2%5D$ 和

$equation?tex=t+%5Cin+%5B1%2F2%2C+1%5D$ 两个子域，求解如下两个方程，

$equation?tex=%5Cint_0%5E%7B1%2F2%7D+R%28t%29+dt+%3D+0%2C%5Cquad%5C%3B+%5Cint_%7B1%2F2%7D+%5E%7B1%7DR%28t%29+dt+%3D+0$

即可得到近似解，

$equation?tex=x%28t%29+%5Capprox+1+-+0.9474%5C+t+%2B+0.3158%5C+t%5E2++$

(3)伽辽金法(Galerkin method)

如果

$equation?tex=R%28t%29%3D0$ 对于任意

$equation?tex=t$ 都成立，显然

$equation?tex=%5Cphi%28t%29R%28t%29%3D0$ 应该对任意检验函数

$equation?tex=%5Cphi%28t%29$ 成立。

根据待定系数个数 n，选取 n 个不同的检验函数

$equation?tex=%5Cphi_i%28t%29%2C%5C+1+%5Cle+i+%5Cle+n$ , 要求对每个

$equation?tex=%5Cphi_i$ , 下式满足

$equation?tex=%5Cint_0%5E1%5Cphi_i%28t%29R%28t%29+dt+%3D+0$

这样，就构造了 n 个方程，刚好可解 n 个待定系数。

检验函数的一种自然选择是

$equation?tex=x%28t%29$ 展开时使用的基底函数

$equation?tex=%7B1%2C+t%5E1%2C+t%5E2%7D$ , 求解

$equation?tex=%5Cint_0%5E1+t+R%28t%29+dt+%3D+0%2C%5Cquad%5C%3B+%5Cint_0%5E1+t%5E2+R%28t%29+dt+%3D+0$

得到近似解，

$equation?tex=x%28t%29+%5Capprox+1+-+0.9143%5C+t+%2B+0.2857%5C+t%5E2++$

(4)最小二乘法 (Least Square Method)

函数极值处一阶导数等于0。求残差平方

$equation?tex=R%5E2%28t%29$ 积分的最小值，可以得到待定系数满足的方程，

$equation?tex=%5Cfrac%7B%5Cpartial%7D%7B%5Cpartial+c_%7B1%7D%7D+%5Cint_%7B0%7D%5E%7B1%7D+R%5E%7B2%7D%28t%29+%5Cmathrm%7Bd%7D+t%3D2+%5Cint_%7B0%7D%5E%7B1%7D+%5Cfrac%7B%5Cpartial+R%7D%7B%5Cpartial+c_%7B1%7D%7D+R%28t%29+%5Cmathrm%7Bd%7D+t%3D0$

对于

$equation?tex=c_1%2C+c_2$ 分别求导数，得到两个方程，

$equation?tex=%5Cint_%7B0%7D%5E%7B1%7D+%5Cfrac%7B%5Cpartial+R%7D%7B%5Cpartial+c_%7B1%7D%7D+R+%5Cmathrm%7B~d%7D+t%3D%5Cfrac%7B3%7D%7B2%7D%2B%5Cfrac%7B7%7D%7B3%7D+c_%7B1%7D%2B%5Cfrac%7B9%7D%7B4%7D+c_%7B2%7D%3D0+%5C%5C+%5Cint_%7B0%7D%5E%7B1%7D+%5Cfrac%7B%5Cpartial+R%7D%7B%5Cpartial+c_%7B2%7D%7D+R+%5Cmathrm%7B~d%7D+t%3D%5Cfrac%7B4%7D%7B3%7D%2B%5Cfrac%7B9%7D%7B4%7D+c_%7B1%7D%2B%5Cfrac%7B38%7D%7B15%7D+c_%7B2%7D%3D0+$

求解可得近似解，

$equation?tex=x%28t%29+%5Capprox+1+-+0.9427%5C+t+%2B+0.3110%5C+t%5E2++$

这四种方法可以用一个统一的式子描述，这个式子跟 Galerkin 给出的式子基本一致，

$equation?tex=%5Cint_%7B0%7D%5E%7B1%7D+w_%7Bj%7D%28t%29+R%28t%29+%5Cmathrm%7Bd%7D+t%3D0%2C+%5Cquad+j%3D1%2C2$

参考文献【1】总结了这四种方法对应的检验函数，如下所示，

有限元与有限差分算法的区别

一般的微分方程都可以写成

$equation?tex=L+f%28x%29+%3D+0$ 形式，其中

$equation?tex=L$ 称为微分算子(Differential Operator)。有限差分和有限元的区别在于有限差分对微分算子进行离散化近似，有限元对待求解函数进行基函数展开近似。

$equation?tex=%5Cunderbrace%7B%5Crm+Differential%5C+Operator%7D_%7B%5Crm+finite%5C+difference%7D+%5Cunderbrace%7Bf%28x%29%7D_%7B%5Crm+finite+%5C+element%7D%3D0$

有限差分算法(finite difference) 用有限差分来近似微分算子，比如一阶微分，

$equation?tex=%7Bdf+%5Cover+dx%7D+%5Capprox+%7Bf%28x%2B%5CDelta+x%29+-+f%28x+-+%5CDelta+x+%29+%5Cover+2%5CDelta+x%7D$

有限元算法(finite element) 用子域上的基函数展开来近似想要求解的方程

$equation?tex=f%28x%29$ ,

$equation?tex=f%28x%29+%3D+%5Csum_%7Bk%3D1%7D%5En+a_k+%5Cphi_k%28x%29$

有了前面加权残差法(Weighted Residual Method)的介绍，相信大家对上面这种基函数展开已经有了比较深刻的理解。有限元算法正是广泛使用了上面的 Subdomain 方法与 Galerkin 方法的思想。

上面例子我们知道解析解，所以选择了泰勒展开时前几个最重要的多项式函数作为基底。如果不知道解的性质，如何选择基底函数？有限元算法应用这么广泛的原因就在于它选择基底函数的方式简单粗暴有用，留在下节详细介绍。

总结

这节介绍了加权残差法，有限元与有限差分的区别。这种使用有限个基底函数展开，再调整待定系数，补偿截断误差的思想，还能在量子力学变分法，量子多体的Ab initio方法，重整化群，密度矩阵重整化群中找到身影，应该说是一种非常有用的计算物理思想。

参考文献：

【1】Advanced Electromagnetic Computation，Dikshitulu K. Kalluri