0
$\begingroup$

I'm trying to implement a C++ code for the evaluation of the solution of an N-body system of ODE. I've started with a 2-body problem just to set the methods properly. I'm comparing the results of a 4 stages Runge-Kutta (RK4), a 4 steps Adams-Bashforth (AB4) and velocity Verlet method. The solutions I get have the correct behaviour for every methods I used, but once I check the conservation of energy for the Verlet method I obtain the following

enter image description here

As far as I understood from the article by Hairer, Lubich and Wanner, Stormer-Verlet method doesn't preserve the energy completely, and this should also apply to velocity Verlet. But I also know that Verlet should be particularly suited to this kind of problems, though I get these behaviours for the RK4 and AB4 energy conservation

enter image description here

enter image description here

which seems to me better than Verlet, they doesn't constantly grow with time at least. So is there something I'm missing about how Verlet method works or could it be just some code problem?

EDIT 1

I have some update and I answer here to Wolfgang Bangerth too. Every picture above results from a time step of $h = 0.1$, using $1000000$ iterations. I realized there was a problem in the code by the way. Now the energy for Verlet method is

enter image description here

which is the same as the previous one but it oscillates more. This oscillation reduce by decreasing the timestep (I tried with time step $h = 0.1, 0.05, 0.025$, related to a fixed number of iterations $N = 1000000$ and about $23, 12, 6$ orbits around the central body respectively) as well as with the code I used before and the slope is almost the same in both cases (about $3 \times 10^{-12}$). Both this picture and the previous one result from $h = 0.1$ for instance.

EDIT 2

Ok, now Verlet seems working properly. Even if the oscillation is clearly not centered on the initial energy and there's also a little modulation which seems to slowly grow with time. I also checked the angular momentum, which should be conserved even better than the energy by Verlet, and I think I can consider it as such (it actually oscillates a bit with a negative slope of order $10^{-21}$).

enter image description here

enter image description here

$\endgroup$
8
  • $\begingroup$ How do the energy plots look if you put them in the same graph? $\endgroup$ – nicoguaro Dec 24 '20 at 17:21
  • 2
    $\begingroup$ @nicoguaro It seems clear from the scales: the 2nd and 3rd ones will look like a flat line if plotted together with the 1st. $\endgroup$ – Federico Poloni Dec 24 '20 at 17:39
  • $\begingroup$ @FedericoPoloni, yes. I just realized that the changes for the last 2 are in the 9th figure. $\endgroup$ – nicoguaro Dec 24 '20 at 17:53
  • 3
    $\begingroup$ How many sub-steps are you taking with the various methods? Asymptotically, for time step zero, all of these methods will converge energy. So the question is how the effort you put in compares to the accuracy you get. It would be interesting to see the Verlet method with four times as many time steps compared to the RK4 method, for example. (Though it is clear that the RK4 will likely win this competition given the graphs you show.) $\endgroup$ – Wolfgang Bangerth Dec 27 '20 at 21:39
  • 2
    $\begingroup$ If the first six digits are constant, then from a numerical perspective the quantity is conserved :-) $\endgroup$ – Wolfgang Bangerth Dec 29 '20 at 5:51

Your Answer

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

Browse other questions tagged or ask your own question.