如何在FiPy中实现耦合PDE的x*dA/dx类项

问题背景 #

我是FiPy新手，已查阅FAQ并研究示例，目标是求解一组耦合偏微分方程（PDEs），相关方程如下：

q1R = TransientTerm(var=A1) == -((g/2)* A4).grad[0] +((g/2)*A3).grad[1]- ((g/2)*A6.grad[0]) - ((g/2)*A5.grad[1]) - ((g/2)*(x-(nx*dx)/2)*A4) + ((g/2)*(y-(ny*dy)/2)*A3) - ((g/2)*(x-(nx*dx)/2)*A6) - ((g/2)*(y-(ny*dy)/2)*A5)- k*((x-(nx*dx)/2)*A1.grad[0]+(y-(ny*dy)/2)*A1.grad[1]) 
eq1I = TransientTerm(var=A2) == ((g/2)* A3).grad[0] +((g/2)*A4.grad[1]+ ((g/2)*A5.grad[0]) - ((g/2)*A6.grad[1]) + ((g/2)*(x-(nx*dx)/2)*A3) + ((g/2)*(y-(ny*dy)/2)*A4) - ((g/2)*(x-(nx*dx)/2)*A5) - ((g/2)*(y-(ny*dy)/2)*A6) - k*((x-(nx*dx)/2)*A2.grad[0]+(y-(ny*dy)/2)*A2.grad[1]) 
eq2R = TransientTerm(var=A3) == -((g/2)* A8).grad[0] +((g/2)*A2.grad[0]- ((g/2)*A7.grad[1]) + ((g/2)*A1.grad[1]) - ((g/2)*(x-(nx*dx)/2)*A2) - ((g/2)*(x-(nx*dx)/2)*A8) - ((g/2)*(y-(ny*dy)/2)*A1) - ((g/2)*(y-(ny*dy)/2)*A7)- k*((x-(nx*dx)/2)*A3.grad[0]+(y-(ny*dy)/2)*A3.grad[1])

eq2I = TransientTerm(var=A4) == -((g/2)* A8).grad[1] +((g/2)*A2.grad[1]+ ((g/2)*A7.grad[0]) - ((g/2)*A1.grad[0]) + ((g/2)*(x-(nx*dx)/2)*A1) - ((g/2)*(y-(ny*dy)/2)*A8) - ((g/2)*(y-(ny*dy)/2)*A2) + ((g/2)*(x-(nx*dx)/2)*A7)- k*((x-(nx*dx)/2)*A4.grad[0]+(y-(ny*dy)/2)*A4.grad[1]) 
eq3R = TransientTerm(var=A5) == -((g/2)*(x-(nx*dx)/2)*A2) - ((g/2)*(x-(nx*dx)/2)*A8) + ((g/2)*(y-(ny*dy)/2)*A1) + ((g/2)*(y-(ny*dy)/2)*A7) + ((g/2)*A2).grad[0] - ((g/2)*A1.grad[1]) - ((g/2)* A8).grad[0] + ((g/2)*A7.grad[1])- k*((x-(nx*dx)/2)*A5.grad[0]+(y-(ny*dy)/2)*A5.grad[1]) 
eq3I = TransientTerm(var=A6) == ((g/2)*(x-(nx*dx)/2)*A1) + ((g/2)*(x-(nx*dx)/2)*A7) + ((g/2)*(y-(ny*dy)/2)*A2) + ((g/2)*(y-(ny*dy)/2)*A8) - ((g/2)*A1).grad[0] - ((g/2)*A2.grad[1]) + ((g/2)* A7).grad[0] + ((g/2)*A8.grad[1])- k*((x-(nx*dx)/2)*A6.grad[0]+(y-(ny*dy)/2)*A6.grad[1]) 
eq4R = TransientTerm(var=A7) == - ((g/2)*(x-(nx*dx)/2)*A4) + ((g/2)*(y-(ny*dy)/2)*A3) - ((g/2)*(x-(nx*dx)/2)*A6) - ((g/2)*(y-(ny*dy)/2)*A5) + ((g/2)* A4).grad[0] - ((g/2)*A3.grad[1]) + ((g/2)*A6.grad[0]) + ((g/2)*A5.grad[1])- k*((x-(nx*dx)/2)*A7.grad[0]+(y-(ny*dy)/2)*A7.grad[1])

eq4I = TransientTerm(var=A8) == + ((g/2)*(x-(nx*dx)/2)*A3) + ((g/2)*(y-(ny*dy)/2)*A4) + ((g/2)*(x-(nx*dx)/2)*A5) - ((g/2)*(y-(ny*dy)/2)*A6) - ((g/2)* A3).grad[0] - ((g/2)*A4.grad[1]) - ((g/2)*A5.grad[0]) + ((g/2)*A6.grad[1])- k*((x-(nx*dx)/2)*A8.grad[0]+(y-(ny*dy)/2)*A8.grad[1]) 
eqn = eq1R & eq1I & eq2R & eq2I & eq3R & eq3I & eq4R & eq4I

已将部分变量初始化为CellVariable（高斯初始条件未给出）：

A1 = CellVariable(name="A1", mesh=mesh, hasOld=True)
A2 = CellVariable(name="A2", mesh=mesh, hasOld=True)
A3 = CellVariable(name="A3", mesh=mesh, hasOld=True)
A4 = CellVariable(name="A4", mesh=mesh, hasOld=True)
x = mesh.x
y = mesh.y

通过循环时间步长求解：

for t in range(2000):
    A1.updateOld()
    A2.updateOld()
    A3.updateOld()
    A4.updateOld()

    eqn.solve(dt=1.e-2)

但得到的解与Mathematica验证的预期结果不符，无法定位问题，现整理疑问如下：

疑问 #

• Q1) 如何在A1的方程中实现x*∂A1/∂x、x*∂A2/∂x、x*∂A2/∂y这类项？
• Q2) 如何确保FiPy识别上述表达式中导数里的x与相乘的x为同一变量？
• Q3) 如何在2D网格上定义x和y？

解答 #

Q1：实现x*∂A1/∂x类项 #

在FiPy中，A1.grad[0]对应∂A1/∂x，直接用x * A1.grad[0]即可表示x*∂A1/∂x，同理x * A2.grad[0]对应x*∂A2/∂x，x * A2.grad[1]对应x*∂A2/∂y。

需注意离散化匹配问题：FiPy中梯度是在面（face）上计算的，而x默认是单元格（cell）中心坐标，直接相乘会存在数值维度不匹配。正确做法是转换为面中心坐标：

# 将单元格中心x转换为面中心x
x_face = mesh.x.faceValue
# 用面坐标计算项
term = k * (x_face - (nx*dx)/2) * A1.grad[0]

Q2：确保x变量一致性 #

只要从同一个网格对象获取x和y（即x = mesh.x，y = mesh.y），FiPy会自动识别它们是网格的坐标变量，与梯度计算使用的坐标完全一致。避免手动定义独立的x变量（如x = np.linspace(...)）即可。

可通过打印type(mesh.x)验证，它应为fipy.variables.cellVariable.CellVariable或fipy.variables.faceVariable.FaceVariable，而非普通numpy数组。

Q3：2D网格定义x和y #

首先创建2D网格，再直接提取坐标变量：

from fipy import Grid2D

# 定义网格参数
nx = 50  # x方向单元格数
ny = 50  # y方向单元格数
dx = 0.1  # x方向单元格尺寸
dy = 0.1  # y方向单元格尺寸

# 创建2D网格
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy)

# 获取单元格中心的x、y坐标
x = mesh.x
y = mesh.y

# 若需要面中心坐标
x_face = mesh.x.faceValue
y_face = mesh.y.faceValue

mesh.x和mesh.y会自动返回对应维度的坐标变量，完全适配FiPy的PDE表达式计算。

额外排查建议 #

1. 变量完整性：方程中用到的A5-A8也需定义为CellVariable并设置初始条件，默认值为0会导致解偏差。
2. 时间步长稳定性：尝试减小dt值，或调整求解器参数（FiPy默认隐式求解，但耦合方程可能需更精细设置）。
3. 边界条件匹配：默认零梯度边界，若Mathematica使用了周期性、Dirichlet等其他边界条件，需补充对应设置：

   # 示例：x方向左右边界设为Dirichlet条件
   A1.constrain(0., mesh.facesLeft)
   A1.constrain(0., mesh.facesRight)
   

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