Avoiding ZeroDivisionError for Runge-Kutta 4(5) solver

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?

First variant

There is no optimal, "first principles" step size selection strategy, so to have a reasonable one has to be sufficient. So as a first improvement make a division by zero impossible, use in the delta computation

(epsilon/(1e-2*epsilon + R))

This will put in a ceiling in the step factor.

Second variant

You have a treatment for when the error is larger than permitted.

However, what is missing is a treatment for when the error is drastically smaller than the error tolerance. So insert a line

if 32*R < epsilon: h *= 2

at the end of the "accept" block. This is for instance important when the solution converges towards an equilibrium as then you want to increase the step size (up to the stability border for explicit methods).

PS

To reduce the amount or severity of catastrophic cancellation you might want to compute the error using the linear combination of the k's with the differences of the step coefficients.

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