Computational Physics Lectures: Ordinary differential equations

Morten Hjorth-Jensen, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

Date: 2017

Differential equations program

• Ordinary differential equations, Runge-Kutta method,chapter 8

• Ordinary differential equations with boundary conditions: one-variable equations to be solved by shooting and Green's function methods, chapter 9

• We can solve such equations by a finite difference scheme as well, turning the equation into an eigenvalue problem. Still one variable. Done in projects 1 and 2.

• If we have more than one variable, we need to solve partial differential equations, see Chapter 10

The material on differential equations is covered by chapters 8, 9 and 10. Two of the final projects deal with ordinary differential equations.

Differential Equations, chapter 8

The order of the ODE refers to the order of the derivative on the left-hand side in the equation

$$ \begin{equation} \frac{dy}{dt}=f(t,y). \label{_auto1} \tag{1} \end{equation} $$

This equation is of first order and $f$ is an arbitrary function. A second-order equation goes typically like

$$ \begin{equation} \frac{d^2y}{dt^2}=f(t,\frac{dy}{dt},y). \label{_auto2} \tag{2} \end{equation} $$

A well-known second-order equation is Newton's second law

$$ \begin{equation} m\frac{d^2x}{dt^2}=-kx, \label{eq:newton} \tag{3} \end{equation} $$

where $k$ is the force constant. ODE depend only on one variable

Differential Equations

partial differential equations like the time-dependent Schr\"odinger equation

$$ \begin{equation} i\hbar\frac{\partial \psi({\bf x},t)}{\partial t}= -\frac{\hbar^2}{2m}\left( \frac{\partial^2 \psi({\bf r},t)}{\partial x^2} + \frac{\partial^2 \psi({\bf r},t)}{\partial y^2}+ \frac{\partial^2 \psi({\bf r},t)}{\partial z^2}\right) + V({\bf x})\psi({\bf x},t), \label{_auto3} \tag{4} \end{equation} $$

may depend on several variables. In certain cases, like the above equation, the wave function can be factorized in functions of the separate variables, so that the Schroedinger equation can be rewritten in terms of sets of ordinary differential equations. These equations are discussed in chapter 10. Involve boundary conditions in addition to initial conditions.

Differential Equations

We distinguish also between linear and non-linear differential equation where for example

$$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t), \label{_auto4} \tag{5} \end{equation} $$

is an example of a linear equation, while

$$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t)-g(t)y^2(t), \label{_auto5} \tag{6} \end{equation} $$

is a non-linear ODE.

Differential Equations

Another concept which dictates the numerical method chosen for solving an ODE, is that of initial and boundary conditions. To give an example, if we study white dwarf stars or neutron stars we will need to solve two coupled first-order differential equations, one for the total mass $m$ and one for the pressure $P$ as functions of $\rho$

$$ \frac{dm}{dr}=4\pi r^{2}\rho (r)/c^2, $$

and

$$ \frac{dP}{dr}=-\frac{Gm(r)}{r^{2}}\rho (r)/c^2. $$

where $\rho$ is the mass-energy density. The initial conditions are dictated by the mass being zero at the center of the star, i.e., when $r=0$, yielding $m(r=0)=0$. The other condition is that the pressure vanishes at the surface of the star.

In the solution of the Schroedinger equation for a particle in a potential, we may need to apply boundary conditions as well, such as demanding continuity of the wave function and its derivative.

Differential Equations

In many cases it is possible to rewrite a second-order differential equation in terms of two first-order differential equations. Consider again the case of Newton's second law in Eq. (eq:newton). If we define the position $x(t)=y^{(1)}(t)$ and the velocity $v(t)=y^{(2)}(t)$ as its derivative

$$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=\frac{dx(t)}{dt}=y^{(2)}(t), \label{_auto6} \tag{7} \end{equation} $$

we can rewrite Newton's second law as two coupled first-order differential equations

$$ \begin{equation} m\frac{dy^{(2)}(t)}{dt}=-kx(t)=-ky^{(1)}(t), \label{eq:n1} \tag{8} \end{equation} $$

and

$$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=y^{(2)}(t). \label{eq:n2} \tag{9} \end{equation} $$

Differential Equations, Finite Difference

These methods fall under the general class of one-step methods. The algoritm is rather simple. Suppose we have an initial value for the function $y(t)$ given by

$$ \begin{equation} y_0=y(t=t_0). \label{_auto7} \tag{10} \end{equation} $$

We are interested in solving a differential equation in a region in space $[a,b]$. We define a step $h$ by splitting the interval in $N$ sub intervals, so that we have

$$ \begin{equation} h=\frac{b-a}{N}. \label{_auto8} \tag{11} \end{equation} $$

With this step and the derivative of $y$ we can construct the next value of the function $y$ at

$$ \begin{equation} y_1=y(t_1=t_0+h), \label{_auto9} \tag{12} \end{equation} $$

and so forth.

Differential Equations

If the function is rather well-behaved in the domain $[a,b]$, we can use a fixed step size. If not, adaptive steps may be needed. Here we concentrate on fixed-step methods only. Let us try to generalize the above procedure by writing the step $y_{i+1}$ in terms of the previous step $y_i$

$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h\Delta(t_i,y_i(t_i)) + O(h^{p+1}), \label{_auto10} \tag{13} \end{equation} $$

where $O(h^{p+1})$ represents the truncation error. To determine $\Delta$, we Taylor expand our function $y$

$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}) + O(h^{p+1}), \label{eq:taylor} \tag{14} \end{equation} $$

where we will associate the derivatives in the parenthesis with

$$ \begin{equation} \Delta(t_i,y_i(t_i))=(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}). \label{eq:delta} \tag{15} \end{equation} $$

Differential Equations

We define

$$ \begin{equation} y'(t_i)=f(t_i,y_i) \label{_auto11} \tag{16} \end{equation} $$

and if we truncate $\Delta$ at the first derivative, we have

$$ \begin{equation} y_{i+1}=y(t_i) + hf(t_i,y_i) + O(h^2), \label{eq:euler} \tag{17} \end{equation} $$

which when complemented with $t_{i+1}=t_i+h$ forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of $O(h^2)$, however the total error is the sum over all steps $N=(b-a)/h$, yielding thus a global error which goes like $NO(h^2)\approx O(h)$.

Differential Equations

To make Euler's method more precise we can obviously decrease $h$ (increase $N$). However, if we are computing the derivative $f$ numerically
by for example the two-steps formula

$$ f'_{2c}(x)= \frac{f(x+h)-f(x)}{h}+O(h), $$

we can enter into roundoff error problems when we subtract two almost equal numbers $f(x+h)-f(x)\approx 0$. Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like. As an example, consider Newton's equation rewritten in Eqs. (eq:n1) and (eq:n2). We define $y_0=y^{(1)}(t=0)$ an $v_0=y^{(2)}(t=0)$. The first steps in Newton's equations are then

$$ \begin{equation} y^{(1)}_1=y_0+hv_0+O(h^2) \label{_auto12} \tag{18} \end{equation} $$

and

$$ \begin{equation} y^{(2)}_1=v_0-hy_0k/m+O(h^2). \label{_auto13} \tag{19} \end{equation} $$

Differential Equations

The Euler method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at $y^{(1)}_1$ using the velocity at $y^{(2)}_0=v_0$. A simple variation is to determine $y^{(1)}_{n+1}$ using the velocity at $y^{(2)}_{n+1}$, that is (in a slightly more generalized form)

$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+h y^{(2)}_{n+1}+O(h^2) \label{_auto14} \tag{20} \end{equation} $$

and

$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2). \label{_auto15} \tag{21} \end{equation} $$

The acceleration $a_n$ is a function of $a_n(y^{(1)}_{n}, y^{(2)}_{n},t)$ and needs to be evaluated as well. This is the Euler-Cromer method.

Differential Equations

Let us then include the second derivative in our Taylor expansion. We have then

$$ \begin{equation} \Delta(t_i,y_i(t_i))=f(t_i)+\frac{h}{2}\frac{df(t_i,y_i)}{dt}+O(h^3). \label{_auto16} \tag{22} \end{equation} $$

The second derivative can be rewritten as

$$ \begin{equation} y''=f'=\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f \label{eq:derivatives} \tag{23} \end{equation} $$

and we can rewrite Eq.\ (eq:taylor) as

$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) +hf(t_i)+ \frac{h^2}{2}\left(\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f\right) + O(h^{3 }), \label{_auto17} \tag{24} \end{equation} $$

which has a local approximation error $O(h^{3 })$ and a global error $O(h^{2})$.

Differential Equations

These approximations can be generalized by using the derivative $f$ to arbitrary order so that we have

$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(f(t_i,y_i)+\dots f^{(p-1)}(t_i,y_i) \frac{h^{p-1}}{p!}) + O(h^{p+1}). \label{_auto18} \tag{25} \end{equation} $$

These methods, based on higher-order derivatives, are in general not used in numerical computation, since they rely on evaluating derivatives several times. Unless one has analytical expressions for these, the risk of roundoff errors is large.

Differential Equations

The most obvious improvements to Euler's and Euler-Cromer's algorithms, avoiding in addition the need for computing a second derivative, is the so-called midpoint method. We have then

$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+\frac{h}{2}\left(y^{(2)}_{n+1}+y^{(2)}_{n}\right)+O(h^2) \label{_auto19} \tag{26} \end{equation} $$

and

$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2), \label{_auto20} \tag{27} \end{equation} $$

yielding

$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n}+\frac{h^2}{2}a_n+O(h^3) \label{_auto21} \tag{28} \end{equation} $$

implying that the local truncation error in the position is now $O(h^3)$, whereas Euler's or Euler-Cromer's methods have a local error of $O(h^2)$.

Differential Equations

Thus, the midpoint method yields a global error with second-order accuracy for the position and first-order accuracy for the velocity. However, although these methods yield exact results for constant accelerations, the error increases in general with each time step.

One method that avoids this is the so-called half-step method. Here we define

$$ \begin{equation} y^{(2)}_{n+1/2}=y^{(2)}_{n-1/2}+h a_{n}+O(h^2), \label{_auto22} \tag{29} \end{equation} $$

and

$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \label{_auto23} \tag{30} \end{equation} $$

Note that this method needs the calculation of $y^{(2)}_{1/2}$. This is done using e.g., Euler's method

$$ \begin{equation} y^{(2)}_{1/2}=y^{(2)}_{0}+h a_{0}+O(h^2). \label{_auto24} \tag{31} \end{equation} $$

As this method is numerically stable, it is often used instead of Euler's method.

Differential Equations

Another method which one may encounter is the Euler-Richardson method with

$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n+1/2}+O(h^2), \label{eq:er1} \tag{32} \end{equation} $$

and

$$ \begin{equation} \label{eq:er2} \tag{33} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$

The Verlet method

Another set of popular algorithms, which are both numerically stable and easy to implement are the Verlet algorithms, with the velocity Verlet method as widely used in for example Molecular dynamics calculations. Consider again a second-order differential equation like Newton's second law, whose one-dimensional version reads

$$ m\frac{d^2 x}{dt^2}= F(x,t), $$

which we rewrite in terms of two coupled differential equations

$$ \frac{dx}{dt}=v(x,t) \hspace{1cm}\mathrm{and}\hspace{1cm} \frac{dv}{dt}=F(x,t)/m=a(x,t). $$

The Verlet method

If we now perform a Taylor expansion

$$ x(t+h) = x(t)+hx^{(1)}(t)+\frac{h^2}{2}x^{(2)}(t)+O(h^3). $$

In our case the second derivative is known via Newton's second law, namely $x^{(2)}(t)=a(x,t)$. If we add to the above equation the corresponding Taylor expansion for $x(t-h)$, we obtain, using the discretized expressions

$$ x(t_i\pm h) = x_{i\pm 1} \hspace{1cm}\mathrm{and}\hspace{1cm} x_i=x(t_i), $$

we arrive at

$$ x_{i+1} = 2x_i - x_{i-1} +h^2x^{(2)}_i+O(h^4). $$

The velocity Verlet method

We note that the truncation error goes like $O(h^4)$ since all the odd terms cancel when we add the two Taylor expansions. We see also that the velocity is not directly included in the equation since the function $x^{(2)}=a(x,t)$ is supposed to be known. If we need the velocity however, we can compute it using the well-known formula

$$ x^{(1)}_i=\frac{x_{i+1}-x_{i-1}}{2h}+O(h^2). $$

We note that the velocity has a truncation error which goes like $O(h^2)$. In for example so-called Molecular dynamics calculations, since the acceleration is normally known via Newton's second law, there is seldomly a need for computing the velocity.

We note also that the algorithm for the position is not self-starting since, for $i=0$ it depends on the value of $x$ at the fictitious value $x_{-1}$.

We can amend this by introducing the velocity Verlet method.

The velocity Verlet method

We have the Taylor expansion for the position given by

$$ x_{i+1} = x_i+hx^{(1)}_i+\frac{h^2}{2}x^{(2)}_i+O(h^3). $$

The corresponding expansion for the velocity is

$$ v_{i+1} = v_i+hv^{(1)}_i+\frac{h^2}{2}v^{(2)}_i+O(h^3). $$

Via Newton's second law we have normally an analytical expression for the derivative of the velocity, namely

$$ v^{(1)}_i = \frac{d^2 x}{dt^2}\vert_{i}= \frac{F(x_i,t_i)}{m}, $$

The velocity Verlet method

If we add to this the corresponding expansion for the derivative of the velocity

$$ v^{(1)}_{i+1} = v^{(1)}_i+hv^{(2)}_i+O(h^2), $$

and retain only terms up to the second derivative of the velocity since our error goes as $O(h^3)$, we have

$$ hv^{(2)}_i\approx v^{(1)}_{i+1}-v^{(1)}_i. $$

We can then rewrite the Taylor expansion for the velocity as

$$ v_{i+1} = v_i+\frac{h}{2}\left( v^{(1)}_{i+1}+v^{(1)}_{i}\right)+O(h^3). $$

The velocity Verlet method

Our final equations for the position and the velocity become then

$$ x_{i+1} = x_i+hv_i+\frac{h^2}{2}v^{(1)}_{i}+O(h^3), $$

and

$$ v_{i+1} = v_i+\frac{h}{2}\left( v^{(1)}_{i+1}+v^{(1)}_{i}\right)+O(h^3). $$

Note well that the term $v^{(1)}_{i+1}$ depends on the position at $x_{i+1}$. This means that you need to calculate the position at the updated time $t_{i+1}$ before the computing the next velocity. Note also that the derivative of the velocity at the time $t_i$ used in the updating of the position can be reused in the calculation of the velocity update as well.

Differential Equations, Runge-Kutta methods

Runge-Kutta (RK) methods are based on Taylor expansion formulae, but yield in general better algorithms for solutions of an ODE. The basic philosophy is that it provides an intermediate step in the computation of $y_{i+1}$.

To see this, consider first the following definitions

$$ \begin{equation} \frac{dy}{dt}=f(t,y), \label{_auto25} \tag{34} \end{equation} $$

and

$$ \begin{equation} y(t)=\int f(t,y) dt, \label{_auto26} \tag{35} \end{equation} $$

and

$$ \begin{equation} y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt. \label{_auto27} \tag{36} \end{equation} $$

Differential Equations, Runge-Kutta methods

To demonstrate the philosophy behind RK methods, let us consider the second-order RK method, RK2. The first approximation consists in Taylor expanding $f(t,y)$ around the center of the integration interval $t_i$ to $t_{i+1}$, that is, at $t_i+h/2$, $h$ being the step. Using the midpoint formula for an integral, defining $y(t_i+h/2) = y_{i+1/2}$ and
$t_i+h/2 = t_{i+1/2}$, we obtain

$$ \begin{equation} \int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto28} \tag{37} \end{equation} $$

This means in turn that we have

$$ \begin{equation} y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \label{_auto29} \tag{38} \end{equation} $$

Differential Equations, Runge-Kutta methods

However, we do not know the value of $y_{i+1/2}$. Here comes thus the next approximation, namely, we use Euler's method to approximate $y_{i+1/2}$. We have then

$$ \begin{equation} y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt} = y(t_i) + \frac{h}{2}f(t_i,y_i). \label{_auto30} \tag{39} \end{equation} $$

This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2.

5 6

< < < ! ! M A T H _ B L O C K

$$ \begin{equation} k_2=hf(t_{i+1/2},y_i+k_1/2), \label{_auto32} \tag{41} \end{equation} $$

with the final value

$$ \begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \label{_auto33} \tag{42} \end{equation} $$

The difference between the previous one-step methods is that we now need an intermediate step in our evaluation, namely $t_i+h/2 = t_{(i+1/2)}$ where we evaluate the derivative $f$. This involves more operations, but the gain is a better stability in the solution.

Differential Equations, Runge-Kutta methods

The fourth-order Runge-Kutta, RK4, which we will employ in the solution of various differential equations below, has the following algorithm

5 9

< < < ! ! M A T H _ B L O C K

$$ k_3=hf(t_i+h/2,y_i+k_2/2)\hspace{0.5cm} k_4=hf(t_i+h,y_i+k_3) $$

with the final result

$$ y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right). $$

Thus, the algorithm consists in first calculating $k_1$ with $t_i$, $y_1$ and $f$ as inputs. Thereafter, we increase the step size by $h/2$ and calculate $k_2$, then $k_3$ and finally $k_4$. The global error goes as $O(h^4)$.

Building a code for the solar system, gravitational force and constants

We start with a simpler case first, the Earth-Sun system in two dimensions only. The gravitational force $F_G$ is

$$ F=\frac{GM_{\odot}M_E}{r^2}, $$

where $G$ is the gravitational constant,

$$ M_E=6\times 10^{24}\mathrm{Kg}, $$

the mass of Earth,

$$ M_{\odot}=2\times 10^{30}\mathrm{Kg}, $$

the mass of the Sun and

$$ r=1.5\times 10^{11}\mathrm{m}, $$

is the distance between Earth and the Sun. The latter defines what we call an astronomical unit AU. From Newton's second law we have then for the $x$ direction

$$ \frac{d^2x}{dt^2}=\frac{F_{x}}{M_E}, $$

and

$$ \frac{d^2y}{dt^2}=\frac{F_{y}}{M_E}, $$

for the $y$ direction.

Building a code for the solar system, force equations

Introducing $x=r\cos{(\theta)}$, $y=r\sin{(\theta)}$ and

$$ r = \sqrt{x^2+y^2}, $$

we can rewrite

$$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$

and

$$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$

for the $y$ direction.

Building a code for the solar system, coupled equations

We can rewrite these two equations

$$ F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x, $$

and

$$ F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y, $$

as four first-order coupled differential equations

7 3

< < < ! ! M A T H _ B L O C K

7 4

< < < ! ! M A T H _ B L O C K

7 5

< < < ! ! M A T H _ B L O C K

$$ \frac{dy}{dt}=v_y. $$

Building a code for the solar system, final coupled equations

The four coupled differential equations

7 7

< < < ! ! M A T H _ B L O C K

7 8

< < < ! ! M A T H _ B L O C K

7 9

< < < ! ! M A T H _ B L O C K

$$ \frac{dy}{dt}=v_y, $$

can be turned into dimensionless equations (as we did in project 2) or we can introduce astronomical units with $1$ AU = $1.5\times 10^{11}$.

Using the equations from circular motion (with $r =1\mathrm{AU}$)

$$ \frac{M_E v^2}{r} = F = \frac{GM_{\odot}M_E}{r^2}, $$

we have

$$ GM_{\odot}=v^2r, $$

and using that the velocity of Earth (assuming circular motion) is $v = 2\pi r/\mathrm{yr}=2\pi\mathrm{AU}/\mathrm{yr}$, we have

$$ GM_{\odot}= v^2r = 4\pi^2 \frac{(\mathrm{AU})^3}{\mathrm{yr}^2}. $$

Building a code for the solar system, discretized equations

The four coupled differential equations can then be discretized using Euler's method as (with step length $h$)

8 4

< < < ! ! M A T H _ B L O C K

8 5

< < < ! ! M A T H _ B L O C K

8 6

< < < ! ! M A T H _ B L O C K

$$ y_{i+1}=y_i+hv_{y,i}, $$

Building a code for the solar system, adding Jupiter

It is rather straightforward to add a new planet, say Jupiter. Jupiter has mass

$$ M_J=1.9\times 10^{27}\mathrm{kg}, $$

and distance to the Sun of $5.2$ AU. The additional gravitational force the Earth feels from Jupiter in the $x$-direction is

$$ F_{x}^{EJ}=-\frac{GM_JM_E}{r_{EJ}^3}(x_E-x_J), $$

where $E$ stands for Earth, $J$ for Jupiter, $r_{EJ}$ is distance between Earth and Jupiter

$$ r_{EJ} = \sqrt{(x_E-x_J)^2+(y_E-y_J)^2}, $$

and $x_E$ and $y_E$ are the $x$ and $y$ coordinates of Earth, respectively, and $x_J$ and $y_J$ are the $x$ and $y$ coordinates of Jupiter, respectively. The $x$-component of the velocity of Earth changes thus to

$$ \frac{dv_x^E}{dt}=-\frac{GM_{\odot}}{r^3}x_E-\frac{GM_J}{r_{EJ}^3}(x_E-x_J). $$

Building a code for the solar system, adding Jupiter

We can rewrite

$$ \frac{dv_x^E}{dt}=-\frac{GM_{\odot}}{r^3}x_E-\frac{GM_J}{r_{EJ}^3}(x_E-x_J). $$

to

$$ \frac{dv_x^E}{dt}=-\frac{4\pi^2}{r^3}x_E-\frac{4\pi^2M_J/M_{\odot}}{r_{EJ}^3}(x_E-x_J), $$

where we used

$$ GM_J = GM_{\odot}\left(\frac{M_J}{M_{\odot}}\right)=4\pi^2 \frac{M_J}{M_{\odot}}. $$

Similarly, for the velocity in $y$-direction we have

$$ \frac{dv_y^E}{dt}=-\frac{4\pi^2}{r^3}y_E-\frac{4\pi^2M_J/M_{\odot}}{r_{EJ}^3}(y_E-y_J). $$

Similar expressions apply for Jupiter. The equations for $x$ and $y$ derivatives are unchanged. This equations are similar for all other planets and as we will see later, it will be convenient to object orient this part when we program the full solar system.

How can I get the initial velocities and positions of the planets

NASA has an excellent site at http://ssd.jpl.nasa.gov/horizons.cgi#top. From there you can extract initial conditions in order to start your differential equation solver. At the above website you need to change from OBSERVER to VECTOR and then write in the planet you are interested in. The generated data contain the $x$, $y$ and $z$ values as well as their corresponding velocities. The velocities are in units of AU per day. Alternatively they can be obtained in terms of km and km/s.

For the first simple system involving the Earth and the Sun, you could just initialize the position with say $x=1$ AU and $y=0$ AU.

How do we develop a larger numerical project

When building up a numerical project there are several elements you should think of, amongst these we take the liberty of mentioning the following.

1. How to structure a code in terms of functions and modules

2. How to read input data flexibly from the command line

3. Find tests and how to write unit tests (test functions). A very good example is Catch, a modern, C++-native, header-only, framework for unit-tests. Try to find suitable tests of the mathematical algorithms as well as tests which reflect the physics of the system.

4. It takes one quick command to let all your code undergo heavy testing

5. How to refactor code in terms of classes (instead of functions only)

6. How to conduct and automate large-scale numerical experiments

More on large numerical projects

The conventions and techniques outlined here will save you a lot of time when you incrementally extend software over time from simpler to more complicated problems.

1. Scale your equations in order to simplify, make for example the equations dimensionless or scale them in convenient units

2. New code is added in a modular fashion to a library (modules)

3. Programs are run through convenient user interfaces

4. Tedious manual work with running programs is automated,

5. Your scientific investigations are reproducible, scientific reports with top quality typesetting are produced both for paper and electronic devices.

A few words on unit testing

A very good and simple framework for unit testing is Catch. Unit Testing is the practice of testing the smallest testable parts, called units, of an application individually and independently to determine if they behave exactly as expected. Unit tests (short code fragments) are usually written such that they can be performed at any time during the development to continually verify the behavior of the code. In this way, possible bugs will be identified early in the development cycle, making the debugging at later stage much easier.

Unit tests with Catch

There are many benefits associated with Unit Testing since it increases confidence in changing and maintaining code. Big changes can be made to the code quickly, since the tests will ensure that everything still is working properly.

• Since the code needs to be modular to make Unit Testing possible, the code will be easier to reuse. This improves the code design.

• Debugging is easier, since when a test fails, only the latest changes need to be debugged.

• Different parts of a project can be tested without the need to wait for the other parts to be available.

• A unit test can serve as a documentation on the functionality of a unit of the code.

A simple Example first, Block tied to a Wall

Our first example is the classical case of simple harmonic oscillations, namely a block sliding on a horizontal frictionless surface. The block is tied to a wall with a spring. If the spring is not compressed or stretched too far, the force on the block at a given position $x$ is

$$ F=-kx. $$

The negative sign means that the force acts to restore the object to an equilibrium position. Newton's equation of motion for this idealized system is then

$$ m\frac{d^2x}{dt^2}=-kx, $$

or we could rephrase it as

$$ \frac{d^2x}{dt^2}=-\frac{k}{m}x=-\omega_0^2x, \label{eq:newton1} \tag{43} $$

with the angular frequency $\omega_0^2=k/m$.

The above differential equation has the advantage that it can be solved analytically with solutions on the form

$$ x(t)=Acos(\omega_0t+\nu), $$

where $A$ is the amplitude and $\nu$ the phase constant. This provides in turn an important test for the numerical solution and the development of a program for more complicated cases which cannot be solved analytically.

Simple Example, Block tied to a Wall

With the position $x(t)$ and the velocity $v(t)=dx/dt$ we can reformulate Newton's equation in the following way

$$ \frac{dx(t)}{dt}=v(t), $$

and

$$ \frac{dv(t)}{dt}=-\omega_0^2x(t). $$

We are now going to solve these equations using the Runge-Kutta method to fourth order discussed previously.

Simple Example, Block tied to a Wall

Before proceeding however, it is important to note that in addition to the exact solution, we have at least two further tests which can be used to check our solution.

Since functions like $cos$ are periodic with a period $2\pi$, then the solution $x(t)$ has also to be periodic. This means that

$$ x(t+T)=x(t), $$

with $T$ the period defined as

$$ T=\frac{2\pi}{\omega_0}=\frac{2\pi}{\sqrt{k/m}}. $$

Observe that $T$ depends only on $k/m$ and not on the amplitude of the solution.

Simple Example, Block tied to a Wall

In addition to the periodicity test, the total energy has also to be conserved.

Suppose we choose the initial conditions

$$ x(t=0)=1\hspace{0.1cm} \mathrm{m}\hspace{1cm} v(t=0)=0\hspace{0.1cm}\mathrm{m/s}, $$

meaning that block is at rest at $t=0$ but with a potential energy

$$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$

The total energy at any time $t$ has however to be conserved, meaning that our solution has to fulfil the condition

$$ E_0=\frac{1}{2}kx(t)^2+\frac{1}{2}mv(t)^2. $$

Simple Example, Block tied to a Wall

An algorithm which implements these equations is included below.

• Choose the initial position and speed, with the most common choice $v(t=0)=0$ and some fixed value for the position.

• Choose the method you wish to employ in solving the problem.

• Subdivide the time interval $[t_i,t_f] $ into a grid with step size

$$ h=\frac{t_f-t_i}{N}, $$

where $N$ is the number of mesh points.

• Calculate now the total energy given by

$$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$

• The Runge-Kutta method is used to obtain $x_{i+1}$ and $v_{i+1}$ starting from the previous values $x_i$ and $v_i$.

  • When we have computed $x(v)_{i+1}$ we upgrade $t_{i+1}=t_i+h$.

  • This iterative process continues till we reach the maximum time $t_f$.

  • The results are checked against the exact solution. Furthermore, one has to check the stability of the numerical solution against the chosen number of mesh points $N$.

Simple Example, Block tied to a Wall

Simple Example, Block tied to a Wall, python code

The following python program performs essentially the same calculations as the previous c++ code.

The second example: The classical pendulum and scaling the equations

The angular equation of motion of the pendulum is given by Newton's equation and with no external force it reads

$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \label{_auto34} \tag{44} \end{equation} $$

with an angular velocity and acceleration given by

$$ \begin{equation} v=l\frac{d\theta}{dt}, \label{_auto35} \tag{45} \end{equation} $$

and

$$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \label{_auto36} \tag{46} \end{equation} $$

More on the Pendulum

We do however expect that the motion will gradually come to an end due a viscous drag torque acting on the pendulum. In the presence of the drag, the above equation becomes

$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \label{eq:pend1} \tag{47} \end{equation} $$

where $\nu$ is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then

$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \label{eq:pend2} \tag{48} \end{equation} $$

with $A$ and $\omega$ two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.

More on the Pendulum

We define

$$ \omega_0=\sqrt{g/l}, $$

the so-called natural frequency and the new dimensionless quantities

$$ \hat{t}=\omega_0t, $$

with the dimensionless driving frequency

$$ \hat{\omega}=\frac{\omega}{\omega_0}, $$

and introducing the quantity $Q$, called the quality factor,

$$ Q=\frac{mg}{\omega_0\nu}, $$

and the dimensionless amplitude

$$ \hat{A}=\frac{A}{mg} $$

More on the Pendulum

We have

$$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$

This equation can in turn be recast in terms of two coupled first-order differential equations as follows

$$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$

and

$$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$

These are the equations to be solved. The factor $Q$ represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude $\hat{A}$ is measured in units of the maximum possible gravitational torque while $\hat{\omega}$ is the angular frequency of the external torque measured in units of the pendulum's natural frequency.

How to start our solar system project

The simplest possible step is to code the Earth-Sun system using Euler's method in two dimensions. The four coupled differential equations can then be discretized using Euler's method as (with step length $h$) are

1 2 2

< < < ! ! M A T H _ B L O C K

1 2 3

< < < ! ! M A T H _ B L O C K

1 2 4

< < < ! ! M A T H _ B L O C K

$$ y_{i+1}=y_i+hv_{y,i}, $$

How to start our solar system project, the simple Euler equations

    /*    
        Euler's method for the Eart-Sun system, simplest possible code
    */ 
    #include <cmath>
    #include <iostream>
    #include <fstream>
    #include <iomanip>
    using namespace  std;
    // output file as global variable
    ofstream ofile;
    // function declarations
    void output( double, double, double, double, double);

    int main(int argc, char* argv[])
    {
    //  declarations of variables
      char *outfilename;
      // Read in output file, abort if there are too few command-line arguments
      if( argc <= 1 ){
        cout << "Bad Usage: " << argv[0] <<
          " read also output file on same line" << endl;
        //    exit(1);
      }
      else{
        outfilename=argv[1];
      }
      ofile.open(outfilename);

      int NumberofSteps = 1000;
      double FinalTime = 1.0;
      double Step = FinalTime/((double) NumberofSteps);
      double time = 0.0;
      // Initial values  x = 1.0 AU and vy = 2*pi
      double pi = acos(-1.0);
      double FourPi2 = 4*pi*pi;
      double x =  1.0; double y =  0.0; double vx = 0.0; double vy = 2.0*pi;
      double r = sqrt(x*x+y*y);

    // now we start solving the differential equations using Euler's method
      while (time <= FinalTime){
        x += Step*vx;
        y += Step*vy;
        vx -= Step*FourPi2*x/(r*r*r);
        vy -= Step*FourPi2*y/(r*r*r);
        r = sqrt(x*x+y*y);
        time += Step;
        output(time, x, y, vx, vy);   // write to file 
      }
      ofile.close();  // close output file
      return 0;
    }   //  End of main function 

    //    function to write out the final results
    void output(double time, double x, double y, double vx, double vy)
    {
      ofile << setiosflags(ios::showpoint | ios::uppercase);
      ofile << setw(15) << setprecision(8) << time;
      ofile << setw(15) << setprecision(8) << x;
      ofile << setw(15) << setprecision(8) << y;
      ofile << setw(15) << setprecision(8) << vx;
      ofile << setw(15) << setprecision(8) << vy << endl;
    }  // end of function output

Where do we go from here?

• Our next step could be to simply replace the brute force coding of the positions and velocities using the vector class discussed in the solar system class example. That is, we continue with the above Euler code and simply replace the variables $x$, $y$, $vx$ and $vy$ with vectors defined in the vec3 class.

• It is important to have a simple working code that is well tested and understood before we move on to new additions. In this part it is useful to introduce possible unit tests.

• The next step is to add a class about for example celestial bodies, as discussed in the simple extension of the Euler method for the Earth-Sun system. This will discussed during the lab sessions.

• Thereafter we can add the full complexity as discussed in the solar system class example to be discussed during the lab sessions.

Adaptive methods

In case the function to integrate varies slowly or fast in different integration domains, adaptive methods are normally used. One strategy is always to decrease the step size. As we have seen earlier, this leads to more CPU cycles and may lead to loss or numerical precision. An alternative is to use higher-order RK methods for example. However, this leads again to more cycles, furthermore, there is no guarantee that higher-order leads to an improved error.

Adaptive methods

Assume the exact result is $\tilde{x}$ and that we are using an RKM method. Suppose we run two calculations, one with $h$ (called $x_1$) and one with $h/2$ (called $x_2$). Then

$$ \tilde{x}=x_1+Ch^{M+1}+O(h^{M+2}), $$

and

$$ \tilde{x}=x_2+2C(h/2)^{M+1}+O(h^{M+2}), $$

with $C$ a constant. Note that we calculate two halves in the last equation. We get then

$$ |x_1-x_2| = Ch^{M+1}(1-\frac{1}{2^M}). $$

yielding

$$ C=\frac{|x_1-x_2|}{(1-2^{-M})h^{M+1}}. $$

We rewrite

$$ \tilde{x}=x_2+\epsilon+O((h)^{M+2}), $$

with

$$ \epsilon = \frac{|x_1-x_2|}{2^M-1}. $$

Adaptive methods

With RK4 the expressions become

$$ \tilde{x}=x_2+\epsilon+O((h)^{6}), $$

with

$$ \epsilon = \frac{|x_1-x_2|}{15}. $$

The estimate is one order higher than the original RK4. But this method is normally rather inefficient since it requires a lot of computations. We solve typically the equation three times at each time step. However, we can compare the estimate $\epsilon$ with some by us given accuracy $\xi$. We can then ask the question: what is, with a given $x_j$ and $t_j$, the largest possible step size $\tilde{h}$ that leads to a truncation error below $\xi$? We want

$$ C\tilde{h} \le \xi, $$

which leads to

$$ \left(\frac{\tilde{h}}{h}\right)^{M+1}\frac{|x_1-x_2|}{(1-2^{-M})}\le \xi, $$

meaning that

$$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$

Adaptive methods

With

$$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$

we can design the following algorithm:

• If the two answers are close, keep the approximation to $h$.

• If $\epsilon > \xi$ we need to decrease the step size in the next time step.

• If $\epsilon < \xi$ we need to increase the step size in the next time step.