避免 Runge-Kutta 4(5) 求解器的 ZeroDivisionError

Question

我正在尝试创建一个 Runge-Kutta 4(5) 求解器来求解初始条件为

y' = 2t

的微分方程

y(0) = 0.5

。这是我到目前为止所拥有的：

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 

这是因为四阶解

u1

和五阶解

u2

非常接近，以至于它们的差异基本上为零，当我计算

delta

时，我得到

1/0

，这显然会导致 ZeroDivisionError。

解决这个问题的一种方法是不计算

delta

并使用更简单的 RK45 版本：

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)

但这虽然有效，但对我来说似乎毫无意义，因为与 RK4 方法相比，它否定了 RK45 方法的自适应步长优势。

有没有办法在不遇到 ZeroDivisionErrors 的情况下保持自适应步长？

Answer 1

您可以尝试捕获错误然后在触发时处理它或使用小数以获得更高的精度。小数中的数字存储为元组，而浮点数是数字的近似值。当高精度很重要时，建议使用小数。

参见例如Error python：[ZeroDivisionError: division by zero] 和浮点任意精度可用吗？

和 decimal 模块文档：decimal — 十进制定点和浮点运算