1

I was under the impression that the algorithm should conserve energy if the system being modelled does. I'm modelling the solar system, which should conserve energy. The program conserves angular momentum and does produce stable orbits, but the total energy (kinetic + gravitational potential) oscillates around some baseline. The oscillations are significant. Are there common reasons why this might happen?

Model assumes planets are point masses, circular orbits (I've also tried elliptical orbits and the energy still oscillates) and uses Newtonian mechanics. I can't think what other features of the program might be affecting the outcome.

If it is just expected that the energy oscillates, what causes that??

user13948
  • 443
  • 1
  • 6
  • 14
  • Please provide the source code of your algorithm as well as an example where you think the algorithm fails to generated the expected output... – Willem Van Onsem Jan 07 '16 at 09:45
  • @WillemVanOnsem I'm not sure it is failing. I wish I could upload the graph of the total energy I have, but I have no idea how. I'll add the section with the algorithm, but I'm pretty sure the algorithm itself is all right, if the energy shouldn't oscillate then it's being caused by some other feature/assumption. – user13948 Jan 07 '16 at 09:49
  • 3
    It depends on an algorithm. Basically, if you solve differential equations numerically and you have some motion integrals in your system, then you should use so-called _structure preserving algorithms_ or _symplectic ones_. Other algorithms are dissipative and do not preserve such structures. –  Jan 07 '16 at 09:50
  • @DariuszSendkowski Yes, velocity verlet is symplectic. This is why I wondered about the energy oscillations. – user13948 Jan 07 '16 at 09:51
  • i don't think velocity verlet is entirely symplectic for circular motion (it is for uniformly accelerated one though). if you use stable starting conditions g = v^2/r and evolve a 2d vector by delta_t, there will be a positional deviation from the circle proportional to g*delta_t^2. unless you are using polar coordinates maybe? also, do you take into account only gravitational force from the sun, or from all other planets as well? – Anton Knyazyev Jan 07 '16 at 14:44
  • @AntonKnyazyev Gravitation force from all planets is accounted for. To calculate acceleration I'm using GMm/r^2 and dividing by m. Is that the same as v^2/r? So if I changed the orbits to ellipses, that might help? Not using polar co-ords, thought about it but stuck with Cartesian. – user13948 Jan 07 '16 at 16:05
  • g=v^2/r is just the simplest case when a point orbits a center (in this case it has an analytical solution). i tried plugging it into the velocity verlet formula, and it deviated from the circle even in this case. with a full n-body gravitation force, i'm pretty sure no numerical integrator would be able to produce perfectly non-deviating results. i'm a bit confused how you change orbits from circles to ellipses though. with a full numerical solution, don't the orbits appear naturally, without you explicitly setting them? – Anton Knyazyev Jan 07 '16 at 16:26
  • so i suggest trying a higher-order integrator (such as runge-kutta), and if the energy deviations almost go away (meaning the the calculations are generally correct), you can re-scale the combined kinetic energy to keep the total energy conserved explicitly – Anton Knyazyev Jan 07 '16 at 16:32
  • @AntonKnyazyev The orbits appear naturally given correct initial conditions. So instead of calculating velocity with v = 2*pi*r/T, I will use Nasa's data on maximum velocity, and set the initial position to the closest approach to the sun and initial velocity to the maximum velocity. That should give elliptical orbits, I suspect. – user13948 Jan 07 '16 at 16:44
  • Is there a name for 'smallest distance from sun'? A technical term? – user13948 Jan 07 '16 at 16:45
  • it's perigee i think – Anton Knyazyev Jan 07 '16 at 17:07
  • @AntonKnyazyev That turned up the right search results. Thank you! – user13948 Jan 07 '16 at 17:15
  • btw, do you also simulate how the sun is affected by the planets? – Anton Knyazyev Jan 08 '16 at 16:38
  • @AntonKnyazyev Yes, planets and the sun orbit their common centre of mass. – user13948 Jan 08 '16 at 22:01
  • @AntonKnyazyev I think all of those comments combined amount to a very good answer, if you post it I'll accept it. – user13948 Jan 10 '16 at 09:50

2 Answers2

2

Look up the Verlet-Störmer paper by Hairer et al. (Geometric numerical integration illustrated by the Störmer/Verlet method). There should be several sources online.

In short, a symplectic integrator preserves a Hamiltonian and thus energy, but it is a modified Hamiltonian. If the method is correctly initialized, the modification is a perturbation of order O(h²), where h is the step size. Incorrect initialization gives a perturbation of O(h), while the observed oscillation should still have an amplitude of O(h²).

Thus the observed pattern of an oscillating energy as computed by the physical formulas is completely normal and expected. An error would be observed if the energy were to (rapidly) deviate from this relatively stable pattern.

An easy, but slightly inefficient method to get from the order 2 Verlet method an order 4 symplectic integrator is to replace 

Verlet(h)

by

Verlet4(h) {
  Verlet(b0*h);
  Verlet(-b1*h);
  Verlet(b0*h);
}

where b0=1/(2-2^(1/3))=1.35120719196… and b1=2*b0-1=1.70241438392…. This is called a "composition method".

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
1

Merged from the comments: For a full gravitational N-body problem, I don't think any numerical integrator will be symplectic. Velocity Verlet isn't symplectic even for a single point orbiting a center (easy to check, since it has a trivial analytical solution with g = v^2/R). So I suggest trying a higher-order integrator (such as Runge-Kutta), and if energy deviations almost go away (meaning the the calculations are generally correct), you can re-scale the combined kinetic energy to keep the total energy conserved explicitly. Specifically, you compute the updated Ekin_actual and Ekin_desired = Etotal_initial - Epotential, and scale all velocities by sqrt(Ekin_desired / Ekin_actual)

Anton Knyazyev
  • 464
  • 4
  • 7
  • This is rather wrong. Verlet is always symplectic. What you mean is that the conserved quantities of the discrete dynamic are not the physical ones of the continuous dynamic but rather modifications of them with the step size as parameter. RK4 however is not symplectic, but preserves first integrals to order 4, which is, for "short" integration intervals, more precise than the error order 2 of Verlet. – Lutz Lehmann Feb 06 '16 at 09:05
  • i don't think *any* verlet is symplectic. as the paper you mentioned (which is very good) explains, it's important to have 'staggered' velocity and position updates, i.e. update velocity 'in the middle of the step'. it's mostly known as 'leapfrog' – Anton Knyazyev Feb 06 '16 at 14:31
  • With the correct initialization, leapfrog and velocity Verlet give the same positions, while the velocities for Leapfrog are in the middle. Velocity Verlet is the composition of both of the symplectic Euler method variants, and thus itself symplectic. Remember, symplectic transformations are about the preservation of the symplectic 2-form, nothing else. The Noether theorems give conservation laws for symmetries of the system. – Lutz Lehmann Feb 06 '16 at 15:59
  • isn't vv just "second order symplectic" though, since it doesn't actually compute a(1/2), but blends it from a0 and a1? meaning it'll fully conserve phase-space area only if a(p) is of second order, but that's not the case with a gravitational n-body problem – Anton Knyazyev Feb 06 '16 at 21:03
  • otoh, i might indeed have misunderstood that. perhaps "2nd order symplectic" just means 2nd order *and* symplectic, regardless of a(pos) dependency. – Anton Knyazyev Feb 06 '16 at 22:27
  • Yes, it is symplectic for any Hamiltionian/conservative system. The gravitational n-body problem falls into this category with H(q,p)=1/2·p^T·M^(-1)·p + V(q) where M is the diagonal mass matrix and V(q) is the sum over all gravitational potential terms. -- 2nd order refers to the error in the trajectory and combined to the error in the first integrals of the system. – Lutz Lehmann Feb 06 '16 at 23:41