FEALPy 求解 Poisson 方程

考虑二维 Poisson 方程 \[ -\Delta u = f\] 给定一个真解为 \[u=\cos\pi x \cos\pi y\] Poisson 方程, 其定义域为 \([0,1]^2\), 只有纯的 Dirichlet 边界条件, 下面演示基于 FEALPy 求解这个算例的过程.

1. 导入创建 pde 模型.

from fealpy.pde.poisson_2d import CosCosData # 导入二维 Poisson 模型实例
pde = CosCosData() # 创建 pde 模型对象

2. 生成网格.

from fealpy.mesh import MeshFactory as MF # 导入网格工厂模块
domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=10, ny=10, meshtype='tri')

3. 建立拉格朗日有限元空间, 代码中 p=1 也可以设为更大正整数, 表示建立 \(p\) 次元的空间.

# 导入 Lagrange 有限元空间
from fealpy.functionspace import LagrangeFiniteElementSpace 
space = LagrangeFiniteElementSpace(mesh, p=1)  # 线性元空间

4. 建立 Dirichlet 边界条件处理对象.

# 导入 Dirichlet 边界处理
from fealpy.boundarycondition import DirichletBC 
bc = DirichletBC(space, pde.dirichlet) # 创建 Dirichlet 边界条件处理对象

5. 创建离散代数系统, 并进行边界条件处理.

uh = space.function() # 创建有限元函数对象
A = space.stiff_matrix() # 组装刚度矩阵对象
F = space.source_vector(pde.source) # 组装右端向量对象
A, F = bc.apply(A, F, uh) # 应用 Dirichlet 边界条件

run serial_construct_matrix with time: 0.0063852430000679306

6. 求解离散系统.

# 导入稀疏线性代数系统求解函数
from scipy.sparse.linalg import spsolve
uh[:] = spsolve(A, F)

7. 计算 L2 和 H1 误差.

L2Error = space.integralalg.error(pde.solution, uh)
H1Error = space.integralalg.error(pde.gradient, uh.grad_value)

8. 画出解函数和网格图像

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
axes = fig.add_subplot(1, 2, 1, projection='3d')
uh.add_plot(axes, cmap='rainbow')
axes = fig.add_subplot(1, 2, 2)
mesh.add_plot(axes)

<matplotlib.collections.PolyCollection at 0x12af35cd0>

output_16_1