Numerical integration using RKF7(8) - different results

Question

For my thesis, I look in trajectories of vehicles through an atmosphere at very high velocities. I have a set of equations of motion, which I propagate using a Runge-Kutta-Fehlberg (RKF) 7(8) numerical integrator, one that has already been developed by the department of my university. The entire code is programmed in C++.

For my vehicle, I need to find a specific bank angle profile that achieves a specific objective (it's not necessary to understand what a bank angle profile is, just that I need one as an input to my equations of motion). The optimization algorithm provides a (random) bank angle profile, which I use as an input to my simulation. Based on the result of that simulation, the optimization algorithm changes the bank angle profile and evaluates the new trajectory. This keeps on repeating until a bank angle profile is found that achieves the set objective. Let's call the corresponding trajectory that results from the bank angle profile the reference trajectory.

Now I also am developing a guidance algorithm. This algorithm requires the different state variables and derived values (such as altitude or atmospheric density) at a specified interval, let's say ever 1.0 seconds. I want to apply this guidance algorithm to the found reference trajectory. Since the RKF7(8) algorithm has adaptive step size control, which is beneficial during the optimization since it reduces the number of function evaluations, it is most certainly not guaranteed that an output at every 1.0 seconds is provided. So, to achieve this output at ever 1.0 seconds, the integration is performed again with the same initial conditions as the reference trajectory, expect that the step size obtained from the RKF7(8) is modified to ensure the output at every 1.0 seconds.

The problem, however, is that the obtained trajectory in case the step size is modified is different from the reference trajectory. My guess was that since I reduce the step size to achieve the output at every 1.0 seconds, I force the accuracy to increase due to the smaller step size. I also tried modifying the tolerances used in the RKF7(8) algorithm, as well as different order methods (RKF4(5), RKF5(6) and Adams-Bashfort-Moulton). I decreased this tolerance (both absolute and relative) to 1e-15 (which is I think more or less the limit since it's close to machine precision), but I still obtain different trajectories. Same thi To give you an impression of the difference, see the figure below:

One thing that might have an influence, but I'm not sure about this, is the high velocity. If I reduce the velocity, there is still a small difference but this does not lead to a noticeable difference in the trajectory. However, for my research, it is essential to have these high velocities.

Does anyone has an explanation for the behavior I showed? If so, what could I do to solve this problem?

Edit Based on the comment of Doug Lipinski, I performed several runs with a fixed step-size integrator, Runge-Kutta 4, for a number of step sizes. See the figure below. From this figure, it appears that the solution does not converge to one solution even though the step sizes is very small.

Answer 1

I think the first step is to confirm which of the solutions is more accurate. If you are using a reference implementation of RKF7(8) that has presumably been validated by others on other problems it seems extremely likely that the error is on your end when you're modifying the time step.

If it is possible (i.e. not too computationally expensive). I would try performing a convergence study with a simple fixed time step integration scheme (say RK4). If refining the time step leads to the fixed time step method converging to one of the two solutions you have then that gives additional verification of the correct solution. If you find that the "reference solution" is incorrect then you need to look at how you're setting your tolerance. Ideally you should be using a relative tolerance, not absolute (see Wolfgang's answer).

Assuming the error occurs somewhere in your modifications to the time step I would suggest a slightly different approach: rather than modifying the time step, interpolate the results to the desired time values. If the plot you have shown is typical of the results a simple polynomial interpolation using nearby points should be very accurate.

Edit:

If successive runs at smaller and smaller fixed step size do not appear to be converging, it may be that your problem is stiff or that the problem exhibits signs of chaos (e.g. sensitive dependence on initial conditions, or in this case small differences in each time step). In the first case you could look at using a different time integrator that is appropriate for stiff ODEs (implicit methods are often used). In the second case I'm not really sure what the best option would be.

Answer 2

If you get qualitatively different results from two ODE integrators, then the time step choice of at least one of them is too large. Which one that is is not immediately obvious to say, but if you play with the time step of an integrator (e.g., making it smaller by a factor of ten) and you get a different result, then the time step was too large; if the results are the same then the time step was likely appropriate.

I would continue to play with the tolerance. 1e-15 is not a priori close to the round-off level. The round-off is around 1e-16, but only when seen relative to the size of the variables you are comparing against. If the variables you are comparing against are on the order 1e-12, then 1e-15 is an entirely inappropriately large tolerance. It's also not clear what the tolerance actually compares to -- the error? the residual? Not knowing what it is in the implementation you use, I would continue to play with it to see what happens.

Answer 3

One remark (seeing that the two answers given thus far come with quite a few comments, I've decided to post this as a separate answer instead):

After you've solved the problem with your ODE integrator, you may wish to consider a Runge-Kutta dense output mechanism for the solution interpolation at $t_k = k$ [sec]. A dense output scheme may yield significant compute-time savings compared to your approach of re-running your integration with smaller time steps, particularly if you need to repeat the integration many times (with different initial conditions), which is likely the case for your "guidance" optimization algorithm.

This is because the dense output approximation comes (more-or-less) "for free" - the polynomial solution interpolant can be built directly from the output of the Runge-Kutta stages. And you will not need to artificially decrease your time steps below those computed by the R-K adaptive-step logic (i.e., you don't sacrifice efficiency).

You will find further information on this forum:

Intermediate values (interpolation) after Runge-Kutta calculation

Another excellent reference is the classic book of Hairer and Wanner.