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

Date: **Oct 5, 2017**

Copyright 1999-2017, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license

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.

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}
$$

$$
\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}
$$

$$
\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.

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.

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.

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}
$$

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

$$
\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.

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}
$$

$$
\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}
$$

$$
\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)$.

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),
$$

$$
\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}
$$

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}
$$

$$
\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}
$$

$$
\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.

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)$.

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}
$$

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

$$
\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}
$$

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).
$$

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

$$
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).
$$

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.

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).
$$

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

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

$$
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).
$$

$$
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.

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}
$$

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}
$$

$$
\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.

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)$.

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.
$$

$$
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.

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.

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}}.
$$

$$
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.
$$

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

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$.

```
In [1]:
```/* This program solves Newton's equation for a block
sliding on a horizontal frictionless surface. The block
is tied to a wall with a spring, and Newton's equation
takes the form
m d^2x/dt^2 =-kx
with k the spring tension and m the mass of the block.
The angular frequency is omega^2 = k/m and we set it equal
1 in this example program.
Newton's equation is rewritten as two coupled differential
equations, one for the position x and one for the velocity v
dx/dt = v and
dv/dt = -x when we set k/m=1
We use therefore a two-dimensional array to represent x and v
as functions of t
y[0] == x
y[1] == v
dy[0]/dt = v
dy[1]/dt = -x
The derivatives are calculated by the user defined function
derivatives.
The user has to specify the initial velocity (usually v_0=0)
the number of steps and the initial position. In the programme
below we fix the time interval [a,b] to [0,2*pi].
*/
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
//#include "lib.h"
using namespace std;
// output file as global variable
ofstream ofile;
// function declarations
void derivatives(double, double *, double *);
void initialise ( double&, double&, int&);
void output( double, double *, double);
void runge_kutta_4(double *, double *, int, double, double,
double *, void (*)(double, double *, double *));
int main(int argc, char* argv[])
{
// declarations of variables
double *y, *dydt, *yout, t, h, tmax, E0;
double initial_x, initial_v;
int i, number_of_steps, n;
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);
// this is the number of differential equations
n = 2;
// allocate space in memory for the arrays containing the derivatives
dydt = new double[n];
y = new double[n];
yout = new double[n];
// read in the initial position, velocity and number of steps
initialise (initial_x, initial_v, number_of_steps);
// setting initial values, step size and max time tmax
h = 4.*acos(-1.)/( (double) number_of_steps); // the step size
tmax = h*number_of_steps; // the final time
y[0] = initial_x; // initial position
y[1] = initial_v; // initial velocity
t=0.; // initial time
E0 = 0.5*y[0]*y[0]+0.5*y[1]*y[1]; // the initial total energy
// now we start solving the differential equations using the RK4 method
while (t <= tmax){
derivatives(t, y, dydt); // initial derivatives
runge_kutta_4(y, dydt, n, t, h, yout, derivatives);
for (i = 0; i < n; i++) {
y[i] = yout[i];
}
t += h;
output(t, y, E0); // write to file
}
delete [] y; delete [] dydt; delete [] yout;
ofile.close(); // close output file
return 0;
} // End of main function
// Read in from screen the number of steps,
// initial position and initial speed
void initialise (double& initial_x, double& initial_v, int& number_of_steps)
{
cout << "Initial position = ";
cin >> initial_x;
cout << "Initial speed = ";
cin >> initial_v;
cout << "Number of steps = ";
cin >> number_of_steps;
} // end of function initialise
// this function sets up the derivatives for this special case
void derivatives(double t, double *y, double *dydt)
{
dydt[0]=y[1]; // derivative of x
dydt[1]=-y[0]; // derivative of v
} // end of function derivatives
// function to write out the final results
void output(double t, double *y, double E0)
{
ofile << setiosflags(ios::showpoint | ios::uppercase);
Ofile << setw(15) << setprecision(8) << t;
ofile << setw(15) << setprecision(8) << y[0];
ofile << setw(15) << setprecision(8) << y[1];
ofile << setw(15) << setprecision(8) << cos(t);
ofile << setw(15) << setprecision(8) <<
0.5*y[0]*y[0]+0.5*y[1]*y[1]-E0 << endl;
} // end of function output
/* This function upgrades a function y (input as a pointer)
and returns the result yout, also as a pointer. Note that
these variables are declared as arrays. It also receives as
input the starting value for the derivatives in the pointer
dydx. It receives also the variable n which represents the
number of differential equations, the step size h and
the initial value of x. It receives also the name of the
function *derivs where the given derivative is computed
*/
void runge_kutta_4(double *y, double *dydx, int n, double x, double h,
double *yout, void (*derivs)(double, double *, double *))
{
int i;
double xh,hh,h6;
double *dym, *dyt, *yt;
// allocate space for local vectors
dym = new double [n];
dyt = new double [n];
yt = new double [n];
hh = h*0.5;
h6 = h/6.;
xh = x+hh;
for (i = 0; i < n; i++) {
yt[i] = y[i]+hh*dydx[i];
}
(*derivs)(xh,yt,dyt); // computation of k2, eq. 3.60
for (i = 0; i < n; i++) {
yt[i] = y[i]+hh*dyt[i];
}
(*derivs)(xh,yt,dym); // computation of k3, eq. 3.61
for (i=0; i < n; i++) {
yt[i] = y[i]+h*dym[i];
dym[i] += dyt[i];
}
(*derivs)(x+h,yt,dyt); // computation of k4, eq. 3.62
// now we upgrade y in the array yout
for (i = 0; i < n; i++){
yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
}
delete []dym;
delete [] dyt;
delete [] yt;
} // end of function Runge-kutta 4

```
In [2]:
```#
# This program solves Newtons equation for a block sliding on
# an horizontal frictionless surface.
# The block is tied to the wall with a spring, so N's eq takes the form:
#
# m d^2x/dt^2 = - kx
#
# In order to make the solution dimless, we set k/m = 1.
# This results in two coupled diff. eq's that may be written as:
#
# dx/dt = v
# dv/dt = -x
#
# The user has to specify the initial velocity and position,
# and the number of steps. The time interval is fixed to
# t \in [0, 4\pi) (two periods)
#
# Note that this is a highly simplifyed rk4 code, intended
# for conceptual understanding and experimentation.
import sys
import numpy, math
#Global variables
ofile = None;
E0 = 0.0
def sim(x_0, v_0, N):
ts = 0.0
te = 4*math.pi
h = (te-ts)/float(N)
t = ts;
x = x_0
v = v_0
while (t < te):
kv1 = -h*x
kx1 = h*v
kv2 = -h*(x+kx1/2)
kx2 = h*(v+kv1/2)
kv3 = -h*(x+kx2/2)
kx3 = h*(v+kv2/2)
kv4 = -h*(x+kx3/2)
kx4 = h*(v+kv3/2)
#Write the old values to file
output(t,x,v)
#Update
x = x + (kx1 + 2*(kx2+kx3) + kx4)/6
v = v + (kv1 + 2*(kv2+kv3) + kv4)/6
t = t+h
def output(t,x,v):
de = 0.5*x**2+0.5*v**2 - E0;
ofile.write("%15.8E %15.8E %15.8E %15.8E %15.8E\n"\
%(t, x, v, math.cos(t),de));
#MAIN PROGRAM:
#Get input
if len(sys.argv) == 5:
ofilename = sys.argv[1];
x_0 = float(sys.argv[2])
v_0 = float(sys.argv[3])
N = int(sys.argv[4])
else:
print "Usage:", sys.argv[0], "ofilename x0 v0 N"
sys.exit(0)
#Setup
ofile = open(ofilename, 'w')
E0 = 0.5*x_0**2+0.5*v_0**2
#Run simulation
sim(x_0,v_0,N)
#Cleanup
ofile.close()

$$
\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}
$$

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

$$
\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}
$$

$$
\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}
$$

$$
\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}).
$$

$$
\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.

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},
$$

**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},
$$

$$
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,
$$

$$
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

9 9

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

1 0 0

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

1 0 1

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

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

1 0 3

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

1 0 4

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

1 0 5

< < < ! ! 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,
$$

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

1 1 0

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

1 1 1

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

1 1 2

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

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

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

$$
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},
$$

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

$$
\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.

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.

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

How to structure a code in terms of functions and modules

How to read input data flexibly from the command line

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.

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

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

How to conduct and automate large-scale numerical experiments

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.

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

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

Programs are run through convenient user interfaces

Tedious manual work with running programs is automated,

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

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},
$$

```
/*
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
```

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.

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}.
$$

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

with

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

$$
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}.
$$

$$
\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.