1

I am trying to solve numerically a problem in orbital motion, using Runge-Kutta 4 method with adaptive step size. Because energy is the most obvious theoretically conservative number in the motion, I want to dymically control the step size of the calculaion such that $m=\frac{\delta E_T}{E_0}$ (the relative error of the energy at the end of the calculation and the starting, analytical energy, the desired accuracy) is equal or smaller than some number.

I tried to follow this logic: first, I calculated the new state of the particle with some $dt$ (starting or given in previous iteration), use it to calculate the current energy and $s=\frac{\delta E_t}{E_0}$ (relative error in current step) and get a new timestep: $dt_{new} = dt(\frac{m}{s})^{0.2}$ (scaling the timestep as we have seen in class), and use the new timestep to calculate the state of the particle, and advance to next iteration. However, calculation attemps so far have been unsuccessful. I noticed that no matter how small initial $dt$ I give, the relative error eventually grows too large and in turn $dt$ collapses to extremly small numbers, and the calculation effectively becomes an infinite loop.

What did I do wrong? Thanks in advance

Frogfire
  • 303
  • 3
  • 9
  • This does not answer you question, but maybe you are interested: There is a different class of numerical integration schemes called "symplectic integrators" (the simplest instance of it is "leapfrog integration"), that exactly preserves the total energy of the system. Thus it might very well be the better for orbital mechanics simulations than Runge-Kutta. – Simon May 02 '21 at 20:12
  • 1
    @Simon : This will not help here, as the energy preservation of symplectic methods breaks down at singularities of the Hamiltonian. What works perfectly well with a pendulum fails if you get too close to the center of a gravity well. Or in other words, the closer the orbit gets to the center, the higher the speed, the shorter the step size needs to be to have enough sampling points in that critical segment. In my answer to the related https://physics.stackexchange.com/questions/633293/numerical-errors I have some estimates for RK4, that should remain relevant for any other explicit method. – Lutz Lehmann May 03 '21 at 16:48

1 Answers1

1

One observation is that you compare the local error over $dt$ to the (desired) global error over time $T$. This makes not much sense. What you want is that the local error is proportional to the global error, roughly $s:m = dt:T$. So you need to change the formula to $$ dt_{\rm new}=dt·\left(\frac{m·dt}{s}\right)^{0.2} ~~\text{ or fully }~~ dt_{\rm new}=dt·\left(\frac{m·dt}{s·T}\right)^{0.2}. $$ The exponent does only give a guess of the immediate optimal step size if the error in the energy has global order 5. However, as a regulation process it will tend towards an optimal step size over a small number of steps even if order and exponent are not that closely related.

Lutz Lehmann
  • 131,652
  • Hello, thanks for the reply. I tried to implement the new $dt$ , but the problem still remains (it gets small very fast until computer rounds to 0). Are you sure that $dt_{new}$ should be proportional to $dt^{1.2}$? if $dt$ is small than this power law would make it even smaller... whatever the case, I am sure this is not the main problem because $dt$ falls to 0 too fast to be explained by this alone – Frogfire May 02 '21 at 19:48
  • There is nothing much more to say with the present information. It could be that the method step has an implementation error, or that the energy function is not correct, or that some typo in the ODE functions makes this non-Hamiltonian. – Lutz Lehmann May 02 '21 at 20:02