Python中有限差分积分结果异常，与RK4差异极大的原因排查

问题 #

我需要求解函数S对应的微分方程，先用四阶龙格-库塔（RK4）方法实现积分，代码如下：

def rk4_step(r, u, h, f, *args):
    k1 = h * f(r, u, *args)
    k2 = h * f(r + h/2, u + k1/2, *args)
    k3 = h * f(r + h/2, u + k2/2, *args)
    k4 = h * f(r + h, u + k3, *args)
    u_new = u + (k1 + 2.0 * k2 + 2.0 * k3 + k4)/6.0
    return u_new

随后用该方法计算得到S[i]的解，结果符合预期：

S[0] = 0
S[1] = 0
for i in range(1, nr_int):
     S[i+1] = rk4_step(r[i], S[i], dr, rhs_H, index, l, 0, 0, 0, H[i], lamb[i], m[i], nu[i], p[i], rho[i], cs2[i])
for i in range(nr_int, nr - 1):
     S[i+1] = rk4s_step(r[i], S[i], dr, rhs_H, index, l, 0, 0, 0, 0, lamb[i], m[i], nu[i], 0, 0, 1)

之后尝试改用有限差分法，基于公式：

S[i] = ( S[i+1] - S[i-1] ) / ( 2 * dr )

编写代码如下：

S[0] = 0
S[1] = 0
for i in range(1, nr_int):
     S[i+1] = S[i-1] + 2 * dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, H[i], lamb[i], m[i], nu[i], p[i], rho[i], cs2[i])
for i in range(nr_int, nr - 1):
     S[i+1] = S[i-1] + 2 * dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, 0, lamb[i], m[i], nu[i], 0, 0, 1)

但运行结果与RK4结果相差约300个数量级，完全错误，请问这种有限差分方法存在什么问题？

问题分析 #

你的有限差分实现存在两个核心问题：

1. 格式稳定性不匹配初值问题 #

你用的是显式中心差分（蛙跳格式），这种多步格式对于常微分方程（ODE）的初值问题是绝对不稳定的：

• RK4是单步方法，每一步只依赖前一个点的信息，稳定性良好，误差不会指数级放大；
• 而显式中心差分每一步依赖前两个点的信息，其数值增长因子的模大于1，每计算一步误差都会被放大，经过多次迭代后，误差会呈指数级爆发，最终导致结果完全偏离真实解。

2. 微分方程的场景适配错误 #

你的原微分方程应为 dS/dr = rhs_H(...)，你写的公式是用中心差分近似S的导数，对应关系本身的代数推导是对的，但这个格式只适合边值问题的求解（比如两点边值问题的差分法），完全不适合你的初值问题场景——初值问题需要从初始点向前递推，这种多步不稳定格式根本无法支撑这类计算。

修正建议 #

如果想改用有限差分法求解初值问题，应使用一阶显式向前差分（虽然精度低于RK4，但稳定性适配初值问题）：

S[0] = 0
for i in range(0, nr - 1):
     # 根据区间选择对应的rhs_H参数
     if i < nr_int:
         S[i+1] = S[i] + dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, H[i], lamb[i], m[i], nu[i], p[i], rho[i], cs2[i])
     else:
         S[i+1] = S[i] + dr * rhs_H(r[i], S[i], index, l, 0, 0, 0, 0, lamb[i], m[i], nu[i], 0, 0, 1)

或者选择隐式差分格式（稳定性更好，但需要解方程组），但对于初值问题，RK4这类单步方法本身就是更合适的选择。

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