I am a software dude working on a personal project for building a spacecraft simulation software. And I've been looking into numerical integration methods, and I am just confused on how to tackle the dynamics of applying forces especially when applied to numerical integration.

Here is my dilemma, the position and velocity of a spacecraft is given, the time step is h, the forces (f) applied on the spacecraft is then calculated (gravity, drag, etc.). Next is we integrate the position (pos = vel) and velocity (accel = force/mass). Now I'm using an RK4 integrator, (1st case) we assume that the forces are uniform all throughout one integration time step. (is this correct?) . (2nd case) the forces should be evaluated inside the integration step (for example, assume k1 to k4 are the values we get throughout RK4 integration, during k2 where t+h/2, use the new position and vel to get a new force, also f2, continue to next step to get k3, evaluate f3, get the new k4 from this and then finally get the new pos and vel from the values k1, k2, k3, k4).

When looking at this, intuitively, 2nd case should be the most accurate right? but case 1 should also be accurate especially for smaller time steps. But now I'm at the point where I am dealing with impulses (e.g. propulsion systems, etc.) and now its getting me confused about this whole dynamics. How can I apply an impulse to my simulations? Does the 1st case assume an impulse is applied and you are just predicting the next time step, and 2nd case was wrong. Or is the 2nd case correct so that when i want to apply an impulse, i should only consider adding force for k1, and zero the forces at k2, k3, k4.

Could anyone please help clarify these stuff for me?

I've been reading through a lot of astrodynamics, simulations, and spacecraft stuff but i still couldn't find my answer and I cant sleep well thinking about it LOL. (Fundamentals of Astrodynamics by David Vallado, Orbital Mechanics for Engineering Students by Howard Curtis, Analytical Mechanics for Space Systems by Schaub, etc.)

  • 1
    $\begingroup$ Propulsion forces are just forces. You treat them like any other force. If they are highly time variant, reduce the step size. In the Shuttle Mission Simulator the propulsion system models (the parts of them that calculated forces) and the dynamics models ran 25 times per second. Other models ran slower, down to as slow as a 1/2 Hz. The SMS also used an RK4 integrator. $\endgroup$ Commented Mar 14 at 1:40
  • 1
    $\begingroup$ "...dealing with impulses (e.g. propulsion systems, etc.)" Assuming that the time step is in (hundreds of) milliseconds, what type of propulsion system are you modelling that switches on and off in a matter of (hundred) milli seconds? Also, is plain RK suitable for stiff ODEs? Usual ODE solvers assume that states involved are smoothly varying. $\endgroup$
    – AJN
    Commented Mar 14 at 13:21

2 Answers 2


Now I'm using an RK4 integrator, (1st case) we assume that the forces are uniform all throughout one integration time step. (is this correct?)

That is not correct.

The point of using a multiple stage integrator such as RK4 is to calculate derivatives (in this case, velocities and accelerations) at the intermediate steps. Assuming the forces are uniform defeats the purpose.

What you should be doing instead is to calculate gravitational acceleration (not force) at the initial stage and at each intermediate stage, mass at the initial stage and at each intermediate stage, thrust (which varies due to thrusters having a build-up and trail-off) at the initial stage and at each intermediate stage, drag at the initial stage and at each intermediate stage, etc.

This means that those intermediate stage derivative function calls come at a computational expense, sometimes a significant expense. The advantage of techniques that take multiple intermediate stages per integration step is that they enable taking larger steps without losing accuracy. However, you shouldn't use RK4 if you are taking one millisecond steps because of flight software or because of high frequency dynamics such as slosh or flex. At such a high frequency, the added computational cost has no benefit with regard to accuracy.

In fact, accuracy tends to improve with larger integration steps, at least up to a point. That inflection point where accuracy and stability degrades with even larger steps depends on the object being integrated (whether the orbit nearly circular, highly elliptical, or hyperbolic). For an object in a nearly circular low Earth orbit, a one second step size using RK4 comes very close to the inflection point that yields the best accuracy for RK4. A lot of space vehicles have flight software that operates at one hertz. This is one of the key reasons RK4 is widely used for LEO spacecraft.

You also asked about impulsive burns. You should treat those as exactly what they are, which is impulsive. Do not treat impulsive burns as forces that extend over a non-zero period of time. Simply change the velocity instantaneously, typically just prior to the initial stage of the RK4 cycle. Note well: Modeling burns as impulsive burns is a low fidelity approach.


Here's a theory I have about how one could implement an impulse with an RK4 integrator.

I believe the formal answer is that the force can be expressed as F(t). The right way to apply an RK4 is to evaluate F(t) at the times associated with the sub intervals when computing k for that sub interval, before combining k to find your final value of the integrated variable (for example the position). So, for example, if you're using the most classic variant of RK4, the times should be evaluated at t0, t0+h/2, t0+h/2, and t0+h, where h is the substep.

An impulse can be formally written as the Dirac delta function, at least in a physics context. I am not totally sure how one would write it in a rocket context but here I am discussing more generally how to integrate it.

F(t) = F0*delta(t-t0), which has the property that at t = t0, the force is infinite, and at t = anything else, it's zero. When integrated in an interval that contains t0, the integral is F0, by definition. The name for the integral of force with respect for time is impulse, so the impulse is F0 if it is integrated over the interval containing t0 at which the impulse occurs.

When computing the k_n for the RK4 substeps within the sub interval, only one of the k_n's will have t0 within its interval. So only one of them will contain an "infinite" value for F(t). But should you really add an infinite value at a single point? That wouldn't be very helpful.

Instead, consider the sum

y' = y + h/6(k1 + 2k2 +2k3+k4)

Suppose you are integrating the function

f(x,t) = G(x,t)+F(x,t) where F is your force

and the interval where t0 occurs is in k2.

Let g1 through g4 be the k's for the g term. The k's for the F term will be zero except for the second term, because the integral of a delta function when t0 is not in that interval is zero. Normally the h accounts for integrating over that sub interval. It's like the base of the rectangle for the area under the curve. I'm perfectly aware that RK4 is not rectangles, more on that later. In this case, the impulse is actually going to replace h*k2, so that term becomes F0. I AM NOT SURE IF THE FACTOR 2/6 ALSO NEEDS TO BE REPLACED.

So, if the constant factor is right,

y' = y0 + h/6(g1 + 2g2 +2 g3 + g4) + 1/3 F0

Question mark on the 1/3. This is if the impulse occurs in the second substep. The factor is 1/6 if it's in the first or fourth, 1/3 if it's in the third. IF THIS SIMPLE APPROACH IS CORRECT.

I would like to add the caveat that every time I have actually worked with a numerical integration routine, forces have been applied over a smoothed interval. So, if an impulse is desired, instead it is convolved with a window function so that the impulse becomes some other very narrow function instead, that is differentiable and possibly also has a continuous second derivative.

To address the fact that the integration of the PDE in the RK4 is not as simple as summing rectangles (the Euler method): The right way to compute the numeric factor in the method I've just proposed is to consider something like f(x+delta x,t+h)-f(x,t) = O(h^4), when the entire process is taken into account. I have not verified that formula with an impulse in F(t) at t0 in the second time step, replacing h k2 for F. I think there could be a constant factor mistake, but dimensionally, this analysis should be correct, I believe.

  • $\begingroup$ I suggest the reference Computational Physics by Mark Newman if you want to derive this. I am currently working, in the sense that I am hoping I will get work any time now, and was not hoping to embark upon a derivation of that nature. But happy to share what I know off the top of my head. $\endgroup$ Commented Mar 14 at 22:25

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.