1
$\begingroup$

I am working in the sphere of orbital simulations, where orbital trajectories are computed using the differential equations describing gravity. Due to the great timescales of orbits, a step size of minutes or even hours is common, that is dt = 60 or dt = 3600. Now here is my problem. Consider the following example differential equation:

def d_dt(v):
    return v

Now consider two steps, one using RK4 and one using Euler's method:

def rk4_step(u0, h=60):
    k1 = h * d_dt(u0)    
    k2 = h * d_dt(u0 + 0.5 * k1)
    k3 = h * d_dt(u0 + 0.5 * k2)
    k4 = h * d_dt(u0 + k3)
    return (k1 + 2*(k2 + k3 ) + k4) / 6

def euler_step(u0, h=60):
    return d_dt(u0) * h

Suppose u0 = 20000, which is in the ballpark of spacecraft orbital speeds. Then one step using RK4 yields a result nearly 10,000x the equivalent step using Euler. I believe this is because RK4 multiplies by dt multiple times in intermediary calculations, whereas Euler's method only multiplies by dt once, so the discrepancy grows quickly.

What I could do is to use smaller time-steps - a step size of 1, that is 1 second, yields just 1.7x the equivalent step using RK4 as compared to Euler. However, this would be a computationally expensive method that would also take much longer. Is there a way to decrease the RK4 method's dependence on time step size without doing so? Help would be much appreciated.

$\endgroup$
2
  • 2
    $\begingroup$ In your code, you're solving an ODE whose exact solution grows exponentially (and that has nothing to do with the application you claim to be interested in). You seem to think that the RK4 solution is worse than the Euler solution, but you don't give a reason for that belief. You haven't compared them to the true solution (which is about 2e30 after one step). $\endgroup$ Commented Nov 2, 2023 at 6:01
  • 2
    $\begingroup$ It will be easier to compare methods by 1. Using a test problem that isn't exponentially growing and 2. non-dimensinalizing your problem to avoid huge numbers and scales in intermediate calculations. $\endgroup$
    – whpowell96
    Commented Nov 2, 2023 at 15:34

2 Answers 2

4
$\begingroup$

In the case of a linear ODE, the effective results of the Euler and RK4 method steps are linear and 4th degree polynomials which are the truncated Taylor expansions of the exponential. As is well-known, these diverge rapidly from each other in their growth behavior, with the exponential having much larger values than the quartic polynomial and this one much larger values than the linear function.

The use of Euler in general should remain educational, it is a test if you understood numerical integration.

With RK4 you get decent results already with about 60 to 100 points per full orbit. Now experiment with higher eccentricity of the orbit to find that the fast parts at the periastrum need a much higher point density to not derail the numerical orbit, while the slow parts at the apoastrum would benefit from a much larger step size.

So a variable-step integrator is advised, and probably additionally a strategy to have independent step sizes for different stellar objects. There exists a well-used "Pleiades" test case that has several close encounters in its observed interval, where such a decoupling would visibly make sense.

$\endgroup$
-3
$\begingroup$

First of all, I'm pretty sure your code is buggy. that said, don't use rk4. use an actually good integrator (e.g. tsit5/vern7) with automatic error control and optionally automatic stiffness detection.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.