410-505 2.1 Solar System Orbits and Runge-Kutta Algorithms Fall 2011

Dynamics of Solar System Objects

The orbits around the Sun of solar system bodies\cite{orbital-elements} are determined by Newton's Law of Universal Gravitation, and Newton's equations of motion. These are systems of ordinary differential equations with solutions specified by initial conditions at some time. The orbital elements of objects involved provide precisely such initial conditions.

Motion of planets around the Sun

Johannes Kepler discovered three laws relating to the motion of planets around the Sun.

1. Planets move in elliptical orbits, with the Sun at one focus.
2. The line joining the planet to the Sun sweeps out equal areas in equal times.
3. The square of the planet's orbital period is proportional to the cube of the semi-major axis of the orbit.

In the figure, $a$ is the semimajor axis, $b$ is the semiminor axis. The Perihelion and Aphelion points lie on the major axis, and are the closet and farthest distances, respectively.

Newton's law of gravity

Isaac Newton proved that Kepler's laws are obeyed by an isolated bound system of any two objects interacting through an inverse square central force. Gravitation is such a force:

$${\bf F}_{12}=-\frac{Gm_1m_2}{ r_{12}^3}{\bf r}_{12}\;,\qquad {\bf F}_{21}=-\frac{Gm_1m_2}{ r_{21}^3}{\bf r}_{21}\;.$$

These forces obey Newton's Third Law: $${\bf F}_{12} = - {\bf F}_{21}\;,$$ which is illustrated in figure above.

Exact solution of the Kepler problem

Since the Sun is by far the most massive object in the Solar System, most discussions use a Copernican coordinate system with the Sun fixed at the origin. However, this system is only approximately inertial. The approximation is an excellent one for comets, and even for a massive planet like Jupiter. If the two bodies have comparable masses, e.g., in a binary star system, then the center of mass system is inertial. The exact solution in this case can be expressed in terms of the reduced mass and relative cordinate $$\mu={m_1m_2\over m_1+m_2}\;, \quad{\bf r}={\bf r}_1-{\bf r}_2\;, \quad m_1{\bf r}_1+m_2{\bf r}_2 = 0\;.$$ If $m_2=M>>m_1=m$ then $\mu\simeq m$, ${\bf r}\simeq {\bf r}_1$, and ${\bf r}_2\simeq0$, which is the Copernican coordinate system.

The Kepler orbit can be written in parametric form $$r(\theta)={a(1-\epsilon^2)\over1-\epsilon\cos\theta}\;,\qquad b=a\sqrt{1-\epsilon^2}\;,$$ where $a$ is the semimajor axis, $\epsilon$ is the {\it eccentricity}, and $b$ the semiminor axis. The orbital speed is given by $$v=\sqrt{G(m_1+m_2)\left({2\over r}-{1\over a}\right)}\;.$$

The conserved energy of the system is given by $$E={1\over2}\mu v^2-{Gm_1m_2\over r}=-{Gm_1m_2\over2a}\;.$$

The period of the orbit is given by Kepler's third law $$T^2={4\pi^2\over G(m_1+m_2)}a^3\;.$$

The orbital angular momentum can be computed using ${\bf L}={\bf r}\times(\mu{\bf v})$.

Of course, even the c.m.\ frame is not strictly inertial if one takes into account the motion of the solar system through the galaxy.

Solving Newton's Equation of Motion Numerically

We need to solve Newton's equation of motion $$\frac{d^2{\bf r}}{dt^2} = - \frac{GM}{r^3} {\bf r}$$ given the initial position and velocity of the planet. This is a second order ordinary differential equation.

For computational purposes, it is convenient to rewrite higher order differential equations as systems of coupled first order equations.

To do this, use the velocity ${\bf v}$ to obtain \begin{eqnarray} \frac{d{\bf r}}{dt} &=& {\bf v} \;, \\ \frac{d{\bf v}}{dt}&=& - \frac{GM}{r^3} {\bf r} \;. \end{eqnarray}

Algorithms for solving such systems are discussed in Chapter 16 of Numerical Recipes\cite{recipes-c16-1}.

Euler's algorithm

A very simple algorithm was introduced by Euler to solve sets of coupled first order differential equations.

Consider the simple equation $$\frac{dx}{dt} = f(x,t) \;.$$

Given an initial value $x(0)$, the differential equation determines an unique solution $x(t)$.

Suppose that this solution is a smooth function with well behaved derivatives. Given $x(t)$ at any time $t$, the value of $x(x+dt)$ at a later time can be approximated by a Taylor expansion $$x(t+dt) = x(t) + \left.\frac{dx}{dt}\right|_{t}dt + {\cal O}(dt^2) = x(t) + f(x,t)dt + {\cal O}(dt^2) \;.$$

Given the initial value $x(0)$ and a step size $dt$, the solution is determined as a sequence of points $$x(0), x(dt) = x(0)+f(x(0),0)dt, x(2dt) = x(dt) + f(x(dt),dt)dt, \ldots$$

The advantage of writing the equations in first order form is now apparent. A general system of ordinary differential equations can be written in first order form by replacing $x$ and $f$ with vectors $$\frac{d{\bf x}}{dt} = {\bf f}({\bf x},t) \;.$$ and Euler's algorithm can be expressed in vector form $${\bf x}(0), {\bf x}(dt) = {\bf x}(0)+{\bf f}({\bf x}(0),0)dt, {\bf x}(2dt) = {\bf x}(dt) + {\bf f}({\bf x}(dt),dt)dt, \ldots$$

Astronomical units

It is a good idea to avoid using very large or very small numerical values in computer programs, in order to avoid potential problems with numerical overflow or underflow. This can usually be done by using a system of units that is natural to the problem being considered.

• Earth's semimajor axis $1$ AU $=1.496\times10^{11}$ m
• Kepler's third law $$G(M_{\rm Sun}+m_{\rm Earth}) = {4\pi^2 a^3\over T^2}$$
• Measure distance in AU and time in years $$G(M_{\rm Sun}+m_{\rm Earth}) = 4\pi^2 \; {\rm AU}^3/{\rm yr}^2$$

Setting Initial Values

The orbital elements of solar system bodies provide initital data that in principle determine their orbits. These include the semimajor axis length $a$ and the eccentricity $\epsilon$. From these two data the initial position and velocity, at the aphelion point for example, can be fixed so as to reproduce the desired Kepler orbit.

• Choose semimajor axis in the $x$ direction and specify the aphelion distance $x(t=0)$ $$x(0) = a(1+\epsilon) \;.$$
• Choose the velocity at aphelion in the positive $y$ direction and set $v_y(t=0)$ equal to its magnitude.
• The speed is determined from the orbit equations: $$v=\sqrt{G(m_1+m_2)\left({2\over r}-{1\over a}\right)}\;.$$
• Solve this equation at $t=0$ for the eccentricity and semimajor axis $$\epsilon = x(0)v_y(0)^2/(G(M+m)) - 1 \;,$$ $$a = x(0)/(1+\epsilon) \;.$$

Determining the period precisely

The period of the orbit can be determined more accurately by interpolating the numerical solution to obtain the time at aphelion.

Consider a set of data points $(x_i,y_i)$ in 2 dimensional space.

Any two points $(x_0,y_0)$, $(x_1,y_1)$ can be connected by a straight line $$y(x)=\frac{(x-x_1)}{(x_0-x_1)}y_0 + \frac{(x-x_0)}{(x_1-x_0)}y_1 \;.$$

A parabola can be drawn to pass through any three points: $$y(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}y_0 + \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}y_1 + \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}y_2 \;.$$

Algorithms for interpolation and extrapolation are discussed in Chapter 3 of Numerical Recipes\cite{recipes-c3-1}.

The Euler-Cromer Algorithm

An alternative to the simple Euler algorithm was introduced by Cromer\cite{cromer-1981}, who found that it tended to conserve energy better in problems with periodic motion.

A. Cromer, "Stable Solutions using the Euler Approximation", Am. J. Phys. 49, 455 (1981), http://dx.doi.org/10.1119/1.12478.

The modification consisted in using an updated value of the velocity of the particle to compute the position: \begin{eqnarray*} {\bf v}(t+dt) &=& {\bf v}(t) + {\bf a({\bf r}(t))}dt \;, \\ {\bf r}(t+dt) &=& {\bf r}(t) + {\bf v}(t+dt)dt \;. \end{eqnarray*} Notice the change from the simple Euler algorithm in the second equation!

Runge-Kutta Methods for ODE Systems

An ordinary differential equation (ODE) has {\it one} independent variable $x$, and one or more dependent variables $y(x),\ldots$.

Physicists are usually interested in a particular solution of a differential equation given a particular set of initial values of the variables.

Mathematicians prefer to consider the family of solutions of the differential equation as shown in the following figure.

Example: $${{\rm d}y\over {\rm d}x}=-y \;,$$ is a linear first order ODE. It is linear in the solution $y$ and first order in derivative ${\rm d}/{\rm d}x$. It describes a {\it family of solutions} or {\it curves} in the $x{-}y$ plane. There is a unique curve passing through each point in the plane. The slope of the curve at any point is given by the differential equation.

To obtain a unique solution of the ODE a {\it condition} must be specified to fix the arbitrary {\it integration constant} which parametrizes the different solutions. For example, an {\it initial condition} on $y(x=0)$ can be chosen.

Second order ODE's and vector notation

Example: The nonlinear damped and driven pendulum equation $${{\rm d}^2\theta\over {\rm d}t^2}+q{{\rm d}\theta\over {\rm d}t} +\sin\theta=b\cos(\omega_0t)\;,$$ is a second order nonlinear ODE with $t$ being the independent variable, which depends on constant parameters $q$ (damping coefficient), $b$ (driving amplitude), and $\omega_0$ (driving angular frequency).

It is always possible to write an ODE of any order as a system of first order ODE's. In this example \begin{eqnarray} {{\rm d}\theta\over{\rm d}t} &=& \omega \;,\\ {{\rm d}\omega\over{\rm d}t} &=& -q\omega-\sin\theta+b\cos(\omega_0t)\;, \\ \end{eqnarray} where $\omega(t)$ is the angular velocity. This system of equations defines a family of curves in $t{-}\theta{-}\omega$ space. To select a particular solution, one can impose initial conditions $\theta(t=0)$, $\omega(t=0)$.

The pedulum equations can be written in compact vector notation $${{\rm d}{\bf y}\over{\rm d}x} \equiv {{\rm d}\over{\rm d}x} \begin{pmatrix}y_1\\ y_2\\\end{pmatrix} \equiv {{\rm d}\over{\rm d}t}\begin{pmatrix}\theta\\\omega\\\end{pmatrix}= \begin{pmatrix}\omega\\ -q\omega-\sin\theta+b\cos(\omega_0t)\\\end{pmatrix} \equiv \begin{pmatrix}f_1(y_1,y_2,x)\\ f_2(y_1,y_2,x)\\\end{pmatrix} \equiv {\bf f}({\bf y}(x),x)\;,$$ where ${\bf y}$ is called the {\it solution vector}, and ${\bf f}$ is called the {\it flow vector} or {\it derivative vector}. To solve the vector differential equation we need to find the unknown solution vector ${\bf y}$ as a function of the independent variable $x$, give a flow vector ${\bf f}$ which is a known function of ${\bf y}(x)$ and $x$.

Euler's algorithm in vector notation

This vector notation is also very convenient for expressing many algorithms used to solve ordinary differential equations. For example, the Euler algorithm from Chapter 1 can be expressed simply $${\bf y}(x+dx)= {\bf y}(x) + dx\;{\bf f}({\bf y}(x),x)\;,$$ where $dx$ is a step size in $x$.

Runge-Kutta methods

Euler's algorithm is not very accurate and also suffers from instabilities. The mathematicians Runge and Kutta\cite{runge-kutta} discovered a class of algorithms that are based on making multiple Euler steps. These {\it Runge-Kutta algorithms} have been found to work very well in many applications. They are described in Numerical Recipes\cite{recipes-c16-1}.

Second order Runge-Kutta algorithm

The simple Euler algorithm uses the slope (tangent to the curve) of the solution $y(x)$ at $x$ to estimate $$y(x+h) \simeq y(x) + h\left.{{\rm d}y\over {\rm d}x}\right|_x =y(x) + hf(y,x) \equiv y(x) + k_1\;.$$

The second order Runge-Kutta algorithm takes a {\it half} Euler step to the point $$x+0.5h\;,\qquad\qquad y(x) + 0.5k_1\;.$$ There is a {\it different} solution curve $y^{(1)}(x)$ passing through this half-way point. The slope at this point is used to compute $$k_2\equiv h\left.{{\rm d}y^{(1)}\over {\rm d}x}\right|_{x+0.5h} =hf(y(x)+0.5k_1,x+0.5h)\;.$$ The second order Runge-Kutta estimate is then computed $$y(x+h)\simeq y(x) + k_2\;.$$ It is easy to show that the error is ${\cal O}(h^3)$. The left hand side can be expanded in a Taylor series: \begin{eqnarray} y(x+h) &=& y(x)+h{{\rm d}y\over{\rm d}x}+{h^2\over2}{{\rm d}^2y\over{\rm d}x^2} +{\cal O}(h^3)\\ &=& y(x) + hf(y,x)+{h^2\over2}{{\rm d}f(y,x)\over{\rm d}x} \\ &=& y(x) + hf(y,x) +{h^2\over2}\left[ {\partial f\over\partial y}{{\rm d}y\over{\rm d}x} + {\partial f\over\partial x}\right] +{\cal O}(h^3)\\ &=& y(x) + hf(y,x) +{h^2\over2}\left[ {\partial f\over\partial y}f + {\partial f\over\partial x} \right]+{\cal O}(h^3)\;. \\ \end{eqnarray} The right hand side can also be expanded $$y(x) + k_2 = y(x) + hf\left(y+{1\over2}k_1,x+{h\over2}\right) =y(x) + hf(y,x) +h\left[ {\partial f\over\partial y}{k_1\over 2} + {h\over2}{\partial f\over\partial x} \right] +{\cal O}(h^3)\;.$$ The two sides agree through terms of ${\cal O}(h^3)$.

Fourth order Runge-Kutta algorithm

Runge and Kutta showed that by combining the results of two additional Euler steps, the error can be reduced to ${\cal O}(h^5)$.

It is also easy to see that these algorithms can be extended to arbitrarily large systems of first order ODE's. The vector form of the fourth order Runge-Kutta algorithm is \begin{eqnarray} {\bf k}_1 &=& h\;{\bf f}\left(\vphantom{{1\over2}}{\bf y}(x),x\right)\\ {\bf k}_2 &=& h\;{\bf f}\left({\bf y}(x)+{1\over2}{\bf k}_1,x+{h\over2}\right)\\ {\bf k}_3 &=& h\;{\bf f}\left({\bf y}(x)+{1\over2}{\bf k}_2,x+{h\over2}\right)\\ {\bf k}_4 &=& h\;{\bf f}\left(\vphantom{{1\over2}}{\bf y}(x)+{\bf k}_3,x+h\right)\\ {\bf y}(x+h) &=&{\bf y}(x) + {{\bf k}_1 + 2{\bf k}_2 + 2{\bf k}_3 + {\bf k}_4\over6} + {\cal O}(h^5)\;. \\ \end{eqnarray}