I'm implementing an RKF4(5) integrator, and I can't figure out if my code is working and I don't understand local truncation error, or if my code isn't working.
I apologize for the size of the code block, but the minimum reproducible example is rather large in this case.
import numpy as np
def RKF45(state, derivative_function, h):
"""
Calculate the next state with the 4th-order calculation, along with a 5th-order error
check.
Inputs:
state: the current value of the function, float
derivative_function: A function which takes a state (given as a float)
and returns the derivative of a function at that point
h: step size, float
"""
k1 = h * derivative_function(state)
k2 = h * derivative_function(state + (k1 / 4))
k3 = h * derivative_function(state + (k1 * (3/32)) + (k2 * (9/32)))
k4 = h * derivative_function(state + (k1 * (1932/2197)) + (k2 * (-7200/2197)) + (k3 * (7296/2197)))
k5 = h * derivative_function(state + (k1 * (439/216)) + (k2 * (-8)) + (k3 * (3680/513)) + (k4 * (-845/4104)))
k6 = h * derivative_function(state + (k1 * (-8/27)) + (k2 * (2)) + (k3 * (-3544/2565)) + (k4 * (1859/4104)) + (k5 * (-11/40)))
y1 = state + ((25/216) * k1) + ((1408/2565) * k3) + ((2197/4101) * k4) - ((1/5)*k5)
y2 = state + ((16/135) * k1) + ((6656/12825) * k3) + ((28561/56430) * k4) - ((9/50) * k5) + ((2/55) * k6)
return(y1, y2)
def integrate_RKF45(t0, tmax, tol, h_init, x_0, df, verbose = False):
"""
integrate a function whose derivative is df from t0 to tmax
t0: starting time
tmax: end time
h_init: initial timestep
x_0: starting position
df: a function which takes x and returns the derivative of a function at x
"""
h = h_init
x_i = x_0
t = t0
while t < tmax:
h = min(h, tmax - t)
y1, y2 = RKF45(x_i, df, h)
err_i = np.abs(y1 - y2)
R = 2 * err_i / h
delta = (tol/R)**(1/4)
if verbose:
print(f"t: {t:0.2e}, dt: {h:0.2e}, x: {x_i:0.2e}, err: {err_i:0.2e}")
if err_i < tol:
t += h
x_i = y1
elif err_i > tol:
h *= delta
return(x_i)
def exponential(x_0, k=1):
"""
A simple test function, this returns the input, so it'll integrate to e^x.
"""
return(k * x_0)
if __name__ == "__main__":
integrate_RKF45(t0 = 0.,
tmax = 0.15,
tol = 1e-4,
h_init = 1e-2,
x_0 = 1.,
df = exponential,
verbose=True)
So, this code works to the extent that it returns an approximation of the integral of whatever function I give it. However, the local truncation errors seem way too big. Running the above code will output:
t: 0.00e+00, dt: 1.00e-02, x: 1.00e+00, err: 3.95e-06
t: 1.00e-02, dt: 1.00e-02, x: 1.01e+00, err: 3.99e-06
t: 2.00e-02, dt: 1.00e-02, x: 1.02e+00, err: 4.03e-06
t: 3.00e-02, dt: 1.00e-02, x: 1.03e+00, err: 4.07e-06
t: 4.00e-02, dt: 1.00e-02, x: 1.04e+00, err: 4.11e-06
t: 5.00e-02, dt: 1.00e-02, x: 1.05e+00, err: 4.16e-06
t: 6.00e-02, dt: 1.00e-02, x: 1.06e+00, err: 4.20e-06
t: 7.00e-02, dt: 1.00e-02, x: 1.07e+00, err: 4.24e-06
t: 8.00e-02, dt: 1.00e-02, x: 1.08e+00, err: 4.28e-06
t: 9.00e-02, dt: 1.00e-02, x: 1.09e+00, err: 4.32e-06
t: 1.00e-01, dt: 1.00e-02, x: 1.11e+00, err: 4.37e-06
t: 1.10e-01, dt: 1.00e-02, x: 1.12e+00, err: 4.41e-06
t: 1.20e-01, dt: 1.00e-02, x: 1.13e+00, err: 4.46e-06
t: 1.30e-01, dt: 1.00e-02, x: 1.14e+00, err: 4.50e-06
t: 1.40e-01, dt: 1.00e-02, x: 1.15e+00, err: 4.55e-06
Where the err value is the difference between the 4th- and 5th- order methods. I was under the impression that an n^th-order method had a local truncation error of order O(dt^(n+1)), which means that the above integration should have errors of around 1e-9 instead of 1e-6.
So is my code wrong or is my understanding wrong? Thanks!
