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.