Runge-Kutta 4(5)求解器：保留自适应步长并避免ZeroDivisionError

问题描述 #

我正在实现一个Runge-Kutta 4(5)求解器，用于求解微分方程y' = 2t，初始条件为y(0) = 0.5。现有代码运行时触发ZeroDivisionError，原因是4阶解u1与5阶解u2差异极小（甚至为0），计算步长调整因子delta时出现除以零的情况。简化版RK45虽能正常运行，但完全丧失了RK45特有的自适应步长优势。需要在保留自适应步长特性的同时避免该错误。

初始实现代码 #

import numpy as np

def rk45(f, u0, t0, tf=100000, epsilon=0.00001, debug=False):
    h = 0.002
    u = u0
    t = t0
    # solution array
    u_array = [u0]
    t_array = [t0]
    
    if debug:
        print(f"t0 = {t}, u0 = {u}, h = {h}")
    while t < tf:
        h = min(h, tf-t)
        k1 = h * f(u, t)
        k2 = h * f(u+k1/4, t+h/4)
        k3 = h * f(u+3*k1/32+9*k2/32, t+3*h/8)
        k4 = h * f(u+1932*k1/2197-7200*k2/2197+7296*k3/2197, t+12*h/13)
        k5 = h * f(u+439*k1/216-8*k2+3680*k3/513-845*k4/4104, t+h)
        k6 = h * f(u-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40, t+h/2)
        u1 = u + 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5
        u2 = u + 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55
        R = abs(u1-u2) / h
        print(f"R = {R}")
        delta = 0.84*(epsilon/R) ** (1/4)
        
        if R <= epsilon:
            u_array.append(u1)
            t_array.append(t)
            u = u1
            t += h
        h = delta * h
        if debug:
            print(f"t = {t}, u = {u1}, h = {h}")
    return np.array(u_array), np.array(t_array)

def test_dydx(y, t):
    return 2 * t

initial = 0.5
sol_rk45 = rk45(test_dydx, initial, t0=0, tf=2, debug=True)

运行报错信息 #

t0 = 0, u0 = 0.5, h = 0.002
R = 5.551115123125783e-14
t = 0.002, u = 0.5000039999999999, h = 0.19463199004973464
R = 0.0

---------------------------------------------------------------------------
ZeroDivisionError 

简化版代码（丧失自适应步长优势） #

import numpy as np

def rk45(f, u0, t0, tf=100000, epsilon=0.00001, debug=False):
    h = 0.002
    u = u0
    t = t0
    # solution array
    u_array = [u0]
    t_array = [t0]
    
    if debug:
        print(f"t0 = {t}, u0 = {u}, h = {h}")
    while t < tf:
        h = min(h, tf-t)
        k1 = h * f(u, t)
        k2 = h * f(u+k1/4, t+h/4)
        k3 = h * f(u+3*k1/32+9*k2/32, t+3*h/8)
        k4 = h * f(u+1932*k1/2197-7200*k2/2197+7296*k3/2197, t+12*h/13)
        k5 = h * f(u+439*k1/216-8*k2+3680*k3/513-845*k4/4104, t+h)
        k6 = h * f(u-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40, t+h/2)
        u1 = u + 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5
        u2 = u + 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55
        R = abs(u1-u2) / h
        
        if R <= epsilon:
            u_array.append(u1)
            t_array.append(t)
            u = u1
            t += h
        else:
            h = h / 2
        if debug:
            print(f"t = {t}, u = {u1}, h = {h}")
    return np.array(u_array), np.array(t_array)

解决方案 #

问题根源 #

测试的微分方程y'=2t是线性方程，RK45的4阶解和5阶解在浮点运算精度下几乎完全一致，误差估计R被舍入为0，导致计算delta时触发除零错误。

修改思路 #

1. 给误差估计R设置极小值下限，避免除以零
2. 限制步长调整因子delta的范围，防止步长过大或过小破坏稳定性
3. 保留RK45自适应步长的核心逻辑，仅在临界情况做保护

修改后的完整代码 #

import numpy as np

def rk45(f, u0, t0, tf=100000, epsilon=1e-5, debug=False):
    h = 0.002
    u = u0
    t = t0
    u_array = [u0]
    t_array = [t0]
    
    if debug:
        print(f"t0 = {t}, u0 = {u}, h = {h}")
    while t < tf:
        h = min(h, tf - t)
        # 计算RK45的各个k值
        k1 = h * f(u, t)
        k2 = h * f(u + k1/4, t + h/4)
        k3 = h * f(u + 3*k1/32 + 9*k2/32, t + 3*h/8)
        k4 = h * f(u + 1932*k1/2197 - 7200*k2/2197 + 7296*k3/2197, t + 12*h/13)
        k5 = h * f(u + 439*k1/216 - 8*k2 + 3680*k3/513 - 845*k4/4104, t + h)
        k6 = h * f(u - 8*k1/27 + 2*k2 - 3544*k3/2565 + 1859*k4/4104 - 11*k5/40, t + h/2)
        # 计算4阶和5阶解
        u1 = u + 25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5
        u2 = u + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50 + 2*k6/55
        
        # 关键修改1：给R设置极小值下限，避免除零
        R = max(abs(u1 - u2) / h, 1e-16)
        if debug:
            print(f"R = {R}")
        
        # 关键修改2：限制delta的范围，防止步长突变
        delta = 0.84 * (epsilon / R) ** (1/4)
        delta = max(0.1, min(5.0, delta))  # 步长最多缩小到1/10，最多放大5倍
        
        if R <= epsilon:
            u_array.append(u1)
            t_array.append(t)
            u = u1
            t += h
        
        h = delta * h
        if debug:
            print(f"t = {t}, u = {u1}, h = {h}")
    
    return np.array(u_array), np.array(t_array)

def test_dydx(y, t):
    return 2 * t

initial = 0.5
sol_rk45 = rk45(test_dydx, initial, t0=0, tf=2, debug=True)

关键修改说明 #

• R的最小值限制：用max(abs(u1-u2)/h, 1e-16)确保R不会为0，避免除零错误，同时1e-16远小于默认的epsilon(1e-5)，不会影响正常的误差判断。
• delta的范围限制：通过max(0.1, min(5.0, delta))限制步长调整的幅度，防止步长突然过大导致稳定性问题，或过小导致计算效率低下。
• 保留自适应逻辑：原有根据误差调整步长的核心逻辑完全保留，仅在极端情况做保护，既解决了除零问题，又保留了RK45自适应步长的优势。

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