Julia ModelingToolkit含Dx(var)*Dx(var)项的PDE求解方法问询

1. 优先保留方程的散度形式，避免手动展开导数乘积 #

电迁移的质量守恒方程核心是 $\frac{\partial c_i}{\partial t} = -\nabla \cdot N_i$，代入通量公式后不要手动展开散度项，让ModelingToolkit自动处理符号运算。直接展开导数乘积容易引入人为错误，还会导致求解器误识别额外依赖变量。示例代码：

using ModelingToolkit, DifferentialEquations

@variables t x
@variables c_cat(..) Φ₂(..)
Dt = Differential(t)
Dx = Differential(x)

# 定义参数
@parameters D_cat z_cat u_cat F

# 构建通量与连续性方程（直接用散度，不手动展开乘积）
N_cat = -D_cat * Dx(c_cat(t,x)) - z_cat * u_cat * F * c_cat(t,x) * Dx(Φ₂(t,x))
eq = Dt(c_cat(t,x)) ~ -Dx(N_cat)

2. 检查变量定义与系统构建 #

• 确保c_cat(t,x)和Φ₂(t,x)都被正确声明为依赖变量，避免将导数项误定义为独立变量。
• 构建PDESystem时，明确列出所有依赖变量和参数，防止求解器遗漏或误识别变量：

sys = PDESystem(
    [eq], 
    [],  # 边界条件后续补充
    [t,x], 
    [c_cat(t,x), Φ₂(t,x)], 
    [D_cat=>1.0, z_cat=>1, u_cat=>0.1, F=>96485]
)

3. 匹配物理合理的边界与初始条件 #

• 电迁移问题中，$\Phi_2$的边界条件需满足电荷守恒，各物种的通量边界条件要和整体电荷平衡一致，避免出现物理矛盾的约束。
• 初始条件设置非奇异的连续值，防止初始时刻导数项出现无穷大或NaN，导致求解器发散。

4. 选择适配非线性PDE的求解器 #

含导数乘积的PDE属于非线性问题，优先使用刚性求解器搭配MethodOfLines离散化：

discretization = MOLFiniteDifference([x=>0:0.1:10], t)
prob = discretize(sys, discretization)
sol = solve(prob, TRBDF2(); abstol=1e-8, reltol=1e-8)

调整abstol和reltol容差，可改善数值不稳定导致的非物理解。

5. 调试变量依赖关系 #

通过以下代码检查系统的状态变量，确认没有额外的未定义依赖变量：

println(ModelingToolkit.get_states(sys))

若出现异常变量，检查方程是否有拼写错误或未声明的符号。

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