N body ODE simulation

Question

I have been staring at the below matlab code, which simulates planetary motion, for a while now and I could really use a fresh set of eyes. I cannot understand why it is giving straight line trajectories for the planets' motion. All the initial data has been taken from Geometric Numerical Integration by Hairer Pg. 13. Thanks for your help.

For the Hamiltonian $$H(q,p) = \frac{1}{2}\sum_{n=1}^N \frac{\|p_n\|^2}{m_n} -G \sum_{n=1}^{N-1} \sum_{j=n+1}^N \frac{m_n m_j}{\|q_n -q_j\|}$$ equation of motion of the planets are $$\frac{dq_n}{dt} = \frac{p_n}{m_n},\\ \frac{dp_n}{dt} = -G\sum_{j \neq n} \frac{m_n m_j}{\|q_n -q_j\|^3}(q_n-q_j)$$ which gives the following Stormer-Verlet method ($i$ is the time index) $$p_n^{i+1/2} = p_n^i -\frac{\Delta t}{2}G\sum_{j \neq n} \frac{m_n m_j}{\|q_n^i -q_j^i\|^3}(q_n^i-q_j^i),\\ q_n^{i+1} = q_n^i + \Delta t \frac{p_n^{i+1/2}}{m_n},\\ p_n^{i+1} = p_n^{i+1/2} -\frac{\Delta t}{2}G\sum_{j \neq n} \frac{m_n m_j}{\|q_n^{i+1} -q_j^{i+1}\|^3}(q_n^{i+1}-q_j^{i+1}).$$

    function NbodyODE(dt, T, N, q0, p0, m, G)
    % This function solves NbodyODE with symplectic Stormer-Verlet method

if nargin < 1, dt = 10; end                 % one time step-size is equal to 10 earth days
if nargin < 2, T = 200000; end              % simulate outer solar system for 200 000 earth days
if nargin < 3, N = 6; end
if nargin < 4, q0 = [0 0 0;                 % initial positions of the planets, the Sun is at (0,0,0)
                    -3.5023653 -3.8169847 -1.5507963;
                     9.0755314 -3.0458353 -1.6483708;
                     8.3101420 -16.2901086 -7.2521278;
                     11.4707666 -25.7294829 -10.8169456;
                     -15.5387357 -25.2225594 -3.1902382
                    ]; end
if nargin < 5, p0 = [0 0 0;                 % initial momenta of the planets
                    0.00565429 -0.00412490 -0.00190589;
                     0.00168318 0.00483525 0.00192462;
                     0.00354178 0.00137102 0.00055029;
                     0.00288930 0.00114527 0.00039677;
                     0.00276725 -0.00170702 -0.00136504
                    ]; end
if nargin < 6, m = [1.00000597682;          % masses of the solar bodies.
                    0.000954786104043;
                     0.000285583733151;
                     0.0000437273164546;
                     0.0000517759138449;
                     1/1.3e8
                    ]; end
if nargin < 7, G = 2.95912208286e-4; end    % Gravitational constant

if (size(q0,1) ~= N || size(p0,1) ~= N || size(m,1) ~= N)
    error('input dimensions mismatch, correct input data')
end

t = 0:dt:T;
L = length(t);
q = zeros(N,3,L); p = q; q(:,:,1) = q0; p(:,:,1) = p0;

[q1, p1] = SV(m, dt, L, N, q, p, G);

for n = 1:N
    x = reshape(q1(n,1,:),L,1); y = reshape(q1(n,2,:),L,1); z = reshape(q1(n,3,:),L,1);
    plot3(x,y,z), hold on
end

function [q, p]= SV(m, dt, L, N, q, p, G)
    Vf = @(m1,m2,x,y) (G*m1*m2*(x-y)/(norm(x-y))^3);

tau = zeros(N,3,L);
for n = 1:N
    for j = [1:n-1 n+1:N]
        tau(n,:,1) = tau(n,:,1) + Vf(m(n), m(j), q(n,:,1), q(j,:,1));
    end
end

for i = 1:L-1
    for n = 1:N
        p(n,:,i+1) = p(n,:,i) - dt/2*tau(n,:,i);
        q(n,:,i+1) = q(n,:,i) + dt/m(n)*p(n,:,i+1);
    end

    for n = 1:N
        for j = [1:n-1 n+1:N]
            tau(n,:,i+1) = tau(n,:,i+1) + Vf(m(n), m(j), q(n,:,i+1), q(j,:,i+1));
        end
    end

    for n = 1:N
        p(n,:,i+1) = p(n,:,i+1) - dt/2*tau(n,:,i+1);
    end
end

Answer 1

On the proper way to implement the Verlet method

You need to first compute all interaction forces, and only after all forces have been computed can you update the position vector. Not doing this and mixing interaction with position updates usually results in an unphysical drift.

You should store the forces or acceleration in an extra array, at the moment you are computing the forces for each computation twice.

p(i+1) = p(i)+f(i) * dt/2
q(i+1) = q(i)+p(i)/m * dt
f(i+1) = forces(q(i+1))
p(i+1) = p(i+1)+f(i+1) * dt/2

One does not need to keep a history of the force terms, one can always use the same array (if one were implementing the method with fine-tuned memory management).

Improved Verlet method - Verlet4

Also try out the "improved" Verlet method where you chain

Verlet(b0*dt); Verlet(-b1*dt); Verlet(b0*dt)

for one time step dt with b0=1/(2-2^(1/3))=1.35120719196 and b1=2*b0-1=1.70241438392. It should be in the book under "composition methods" or "methods obtained by composition". It gives you a 4th order method for little extra programming effort.

---

How to initialize and avoid drift of the system

Study the implementation in http://www.boost.org/doc/libs/1_54_0/libs/numeric/odeint/examples/solar_system.cpp, especially the initialization. The total impulse of the system is not zero. The second set of vectors is not the initial impulse but the initial velocity, so they need to be converted to impulses. Do this via

pmean = zeros(1,3); qmean = pmean; tmass = 0;
for n=1:N
    pmean = pmean + p0(n,:)*mass(n);
    qmean = qmean + q0(n,:)*mass(n);
    tmass = tmass + mass(n);
end
qmean = qmean / tmass; pmean = pmean / tmass
for n=1:N
    q0(n,:) = q0(n,:) - qmean;
    p0(n,:) = (p0(n,:) - pmean)*mass(n);
end