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 threebodyproblem. Newton solved the twobodyproblem (i.e. two point masses moving in their mutual gravitational field), but was unable to solve the threebodyproblem. 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.
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 3D Cartesian coordinates we can decompose the force into three components F = F(F_{x},F_{y},F_{z}), where
F_{x} = −∂Φ/∂x
F_{y} = −∂Φ/∂y
F_{z} = −∂Φ/∂z.
and note that these are partial derivatives.
The equations of motion of the star are (force = mass x acceleration)
d^{2}x/dt^{2} = −∂Φ/∂x,
with similar forms in y and z.
Now we will consider a star at the position
P(x_{0},y_{0},z_{0}) at time t = 0 with velocities
(dx_{0}/dt, dy_{0}/dt, dz_{0}/dt) which we denote by
(u_{0},v_{0},w_{0}).
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(x_{1},y_{1},z_{1}) after a short time Δt based on the velocity of the star at P:
x_{1} = x_{0} + u_{0} Δt
y_{1} = y_{0} + v_{0} Δt
z_{1} = z_{0} + w_{0} Δ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
u_{1} = u_{0} + F_{x}Δt
v_{1} = v_{0} + F_{y}Δt
w_{1} = w_{0} + F_{z}Δ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 RungeKutta 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 2D.
The simple 2D harmonic oscillator potential is
Φ(r) = &Omega^{2}r^{2}/2 + const,
and the equation of motion in x is (where r = r(x,y) = sqrt(x^{2} + y^{2}))
d^{2}x/dt^{2} = F_{x} = −dΦ/dr . dr/dx = −&Omega^{2}r . x/r = −&Omega^{2}x
with a similiar equation in y.
We insert this form into the equations above and solve the equations numerically. The equations in the Xcoordinate for the nth step in the integration are:
x_{n+1} = x_{n} + u_{n}Δt
u_{n+1} = u_{n} − Ω^{2}x_{n}Δt.
A similar form holds for the ycoordinate.
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.
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, integrateshmorbit.gnp, shmorbit.gnp integrateshmorbit.gnp,shmorbit.gnp
# reset all paramters
# starting position in x, u (u=dx/dt) at time t =
0.0
# step (loop) counter
# time step
# total time of integration
# SHM force law
set xrange [15:15]
print "the real orbit is plotted with o"
# start running the integration

shmtime.gnp
# reset all paramters
# starting position in x, u (u=dx/dt) at time t =
0.0
# step (loop) counter
# time step
# total time of integration
# SHM force law
set xrange [0:2.1]
print "the difference between the real orbit and the
"
# start running the integration

integrateshmorbit.gnp
t = t + dt
epsilon_x = pi/2
xt = Ax*sin(omega*tepsilon_x)
# plot the point with "o" at the position after the
"at"
# dummy statement (do not remove) to fool gnuplot
# this is the actual loop

integrateshmtime.gnp
t = t + dt
epsilon_x = pi/2
xt = Ax*sin(omega*tepsilon_x)
# plot the difference between the true and estimated
r_est = sqrt(x**2+y**2)
# dummy statement (do not remove) to fool gnuplot
# this is the actual loop

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!
i = 0
omega = 5! OPEN "orbit.out" FOR OUTPUT AS #1 PRINT #1, t, x, y FOR i = 1 TO imax t = t + 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 upperleft panel), showing that the solution with N=200 is still slightly different to the exact 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(t_{0}) = x_{0} x'(t_{0}) = x'_{0}.
Here, x'' is the 2nd derivative with respect to time t, x' is the first derivative, and the initial position is x_{0} and initial velocity is x'_{0} at time t = t_{0}.
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) + dt^{2} x''(t)/2 + dt^{3} x'''(t)/6 + ...
If we retain second order terms only, we get the following approximations:
x(t+dt) = x(t) + dt x'(t) + dt^{2} 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 RungeKutta method). How much more accurate is the orbit for the simple harmonic case? 
&Phi(r) = −GM/r
and the force acting at r is
F(r) = −GM/r^{2}
in 3D cartesian coordinates the force has components F_{x}, F_{y}, F_{z} where
F_{x} = dΦ/dx . dx/dr = −2GMx/r^{3}
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 initial values of variables t = 0!
REM step (loop) counter
REM number of points along the orbit
REM time step
REM total time of integration
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)
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. 
A yet more accurate integrator is the RKN, or RungeKuttaNystrom method. It extends the (2nd order) RK method above, and is a fourth ourder method, taking into account all terms to dt^{4}. From the step n to step n+1, the following quantities are evaluated:
A_{n} = dt f(t,x_{n},x'_{n})/2
β_{n} = dt (x'_{n} + A_{n}/2)/2
B_{n} = dt f(t+dt/2,x_{n}+β_{n},x'_{n}+A_{n})/2
C_{n} = dt f(t+dt/2,x_{n}+β_{n},x'_{n}+B_{n})/2
δ_{n} = dt(x'_{n} + C_{n})
D_{n} = dt f(t+dt,x_{n}+δ_{n},x'_{n}+2C_{n})/2.
This is more general than we need, since the force on the particle depends only on its position, so these reduce to :
A_{n} = dt f(x_{n})/2
β_{n} = dt (x'_{n} + A_{n}/2)/2
B_{n} = dt f(x_{n}+β_{n})/2
C_{n} = dt f(x_{n}+β_{n})/2
δ_{n} = dt(x'_{n} + C_{n})
D_{n} = dt f(x_{n}+δ_{n})/2.
The steps involved in getting from the from nth step to the (n+1)th step in the integration are:
x_{n+1} = x_{n} + dt (v_{n} + K_{n})
where K_{n} = (A_{n} + B_{n} + C_{n})/3
v_{n+1} = v_{n} + J_{n}
where J_{n} = (A_{n} + 2B_{n} + 2C_{n} + D_{n})/3.
Exercise #3: This is a more challenging exercise. Write a program which solves the Kepler orbit using the RungeKuttaNystrom 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 RungeKutta integrators? 