python求解偏微分方程组 python求解多元偏微分方程组

$-\nabla^{2} u(\boldsymbol{x})=f(\boldsymbol{x}), \quad \boldsymbol{x} \text { in } \Omega,\qquad(1)$

$u(\boldsymbol{x})=u_{\mathrm{D}}(\boldsymbol{x}), \quad \boldsymbol{x} \text { on } \partial \Omega,\qquad(2)$

其中，$u=u(\boldsymbol{x})$ 为未知函数，$f=f(\boldsymbol{x})$ 为规定函数，$\nabla^{2}$ 为拉普拉斯算子（常写为 $\nabla$），$\Omega$ 为空间域，$\partial \Omega$ 为 $\Omega$的边界。 泊松问题，包括偏微分方程 $-\nabla^{2} u=f$ 和 $\partial \Omega$ 上的边界条件 $u=u_{\mathrm{D}}$，是边界值问题的一个例子，在开始使用 FEniCS 解决它之前必须精确说明。

在坐标为 x 和 y 的二维空间中，我们可以写出泊松方程为

$-\frac{\partial^{2} u}{\partial x^{2}}-\frac{\partial^{2} u}{\partial y^{2}}=f(x, y),\qquad(3)$

未知数 $\boldsymbol{u}$ 现在是两个变量的函数，$u=u(x, y)$，在二维域 $\Omega$

泊松方程出现在许多物理环境中，包括热传导、静电、物质扩散、弹性杆的扭曲、无粘性流体流动和水波。 此外，该方程出现在更复杂的偏微分方程系统的数值分裂策略中，尤其是 Navier-Stokes 方程。

求解边界值问题包括以下步骤：

有限元变分公式

FEniCS 基于有限元方法，它是 PDE 数值解的通用且高效的数学机器。 有限元方法的起点是以变分形式表示的偏微分方程。

将 PDE 转化为变分问题的基本方法是将 PDE 乘以函数 $v$，在域 $\Omega$ 上对所得方程进行积分，并通过具有二阶导数的部分项进行积分。 乘以 PDE 的函数 $v$ 称为测试函数。 要逼近的未知函数 $u$

在本例中，我们首先将泊松方程乘以测试函数 $v$ 并在 $\Omega$

$-\int_{\Omega}\left(\nabla^{2} u\right) v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x,\qquad(4)$

我们在这里让 dx 表示域 Ω 上积分的微分元素。 我们稍后将让 ds 表示在 Ω 边界上积分的微分元素。

当我们推导出变分公式时，一个常见的规则是我们尽量保持 $u$ 和 $v$ 的导数的阶数尽可能小。 在这里，我们有 $u$ 的二阶空间导数，可以通过应用部分积分技术将其转换为 $u$ 和 $v$

$-\int_{\Omega}\left(\nabla^{2} u\right) v \mathrm{~d} x=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x-\int_{\partial \Omega} \frac{\partial u}{\partial n} v \mathrm{~d} s,\qquad(5)$

其中 $\frac{\partial u}{\partial n}=\nabla u \cdot n$ 是$u$ 在边界外法线方向$n$

变分公式的另一个特点是测试函数 $v$ 需要在解 $u$ 已知的边界部分消失。 在当前问题中，这意味着整个边界 $\partial \Omega$ 上的 $v=0$。 因此，等式（5）右边的第二项消失了。 从等式（4）和（5）可以得出

$\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x,\qquad(6)$

抽象有限元变分公式

事实证明，为变分问题引入以下规范符号是很方便的：找到 $u \in V$

$a(u, v)=L(v) \quad \forall v \in \hat{V},\qquad(9)$

对于泊松方程，我们有：

$a(u, v)=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x,\qquad(10)$

$L(v)=\int_{\Omega} f v \mathrm{~d} x,\qquad(11)$

实现代码（Python）

from fenics import *
# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, ’P’, 1)
# Define boundary condition
u_D = Expression(’1 + x[0]*x[0] + 2*x[1]*x[1]’, degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bc)
...

详情参阅 - 亚图跨际