FEALPy 有限元求解 Poisson 方程
考虑二维 Poisson 方程
−
Δ
u
=
f
-\Delta u = f
−Δu=f
给定一个真解为
u
=
cos
π
x
cos
π
y
u=\cos\pi x \cos\pi y
u=cosπxcosπy
Poisson 方程, 其定义域为
[
0
,
1
]
2
[0,1]^2
[0,1]2, 只有纯的 Dirichlet 边界条件, 下面演示基于 FEALPy 求解这个算例的过程.
- 导入创建 pde 模型.
from fealpy.pde.poisson_2d import CosCosData # 导入二维 Poisson 模型实例
pde = CosCosData() # 创建 pde 模型对象
- 生成网格.
from fealpy.mesh import MeshFactory as MF # 导入网格工厂模块
domain = [0, 1, 0, 1]
mesh = MF.boxmesh2d(domain, nx=10, ny=10, meshtype='tri')
- 建立拉格朗日有限元空间, 代码中
p=1
也可以设为更大正整数, 表示建立 p p p 次元的空间.
# 导入 Lagrange 有限元空间
from fealpy.functionspace import LagrangeFiniteElementSpace
space = LagrangeFiniteElementSpace(mesh, p=1) # 线性元空间
- 建立 Dirichlet 边界条件处理对象.
# 导入 Dirichlet 边界处理
from fealpy.boundarycondition import DirichletBC
bc = DirichletBC(space, pde.dirichlet) # 创建 Dirichlet 边界条件处理对象
- 创建离散代数系统, 并进行边界条件处理.
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
- 求解离散系统.
# 导入稀疏线性代数系统求解函数
from scipy.sparse.linalg import spsolve
uh[:] = spsolve(A, F)
- 计算 L2 和 H1 误差.
L2Error = space.integralalg.error(pde.solution, uh)
H1Error = space.integralalg.error(pde.gradient, uh.grad_value)
- 画出解函数和网格图像
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>