Avoiding ZeroDivisionError for Runge-Kutta 4(5) solver

Question

I am trying to create a Runge-Kutta 4(5) solver to solve the differential equation y' = 2t with the initial condition y(0) = 0.5. This is what I have so far:

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)

When I run it, I get this:

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 

This is because the 4th order solution u1 and 5th order solution u2 are so close together that their difference is essentially zero, and when I calculate delta I get 1/0 which obviously results in a ZeroDivisionError.

One way to solve this is to not calculate delta and use a much simpler version of 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)

But this, while it works, seems incredibly pointless to me because it negates the adaptive step size advantage of a RK45 method as compared to an RK4 method.

Is there any way to preserve an adaptive step size while not running into ZeroDivisionErrors?

Answer 1

You can try to catch the error and then handle it when triggered or use decimals for higher precision. The digits in decimals are stored as tuples vs floats being an approximation of a number. When high precision is important, it's recommended to use decimals.

---

try/except

try:
    z = x / y
except ZeroDivisionError:
    z = 0

---

decimal

from decimal import *

getcontext().prec = 6
Decimal(1) / Decimal(7)

Output: Decimal('0.142857')

from decimal import *

getcontext().prec = 28
Decimal(1) / Decimal(7)

Output: Decimal('0.1428571428571428571428571429')

---

For examples see Error python : [ZeroDivisionError: division by zero] and Is floating point arbitrary precision available?

And decimal module docs: decimal — Decimal fixed point and floating point arithmetic