5
$\begingroup$

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:

Altitude as a function of time Altitude as a function of time

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.

Trajectory propaged using RK4 for different step sizes

Improve this question
$\endgroup$
5
$\begingroup$

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.

Improve this answer
$\endgroup$
8
  • $\begingroup$ I performed several simulation runs with the same initial conditions, but using RK4 instead for several different values of the time. The results are shown in this figure. As you can see, the results are still completely different from the runs with RKF7(8), and there is also no sign of convergence. I could try to reduce the step size even further, but the run with the smallest step size already took well over an hour to complete. $\endgroup$ – theflyingaviator Oct 19 '15 at 8:56
  • $\begingroup$ With respect to the tolerances, I can both set absolute and relative tolerances. I initially tried the small (e-10, e-15) values for the relative tolerances and 1.0 and 0.1 meters for the position and velocity respectively. I then decreased this value for the absolute tolerance as it was giving poor results. $\endgroup$ – theflyingaviator Oct 19 '15 at 9:02
  • $\begingroup$ The way I modify the step size is as follows. The integrator has a function to determine the next step size based on the provided tolerances. As this step size exceeds a certain interval (let's say 1 sec.), I reduce the time step such that the next time will have an integer value (eg, 156.0 s, 157.0 s. etc.). By reducing this time step, I artificially decrease the tolerance, which should by definition result in a more accurate solution and should therefore not be any different from the original solution. I rather not use the interpolation, since this might introduce rounding errors (possibly) $\endgroup$ – theflyingaviator Oct 19 '15 at 9:06
  • $\begingroup$ @Propuls1on See the edit to my answer. You should also probably edit your question to include that new plot so others can easily see the new information you have. $\endgroup$ – Doug Lipinski Oct 19 '15 at 11:21
  • $\begingroup$ @Propuls1on: If you don't see convergence, then there are only two options: (i) the time step is still not small enough, regardless of the fact that your computation takes awkwardly long, (ii) there is a bug somewhere. $\endgroup$ – Wolfgang Bangerth Oct 19 '15 at 13:10
5
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ As far as I understood it, the integrator implementation has two different tolerances: the absolute and relative tolerance. These tolerances are used to scale the step size: if the difference between the higher order method (which is assumed to be more accurate and is used as an "equivalent" of the exact solution) and the lower order method is within the tolerances, the step size was appropriate or could even be increased. If the difference is not inside the tolerance, the step size should be decreased. $\endgroup$ – theflyingaviator Oct 19 '15 at 6:45
  • $\begingroup$ The variables that are integrated are contained in the state vector (so X, Y, Z, Xdot, Ydot, Zdot). The position components are in the order 10e6 and the velocity components are in the order 10e3. $\endgroup$ – theflyingaviator Oct 19 '15 at 6:47
2
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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