有限元方法软件scikit-fem的两个例子

2023-07-08 10:07 作者:quantumtopology 0人读过 | 我要投稿

最近在使用有限元方法。有限元软件有很多，最出名的应该是fenics，它功能非常多，功能包装得也很好，使用方便。唯一的缺点是太大了，如果只是写个小软件去安装它就显得非常臃肿。想象一下，你使用外部库比你自身的程序大的多，而这个外部库还是不常用的。fenics在wsl上还不太好安装。因为windows下认为两个名字相同但大小写不同的文件是同一文件，但安装fenics时需要安装一个和已有库同名的软件包，所以就安装不上。后来找来找去找到一个名为scikit-fem的包。它很小，看它官网上帮助文档发现也能解很多问题。但缺点是这帮助文档写的太简略了！有好多代码和函数需要猜它什么意思。这了记录两个典型例子，可能会对使用scikit-fem的人有所帮助。

一维波动方程

该例子是ex44.py中的例子。对于方程 $u_%7Btt%7D%20%3D%20c%5E2u_%7Bxx%7D$ ，它需要转化成一个一阶方程组求解。所以这个例子对于求解方程组也很有帮助。

$%5Cbegin%7Bcases%7D%0A%5Cfrac%7B%5Cpartial%20u_1%7D%7B%5Cpartial%20t%7D%20%3D%20-u_2%20%5C%5C%0A%5Cfrac%7B%5Cpartial%20u_2%7D%7B%5Cpartial%20t%7D%20%3D%20-c%5E2%20%5Cfrac%7B%5Cpartial%5E2%20u_1%7D%7B%5Cpartial%20x%5E2%7D%0A%5Cend%7Bcases%7D$

使用Crank-Nicolson差分格式

$%5Cbegin%7Bcases%7D%0Au_1%5E%7Bk%2B1%7D%20-%20u_1%5Ek%20%26%3D%20-0.5%5CDelta%20t(u_2%5E%7Bk%2B1%7D%2Bu_2%5Ek)%20%5C%5C%0Au_2%5E%7Bk%2B1%7D%20-%20u_2%5Ek%20%26%3D%20-0.5c%5E2%5CDelta%20t%20%5Cleft(%20%5Cfrac%7B%5Cpartial%5E2%20u_1%5E%7Bk%2B1%7D%7D%7B%5Cpartial%20x%5E2%7D%20%20%20%2B%20%5Cfrac%7B%5Cpartial%5E2%20u_1%5E%7Bk%7D%7D%7B%5Cpartial%20x%5E2%7D%20%20%20%5Cright)%0A%5Cend%7Bcases%7D$

把其中第二个方程转化成变分形式（第一个方程不用转化，因为没有导数项），并整理得到

$%5Cbegin%7Bcases%7D%0Au_1%5E%7Bk%2B1%7D%2B0.5%5CDelta%20t%20u_2%5E%7Bk%2B1%7D%20%26%3D%20u_1%5E%7Bk%7D%20-%200.5%5CDelta%20t%20u_2%5Ek%20%5C%5C%0A%5Cint%20u_2%5E%7Bk%2B1%7D%20v%20-%200.5c%5E2%5CDelta%20t%5Cint%20%5Cfrac%7B%5Cpartial%20u_1%5E%7Bk%2B1%7D%7D%7B%5Cpartial%20x%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%20%26%3D%20%5Cint%20u_2%5Ek%20v%20%2B%200.5c%5E2%5CDelta%20t%20%5Cint%20%5Cfrac%7B%5Cpartial%20u_1%5E%7Bk%7D%7D%7B%5Cpartial%20x%7D%5Cfrac%7B%5Cpartial%20v%7D%7B%5Cpartial%20x%7D%0A%5Cend%7Bcases%7D$

现在看如下代码块

可以看到矩阵A和B的形式分别为

$%5Cbegin%7Bpmatrix%7D%0A%26I%20%260.5%5Cmathrm%7Bd%7DtI%20%5C%5C%0A%26-c%5E2L%200.5%5Cmathrm%7Bd%7Dt%20%26M%0A%5Cend%7Bpmatrix%7D%2C%20%5Cquad%0A%5Cbegin%7Bpmatrix%7D%0A%26I%20%26-0.5%5Cmathrm%7Bd%7DtI%20%5C%5C%0A%26c%5E2L%200.5%5Cmathrm%7Bd%7Dt%20%26M%0A%5Cend%7Bpmatrix%7D$

分别对应上面推导公式的右端和左端。时间t^k时刻的u乘以B，得到Bu^k，然后求解Au^{k+1} = Bu^k，得到u^{k+1}。

2. 非线性方程

非线性方程转化成变分形式会得到非线性泛函。此时把用基组展开的试探函数带入到变分形式中将得到关于展开系数的非线性方程。此时用牛顿迭代法去求解即可。在使用scikit-fem时，我们需要把变分形式线性化交给程序求解。一个泛函 $F%5Bu%5D$ 泛函导数的定义为

$%5Clim_%7Bt%5Crightarrow%200%7D%20%5Cfrac%7BF%5Bu%2Bt%5Cdelta%20u%5D%20-%20F%5Bu%5D%7D%7Bt%7D%3D%5Cfrac%7B%5Cdelta%20F%7D%7B%5Cdelta%20u%7D%5Cdelta%20u$

引入参数t是为了能和普通函数求导一样求出泛函导数，其实是不必要的。这里选用的例子是ex10.py。对于极小曲面问题（帮助文档里写错了，掉了根号）

$%5Cnabla%5Ccdot%5Cleft(%20%5Cfrac%7B%5Cnabla%20u%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%20%5Cnabla%20u%20%7C%7C%5E2%7D%7D%5Cright)%20%3D%200$

它的变分形式为

$F(u%2Cv)%20%3D%20%5Cint%20%5Cfrac%7B%5Cnabla%20u%5Ccdot%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%5Cmathrm%7Bd%7D%5COmega%20%20%3D%200$

为了求解泛函导数，把如下泛函在t = 0处展开到一阶（为了方便省略积分号）

$%5Cfrac%7B%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%7D%7D%20%3D%20%5Cfrac%7B%5Cnabla%20u%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%2B%20J%5Ccdot%20t$

$J%3D%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%7D%7D%20%0A-%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%20%5Cnabla%20v%7D%7B(1%20%2B%20%7C%7C%5Cnabla%20(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2)%5E%7B3%2F2%7D%7D%20%5Cfrac%7B1%7D%7B2%7D(2%5Cnabla%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u%2B2t%20%5Cnabla%20%5Cdelta%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u)$

注意这里使用了

$%7C%7C%5Cnabla(u%20%2B%20t%5Cdelta%20u)%7C%7C%5E2%20%3D%20%5Cnabla(u%20%2B%20t%5Cdelta%20u)%5Ccdot%20%5Cnabla(u%20%2B%20t%5Cdelta%20u)%5C$

当t趋于0时

$J(t%5Crightarrow%200)%3D%20%5Cfrac%7B%5Cnabla%20(u%20%2B%20%5Cdelta%20u)%20%5Cnabla%20v%7D%7B%5Csqrt%7B1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2%7D%7D%20%0A-%20%5Cfrac%7B%5Cnabla%20u%20%5Cnabla%20v%7D%7B(1%20%2B%20%7C%7C%5Cnabla%20u%7C%7C%5E2)%5E%7B3%2F2%7D%7D%20%5Cfrac%7B1%7D%7B2%7D(2%5Cnabla%20u%5Ccdot%20%5Cnabla%20%5Cdelta%20u)$