如何使用FiPy声明二维边界条件以复现Matlab偏微分方程组求解示例

核心错误点 #

• 维度理解错误：你要复现的Matlab示例是1D空间+时间的瞬态PDE问题，无需将时间维度硬编码为2D网格的y轴，FiPy的TransientTerm原生支持时间维度演化，你当前用2D网格模拟直接把问题本质搞错了。
• 边界条件设置错误：原问题的边界条件对应空间x的左右两端，你之前既搞反了约束的位置，又误将时间维度的边界当成了空间边界。
• 初始值赋值逻辑错误：双重循环计算cell位置的逻辑不符合FiPy的网格存储规则（FiPy的2D网格cell编号以y为最优先维度遍历），导致初始值赋值完全错误，且setValue原生支持基于网格坐标的条件判断，无需手动生成mask。
• 结果索引逻辑错误：后续提取u1、u2值做可视化的索引逻辑和网格存储规则不匹配，即使计算正确可视化也会出错。

修正实现思路 #

1. 改用1D空间网格，通过时间步迭代获取不同时刻的解，将每步的解保存后拼接为二维数组用于曲面绘制
2. 直接用网格内置坐标属性设置初始值和边界条件，无需手动循环生成mask
3. 边界条件和原问题对齐：
   • x=0（左边界）：u1=1，u2的法向梯度为0
   • x=1（右边界）：u2=0，u1的法向梯度为0

修正后核心代码片段 #

from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm
from fipy.tools import numerix
import numpy as np
import plotly.graph_objects as go

# 配置参数
nx = 100 # 空间x方向网格数
Lx = 1.0 # 空间x范围[0,1]
total_time = 2.0 # 总时间
dt = 0.002 # 时间步长
n_steps = int(total_time / dt)

# 1D空间网格
mesh = Grid1D(nx=nx, Lx=Lx)
u1 = CellVariable(name="u1", hasOld=True, mesh=mesh, value=0.0)
u2 = CellVariable(name="u2", hasOld=True, mesh=mesh, value=0.0)

# 初始值设置，无需循环
u1.setValue(1.0, where=mesh.x < Lx/nx) # x=0处u1初始为1

# 边界条件设置
u1.constrain(1.0, where=mesh.facesLeft)
u2.constrain(0.0, where=mesh.facesRight)
u1.faceGrad.constrain(0.0, where=mesh.facesRight)
u2.faceGrad.constrain(0.0, where=mesh.facesLeft)

# 方程定义和原代码一致
eqn0 = TransientTerm(var=u1) == DiffusionTerm(0.024, var=u1) - (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
eqn1 = TransientTerm(var=u2) == DiffusionTerm(0.170, var=u2) + (numerix.exp(5.73*(u1-u2))-numerix.exp(-11.46*(u1-u2)))
eqn = eqn0 & eqn1

# 存储每一步结果
u1_history = np.zeros((n_steps, nx))
u2_history = np.zeros((n_steps, nx))

# 时间迭代
for t_idx in range(n_steps):
    u1.updateOld()
    u2.updateOld()
    eqn.sweep(dt=dt)
    u1_history[t_idx] = u1.value
    u2_history[t_idx] = u2.value

# 可视化
x = np.linspace(0, Lx, nx)
t = np.linspace(0, total_time, n_steps)

fig = go.Figure(data=[go.Surface(z=u1_history, x=x, y=t)])
fig.update_layout(scene = dict(
                    xaxis_title='X',
                    yaxis_title='T',
                    zaxis_title='u1'),
                    title='u1',
                    width=600,
                    height=600
                  )
fig.show()

fig = go.Figure(data=[go.Surface(z=u2_history, x=x, y=t)])
fig.update_layout(scene = dict(
                    xaxis_title='X',
                    yaxis_title='T',
                    zaxis_title='u2'),
                    title='u2',
                    width=600,
                    height=600
                  )
fig.show()

额外优化说明 #

setValue方法的where参数直接支持布尔类型的numpy数组，你可以直接通过mesh.x、mesh.y获取所有cell的坐标，拼接任意判断条件即可实现精准赋值，完全无需手动遍历计算位置索引。

内容的提问来源于stack exchange，提问作者user16965076