Runge-Kutta-Fehlberg $4(5)$ method adaptive size $h$ - iterating too much

Question

I tried to solve baryocentric two body problem with Runge-Kutta-Fehlberg with adaptive size method.

I have two differential equations

$$ \frac{d^2x}{dt^2} = - \frac{\mu}{(x^2 + y^2)^{3/2}} $$

$$ \frac{d^2y}{dt^2} = - \frac{\mu}{(x^2 + y^2)^{3/2}} $$

My initial conditions are

$$ \mu = 3.98 x 10^5 km^3/s^2 $$ $$ x_0 = -2500 km $$ $$ y_0 = -5500 km $$ $$ v_{x0} = 7.5 km/s $$ $$ v_{y0} = 0.5km/s $$

Relative Error = $$10^{-2}$$

Thanks to @Lutz Lehmann I fixed the error and edited the code.

import numpy as np
from scipy import integrate, linalg
import matplotlib.pyplot as plt

# Constants
mu          = 3.98 *10**5

lowerBoundT = 0
upperBoundT = 100000

relHata = 10**(-2)
#                   X0          Y0          VX0         VY0
initialValues = [   -2500,      -5500,      7.5,        0.5     ]

def two_body(t, p):
    x, y, vx, vy = p
    return np.array([vx, vy, -mu * x / (x**2 + y**2)**1.5, -mu * y / (x**2 + y**2)**1.5])

def rkf45Vec(fVec,initialT, initialY, tStop, n, relHata):
    tAll = np.array([initialT])
    yAll = np.array([[initialY]])

    h = (tStop - initialT) / n
    y = initialY
    t = initialT

    def calculateK(t, y, h):
        h_new = h
        h_prev = h
        k = 0
        while True:
            k0 = h_new * fVec(t            , y)
            k1 = h_new * fVec(t + 1/4*h    , y + 1/4*k0)
            k2 = h_new * fVec(t + 3/8*h    , y + 3/32*k0 + 9/32*k1)
            k3 = h_new * fVec(t + 12/13*h  , y + 1932/2197*k0 - 7200/2197*k1 + 7296/2197*k2)
            k4 = h_new * fVec(t + h        , y + 439/216*k0 - 8*k1 + 3680/513*k2 - 845/4104*k3)
            k5 = h_new * fVec(t + 1/2*h    , y - 8/27*k0 + 2*k1 - 3544/2565*k2 + 1859/4104*k3 - 11/40*k4)
            k = 25/216*k0 + 1408/2565*k2 + 2197/4104*k3 - 1/5*k4

            # Compute the error estimate and the new step size
            delta = np.abs(k0/360 - 128*k2/4275 - 2197*k3/75240 + k4/50 + 2*k5/55)
            s = 0.84*(relHata/delta)**0.25
            maxDelta = np.max(s)
            h_new = h * np.min(s)

            # Update the solution
            if maxDelta > relHata:
                break
            h_prev = h_new

        return k, h_new, h_prev

    iterationCount = 0
    while t < tStop:
        iterationCount += 1
        k, h_new, h_prev = calculateK(t,y,h, 1)

        y = y + k
        t = t + h_prev

        h = h_new
        tAll = np.append(tAll, t)
        yAll = np.append(yAll, y)
    yAll = yAll.reshape((len(tAll), len(initialY)))
    print(f'rkf45Vec step count={n}  {iterationCount} times iterated')
    return tAll, yAll

rkf45Vec_100 = rkf45Vec(two_body, lowerBoundT, initialValues, upperBoundT, 100, relHata)

My differentational equations

Answer 1

You need to use the appropriate step size formula for the method you use. The error tolerance relHata is a bound or guide for the unit step error. The conversion from the local step error to the unit step error is to divide by $h$. This extrapolates the error from the interval of length $h$ to an interval of length $1$.

Also, for the relative error to be relative, you have to scale it by the size of the component, either the current size or a moving average of it.

Remaining in your code structure, I would organize the step size control code as

while True:
    k4, k5, err = RKF45_offsets(f,t,y,h_new)
    scale = abs(y)+abs(k4)
    rho = relHata*h_new*min(scale/abs(err))
    h_prev = h_new
    h_new = 0.84 * rho**0.25 * h_prev
    if 0.01 < rho < 1.2:
        break
return k4, h_new, h_prev