# Module 6

## Numerical Integration of Orbits

In the previous lectures we saw that the orbits of particles under the influence of central forces have certain simple properties. The orbit is restricted to a plane and we can write down simple equations for the conservation of energy and angular momentum. Orbits are either bound or unbound, depending on whether they have a kinetic energy greater than the potential energy. Some orbits are closed and always return to their starting point, while others are open.

We looked at two examples of central forces, the central harmonic oscillator and the Keplerian force field.  For both of these cases the orbits can be found exactly (analytically) as a function of time, which is to say that we can solve for the orbit. This is not so in general however. Usually, for any realistic potential field, the orbit is complicated and does not admit of simple or general solutions, even if the potential does not change with time.

The most famous case of a simple system leading to very complicated orbits is called the three-body-problem. Newton solved the two-body-problem (i.e. two point masses moving in their mutual gravitational field), but was unable to solve the three-body-problem. Many people have worked on this problem since Newton, and it is in fact still a very active field even today, three hundred years later! The arrangement of three bodies turns out to be very sensitive to the initial positions and/or velocities of the three bodies and displays what is now known as chaotic behaviour.

In this lecture we will examine the numerical techniques used to solve for orbits when the analytic solution is not known.  This is the way that orbits for spacecraft (such as those that went to the Moon) and stars in the potential field of the Galaxy must be calculated.

## Equations of Motion

Our aim is to determine the orbit of a star in three dimensional space in a potential field Φ(x,y,z). We know the position and velocity of the star at a given time t = t0.

The potential is time independent i.e. ∂Φ/∂t = 0.

Let the force (per unit mass) acting on the star be given by F(x,y,z) at the position r(x,y,z) so that

F = −dΦ(r)/dr

In 3-D Cartesian coordinates we can decompose the force into three components F = F(Fx,Fy,Fz), where

Fx = −∂Φ/∂x

Fy = −∂Φ/∂y

Fz = −∂Φ/∂z.

and note that these are partial derivatives.

The equations of motion of the star are (force = mass x acceleration)

d2x/dt2 = −∂Φ/∂x,

with similar forms in y and z.

Now we will consider a star at the position P(x0,y0,z0) at time t = 0 with velocities (dx0/dt, dy0/dt, dz0/dt) which we denote by (u0,v0,w0). The true (unknown) path of the star is shown in green. At P we can use the velocities and position to project the stellar path forward to Q after a time Δt. The actual position of the star after this time is R, and there will always be a small difference (error) between R and Q. In many situations, as we shall see, by choosing Δt small enough one can make this error arbitrarily small, and compute an orbit to any desired accuracy.

To a first approximation, we can compute that the star will move to a new position Q(x1,y1,z1) after a short time Δt based on the velocity of the star at P:

x1 = x0 + u0 Δt

y1 = y0 + v0 Δt

z1 = z0 + w0 Δt.

Furthermore, we can compute the change in the velocities during the interval Δt based on the acceleration due to the force acting on the star at P

u1 = u0 + FxΔt

v1 = v0 + FyΔt

w1 = w0 + FzΔt.

We now have an estimate for the position and velocity of the star at Q. One can repeat (iterate) this procedure to predict the position and velocity of the star as a function of time. Mathematically, this is termed an initial value problem. We solve the equations of motion for a single case by integrating the equations numerically rather than finding the general analytical solution.

The method described above is called Euler's method. A large range of physical problems can be solved using methods based on this technique, and as a consequence a lot of effort has gone into developing fast and accurate versions of the technique; basically, higher order approximations than Euler's linear method are used, and they are known as Runge-Kutta methods. They are much more accurate than the Euler method, which is used here merely to demonstrate the idea. We'll look at those after applying Euler's method to a simple case, simple harmonic motion in 2-D.

# Numerical solutions to stellar orbits

## Simple Harmonic Oscillator

The simple 2-D harmonic oscillator potential is

Φ(r) = &Omega2r2/2 + const,

and the equation of motion in x is (where r = r(x,y) = sqrt(x2 + y2))

d2x/dt2 = Fx = −dΦ/dr . dr/dx = −&Omega2r . x/r = −&Omega2x

with a similiar equation in y.

We insert this form into the equations above and solve the equations numerically. The equations in the X-coordinate for the nth step in the integration are:

xn+1 = xn + unΔt

un+1 = un − Ω2xnΔt.

A similar form holds for the y-coordinate.

These equations give us a scheme for "numerically integrating" the orbit. Two examples of the numerical procedure follow, firstly using GNUPLOT and then using QBASIC.

## Numerical solutions using GNUPLOT

The left column in the table below shows GNUPLOT input files which plot the motion of a simple harmonic oscillator in the (X,Y) plane. A simple Euler method is used to numerically integrate the equations of motion, while the exact solution is also plotted for comparison (you could verify that the exact solution is the right one from the given initial condidtions of the orbit, which are that at t=0 the star starts at (X,Y) = (5.0,0.0) with velocities (U,V) = (0.0,50.0).

The right hand column column plots out the difference between the exact solution and the numerical solution as a funtion of time, t. You can experiment with using a different number of steps (change imax in the program). This has the effect of making the step size larger or smaller (the step size is denoted by the variable dt. The larger the number of steps the longer the solution takes to compute and the more accurately it matches the analytical solution.

Here are the files:

shmorbit.gnp, integrate-shmorbit.gnp, shmorbit.gnp integrate-shmorbit.gnp,

 shmorbit.gnp # reset all paramters reset clear # starting position in x, u (u=dx/dt) at time t = 0.0 # starting position in y, v (v=dy/dt) t = 0.0 x = 5.0  u = 0.0 y = 0.0 v = 50.0 # step (loop) counter i = 0 # number of points along the orbit imax = 200 # time step dt = 2.0/imax # total time of integration tmax = imax*dt # SHM force law omega = 5.0 fx(x) = -x*omega**2  fy(y) = -y*omega**2 set xrange [-15:15] set yrange [-15:15] print "the real orbit is plotted with o" print "and the estimated orbit with +" # start running the integration load 'integrate-shmorbit.gnp' shmtime.gnp # reset all paramters reset clear # starting position in x, u (u=dx/dt) at time t = 0.0 # starting position in y, v (v=dy/dt) t = 0.0 x = 5.0  u = 0.0 y = 0.0 v = 50.0 # step (loop) counter i = 0 # number of points along the orbit imax = 200 # time step dt = 2.0/imax # total time of integration tmax = imax*dt # SHM force law omega = 5.0 fx(x) = -x*omega**2  fy(y) = -y*omega**2 set xrange [0:2.1] set yrange [-1.5:1.5] print "the difference between the real orbit and the " print "estimated orbit is plotted with x" # start running the integration load 'integrate-shmtime.gnp' integrate-shmorbit.gnp t = t + dt x = x + dt*u u = u + dt*fx(x) y = y + dt*v v = v + dt*fy(y)  epsilon_x = -pi/2 epsilon_y = -pi/2 Ax = 5.0 Ay = 10.0 xt = Ax*sin(omega*t-epsilon_x) yt = -Ay*cos(omega*t-epsilon_y) # plot the point with "o" at the position after the "at" set label 'o' at xt,yt center set label '+' at x,y center # dummy statement (do not remove) to fool gnuplot  # into plotting a single point plot -5000 # this is the actual loop  if (t < tmax) reread integrate-shmtime.gnp t = t + dt x = x + dt*u u = u + dt*fx(x) y = y + dt*v v = v + dt*fy(y)  epsilon_x = -pi/2 epsilon_y = -pi/2 Ax = 5.0 Ay = 10.0 xt = Ax*sin(omega*t-epsilon_x) yt = -Ay*cos(omega*t-epsilon_y) # plot the difference between the true and estimated # radius of the star as a function of time r_est  = sqrt(x**2+y**2) r_true = sqrt(xt**2+yt**2) diff = r_est - r_true set label 'x' at t,diff # dummy statement (do not remove) to fool gnuplot  # into plotting a single point plot -5000 # this is the actual loop  if (t < tmax) reread

## Numerical solutions QBASIC

Below is a version of the above program in the programming language BASIC. You should be able to run this on any PC using QBASIC, which is included on all systems. Use the command QBASIC under DOS, load the program using the "new" command in the "file" menu, and then RUN it. The time, and the X and Y positions of the star are saved to a file (called "orbit.out"). You can adjust the stepsize by increasing or decreasing the number of total steps through the variable imax. The results in the output file can be plotted using GNUPLOT. As an exercise you could modify the program to also print out the analytical orbital solution and plot both the numerical and analytical solutions using GNUPLOT.

Note that using QBASIC is much faster than using the GNUPLOT programs above, but GNUPLOT lets you watch the orbit in real time. Much faster still is to use a compiled language, such as fortran, C, C++, VisualBasic etc. etc. Such languages will compute orbits of this type hundreds of times faster than the interpreted languages which are used here.

 REM program to calculate the orbit of a star in  REM a harmonic potential x = 5! u = 0! y = 0! t = 0! v = 50! i = 0 imax = 200 dt = 2! / imax omega = 5! OPEN "orbit.out" FOR OUTPUT AS #1 PRINT #1, t, x, y FOR i = 1 TO imax t = t + dt x = x + u * dt y = y + v * dt u = u - x * omega * omega * dt v = v - y * omega * omega * dt PRINT #1, t, x, y NEXT i END Results of a calculation of a simple harmonic orbit with an exact analytical solution are shown below. The star starts at x=5.0, y=0.0 with velocities u=0.0, v=50.0. The force law is F(r) = −Ω2 r with Ω = 5.0. The exact solution is shown by the green line and has the form

x(t) = 5 sin(Ωt+π/2)

y(t) = −10 cos(Ωt+π/2)

The lower left panel shows an iterative solution with 10 steps per time unit (Δt = 0.1), clearly a very poor solution. As N increases to N=20 and N=200 the solutions get closer and closer to the exact solution. The upper right panel shows a zoomed region (indicated by the red box in the upper-left panel), showing that the solution with N=200 is still slightly different to the exact solution.

## Runge-Kutta numerical solution

Consider an equation of motion of the form

x'' = f(t,x,x'),

and that we want to solve for x(t) from initial conditions :

x(t0) = x0     x'(t0) = x'0.

Here, x'' is the 2nd derivative with respect to time t, x' is the first derivative, and the initial position is x0 and initial velocity is x'0 at time t = t0.

Let us expand the position coordinate x as a Taylor series in time, t, for a small step in time, dt.

x(t+dt) = x(t) + dt x'(t) + dt2 x''(t)/2 + dt3 x'''(t)/6 + ...

If we retain second order terms only, we get the following approximations:

x(t+dt) = x(t) + dt x'(t) + dt2 x''(t)/2
x'(t+dt) = x'(t) + dt x''(t).

Comparison with the Euler method above will show that we have improved it by including all the second order terms. There is actually only one change to be made to the Euler method - the equation which computes the change in the position includes new term, while the equation for the change in the velocity is the same as in the Euler case.

 Exercise #1: Modify the program to take into account second order terms (using the Runge-Kutta method). How much more accurate is the orbit for the simple harmonic case?

## 2) Orbits in a Kepler Potential

The Keplerian potential is

&Phi(r) = −GM/r

and the force acting at r is

F(r) = −GM/r2

in 3-D cartesian coordinates the force has components Fx, Fy, Fz where

Fx = dΦ/dx . dx/dr = −2GMx/r3

with similiar expressions in the other two components.

Inserting these into Eqns (1) above we derive the numerical form to be solved.

Here is a BASIC program to solve to orbit.

 REM program kepler.  REM This program integrates the equations  REM of motion of a particle REM in a Keplerian force field. REM Uses the not very accurate Euler method REM so that the number of points along the orbit REM needs to be large. REM initial values of variables       t = 0!       x = 5!       u = 0!       y = 0!       v = 20! REM step (loop) counter       i = 0 REM number of points along the orbit       imax = 10000 REM time step       dt = 2! / imax REM total time of integration       tmax = imax * dt REM keplerian force law       c = 3000!       OPEN "c:/gnuplot3.7/kepler.dat" FOR OUTPUT AS #1 REM   OPEN "kepler.dat" FOR OUTPUT AS #1         r = SQR(x * x + y * y)       PRINT #1, t, x, y, r       FOR i = 1 TO imax          r = SQR(x * x + y * y)          t = t + dt          x = x + dt * u          u = u - dt * x * c / r ^ 3          y = y + dt * v          v = v - dt * y * c / r ^ 3          PRINT #1, t, x, y, r       NEXT i       END

 Exercise #2: This program is quite interesting to experiment with. What happens with small changes to the program may surprise you. Start by changing the initial position or velocity of the star; after that, experiment with different step sizes. How sensitive is the procedure to the adopted step size? What indications are there that an orbit may have been computed completely incorrectly? Recalling that certain quantities are conserved as a particle moves along its orbit, suggest a method to help decide whether an orbital calculation is reliable or not.

## Runge-Kutta-Nystrom numerical solution

A yet more accurate integrator is the RKN, or Runge-Kutta-Nystrom method. It extends the (2nd order) RK method above, and is a fourth ourder method, taking into account all terms to dt4. From the step n to step n+1, the following quantities are evaluated:

An = dt f(t,xn,x'n)/2

βn = dt (x'n + An/2)/2

Bn = dt f(t+dt/2,xnn,x'n+An)/2

Cn = dt f(t+dt/2,xnn,x'n+Bn)/2

δn = dt(x'n + Cn)

Dn = dt f(t+dt,xnn,x'n+2Cn)/2.

This is more general than we need, since the force on the particle depends only on its position, so these reduce to :

An = dt f(xn)/2

βn = dt (x'n + An/2)/2

Bn = dt f(xnn)/2

Cn = dt f(xnn)/2

δn = dt(x'n + Cn)

Dn = dt f(xnn)/2.

The steps involved in getting from the from nth step to the (n+1)th step in the integration are:

xn+1 = xn + dt (vn + Kn)

where Kn = (An + Bn + Cn)/3

vn+1 = vn + Jn

where Jn = (An + 2Bn + 2Cn + Dn)/3.

 Exercise #3: This is a more challenging exercise. Write a program which solves the Kepler orbit using the Runge-Kutta-Nystrom method. You can modify the QBASIC program, or write your own in some other language if you want. How much more accurate/stable is this method compared to the Euler and Runge-Kutta integrators?