如何在Python中手动实现隐式数值方法求解刚性微分方程？

首先需要澄清：你写的bdf2_step代码结构其实是隐式Runge-Kutta（IRK）（比如Radau IIA 3阶方法），而非向后微分公式BDF2（BDF是线性多步法，依赖前n步的解）。隐式方法的核心问题是需要求解非线性方程组——因为k1/k2/k3同时出现在等式的两边，无法像显式RK那样顺序计算，必须通过迭代法求解。

方案1：修正隐式Runge-Kutta（Radau IIA）的实现 #

我们可以用非线性方程组求解器（比如牛顿法）来计算k1/k2/k3，具体步骤如下：

代码实现 #

import numpy as np
from scipy.optimize import root

def irk_radau_step(r, u, h, f, *args):
    # Radau IIA 3阶方法的系数（你的原始参数）
    sqrt6 = np.sqrt(6)
    a = np.array([
        [(88-7*sqrt6)/360, (296-169*sqrt6)/1800, (-2+3*sqrt6)/225],
        [(296+169*sqrt6)/1800, (88+7*sqrt6)/360, (-2-3*sqrt6)/225],
        [(16-sqrt6)/36, (16+sqrt6)/36, 1/9]
    ])
    b = np.array([(16-sqrt6)/36, (16+sqrt6)/36, 1/9])
    c = np.array([(4-sqrt6)/10, (2+sqrt6)/10, 1])
    
    # 定义残差函数：res(k) = k - h*f(r + c*h, u + a@k)
    def residual(k):
        k_col = k.reshape(-1, 1)  # 转为列向量方便矩阵运算
        u_k = u + a @ k_col
        # 计算每个节点的f值
        f_vals = np.array([
            f(r + ci*h, uk.flatten(), *args) 
            for ci, uk in zip(c, u_k)
        ]).reshape(-1, 1)
        return (k_col - h * f_vals).flatten()
    
    # 用显式RK4的结果作为初始猜测，提高求解成功率
    k1_guess = h * f(r, u, *args)
    k2_guess = h * f(r + h/2, u + k1_guess/2, *args)
    k3_guess = h * f(r + h/2, u + k2_guess/2, *args)
    k_initial = np.concatenate([k1_guess, k2_guess, k3_guess])
    
    # 求解非线性方程组
    sol = root(residual, k_initial)
    if not sol.success:
        raise RuntimeError(f"隐式RK求解失败: {sol.message}")
    
    # 计算新的u值
    k_sol = sol.x.reshape(-1, 1)
    u_new = u + b @ k_sol
    return u_new.flatten()

方案2：实现真正的BDF2方法（线性多步法） #

BDF2是二阶隐式线性多步法，依赖前两步的解，公式为：
$$\frac{3u_{n+1} - 4u_n + u_{n-1}}{2h} = f(r_{n+1}, u_{n+1})$$
同样需要求解关于$u_{n+1}$的非线性方程，代码实现如下：

代码实现 #

import numpy as np
from scipy.optimize import root

def bdf2_step(r, u_curr, u_prev, h, f, *args):
    r_next = r + h
    
    # 定义残差函数：res(u_new) = 3u_new -4u_curr +u_prev -2h*f(r_next, u_new)
    def residual(u_new):
        return 3*u_new - 4*u_curr + u_prev - 2*h*f(r_next, u_new, *args)
    
    # 用欧拉法结果作为初始猜测
    u_initial = u_curr + h*f(r, u_curr, *args)
    
    # 求解非线性方程
    sol = root(residual, u_initial)
    if not sol.success:
        raise RuntimeError(f"BDF2求解失败: {sol.message}")
    
    return sol.x

BDF2的使用方式 #

因为BDF2是多步法，需要先用显式方法（比如你的RK4）计算第一步的解，再迭代：

def solve_with_bdf2(f, r0, u0, h, num_steps):
    r = r0
    # 初始化前两步：用RK4计算第一步
    u_prev = u0
    u_curr = rk4_step(r, u0, h, f)
    r += h
    results = [(r0, u0), (r, u_curr)]
    
    # 后续步骤用BDF2
    for _ in range(num_steps - 2):
        u_new = bdf2_step(r, u_curr, u_prev, h, f)
        r += h
        u_prev, u_curr = u_curr, u_new
        results.append((r, u_new))
    
    return results

关键注意事项 #

• 隐式方法必须求解非线性方程组：这是你之前代码报错的核心原因，直接赋值k1会导致循环引用，必须通过迭代法求解。
• 初始猜测很重要：用显式方法的结果作为初始值，能大幅提高非线性求解器的收敛速度和成功率。
• 关于scipy solve_ivp的问题：如果要用该函数，必须指定method='BDF'或method='Radau'（默认是RK45，不适合刚性系统），同时调整atol/rtol参数控制误差，避免步长波动。

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