Python数值和偏微分方程解

Python数值分析

线性方程组

使用numpy.linalg.solve求解以下方程式

$$\begin{array}{rl}4x_1+3x_2-5x_3& =2\\ -2x_1-4x_2+5x_3& =5\\ 8x_1+8x_2& =-3\end{array}$$

import numpy as np

A = np.array([[4, 3, -5], 
              [-2, -4, 5], 
              [8, 8, 0]])
y = np.array([2, 5, -3])

x = np.linalg.solve(A, y)
print(x)

输出：

[ 2.20833333 -2.58333333 -0.18333333]

手工计算得出的结果与上述方法相同。 在后台，求解器实际上是在进行LU分解以获得结果。 您可以检查该函数的帮助，它需要输入矩阵为正方形且为全秩，即所有行（或等效地，列）必须线性独立。

片段14

尝试使用矩阵求逆法求解上述方程

A_inv = np.linalg.inv(A)

x = np.dot(A_inv, y)
print(x)

输出：

[ 2.20833333 -2.58333333 -0.18333333]

我们还可以使用scipy包获得用于LU分解的$L$和$U$矩阵

片段15

获取上述矩阵A的$L$和$U$

from scipy.linalg import lu

P, L, U = lu(A)
print('P:\n', P)
print('L:\n', L)
print('U:\n', U)
print('LU:\n',np.dot(L, U))

输出：

P:
 [[0. 0. 1.]
 [0. 1. 0.]
 [1. 0. 0.]]
L:
 [[ 1.    0.    0.  ]
 [-0.25  1.    0.  ]
 [ 0.5   0.5   1.  ]]
U:
 [[ 8.   8.   0. ]
 [ 0.  -2.   5. ]
 [ 0.   0.  -7.5]]
LU:
 [[ 8.  8.  0.]
 [-2. -4.  5.]
 [ 4.  3. -5.]]

我们可以看到我们得到的$L$和$U$与我们在上述中手工得到的有所不同。 您还将看到LU函数返回的置换矩阵$P$。 此置换矩阵记录了我们如何更改方程式的顺序，以简化计算目的（例如，如果第一行中的第一个元素为零，则它不能成为枢轴方程，因为您不能将其他行中的第一个元素转换为 零。因此，我们需要切换方程式的顺序以获得新的枢轴方程式）。 如果将$P$与$A$相乘，您会发现此置换矩阵会逆转这种情况下方程的顺序。

将$P$和$A$相乘，看看置换矩阵对$A$有什么影响

print(np.dot(P, A))

输出：

[[ 8.  8.  0.]
 [-2. -4.  5.]
 [ 4.  3. -5.]]

插值和曲线拟合 | 方程根 | 微分 | 积分 | 初始值问题 | 两边值问题 | 对称矩阵特征值问题 | 优化

Python解偏微分方程

求解泊松方程

数学问题公式化：

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

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

其中，$u=u(\mathbf{x})$ 为未知函数，$f=f(\mathbf{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),(3)$$

未知数 $\mathbf{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 }fv\mathrm{d}x,(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,(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 }fv\mathrm{d}x,(6)$$

抽象有限元变分公式

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

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

对于泊松方程，我们有：

$$a(u,v)={\int }_{\Omega }\nabla u\cdot \nabla v\mathrm{d}x,(10)$$

$$L(v)={\int }_{\Omega }fv\mathrm{d}x,(11)$$

实现代码（Python）

# 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)
# Plot solution and mesh
plot(u)
plot(mesh)
# Save solution to file in VTK format
vtkfile = File(’poisson/solution.pvd’)
vtkfile << u
# Compute error in L2 norm
error_L2 = errornorm(u_D, u, ’L2’)
# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)

import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))
# Print errors
print(’error_L2 =’, error_L2)
print(’error_max =’, error_max)
# Hold plot
interactive()

有限元求解器库 | 子域和边界条件 | 改进泊松求解器

源代码

详情参阅 - 亚图跨际