10 Ordinary Differential Equations (ODEs)

19 Feb 2019

Most physical phenomena are ultimately described by a relationship between changing quantities, resulting in differential equations. If such an equation only contains one independent variable (such as time) and hence only total derivatives (and no partial derivatives) we classify it as an ordinary differential equation (ODE). An ODE of order $n$ contains no derivatives higher than the $n$-th derivative:

$$F(t, y, y^{(1)}, \dots, y^{(n)}) = 0$$

The dependent variable $y = y(t)$ is a function of the independent variable $t$ and $y^{(n)} = \frac{d^{n}y}{dt^{n}}$ is the $n$-th derivative with respect to $t$. A linear ODE only contains first powers of $y$ or its derivatives. A non-linear ODE may contain higher powers. There often exist methods to solve linear ODEs analytically but this is impossible for most non-linear ODEs. Solutions of an ODE are fixed by the initial conditions¹, e.g., $y_0 = y(t_{0})$ and similar for all higher derivatives. For an ODE of order $n$, exactly $n$ initial conditions are needed.

Using ODE integration algorithms (integrators) we can solve linear and non-linear ODEs of any order numerically. The basic idea is to start with the initial conditions and then propagate $y$ for a small step $h$ to numerically compute $y(t_0 + h)$ and all its derivatives. By repeating the process, one extrapolates from the initial condition to any "later" value of $t$.

We will first study the simplest integrator, the Euler algorithm, and then look at a more robust class of integrators, the Runge-Kutta schemes.

Introductory example: The bouncing ball

We will introduce the basic ideas for solving ODEs by looking at a very simple physical example: the bouncing ball with Newton's equations of motion

$$\frac{d^2 y}{dt^2} = -g$$

where $g$ is the constant acceleration due to gravity² and $y(t)$ is the position of the ball as a function of time (its trajectory).

The forward Euler scheme for any first order ODE

$$\frac{dy}{dt} = f(y, t)$$

is

$$y(t + h) = y(t) + h f(y(t), t).$$

In order to solve the original 2nd order equation of motion we make use of the fact that one $n$-th order ODE can be written as $n$ coupled first order ODEs, namely

$$\begin{align} \frac{dy}{dt} &= v\\ \frac{dv}{dt} &= -g. \end{align}$$

Solve each of the first order ODEs with the Euler algorithm:

$$\begin{align} y(t + h) &= y(t) + h v(t)\\ v(t + h) &= v(t) - h g. \end{align}$$

In class we developed a simple simulation for free fall and for the bouncing ball.

• bounce_ball_numpy.ipynb: simulation and graphing with matplotlib
• bounce_ball_visual.ipynb: visualization with vpython inside the Jupyter notebook using vpython-jupyter

ODEs, standard form, and principles of integrators

The first part of this lesson introduces the formalism and the basic ideas behind integration algorithms. The Jupyter notebook 10-ODEs.ipynb contains the lecture notes.

We then revisit the Euler algorithm and introduce much better algorithms, namely the Runge-Kutta schemes and the Verlet integrator.

Accuracy and performance of Euler and Runge-Kutta integrators

The second part shows how to assess the accuracy of integration algorithms and analyzes in more depth the simple Euler integrator and the Runge-Kutta methods, namely rk2 and rk4. As examples we use the simple 1D harmonic oscillator with potential energy function $U(x) = \frac{1}{2} k x^2$, an anharmonic oscillator ($U(x) = \frac{1}{2} k x^2 (1-\frac{2}{3} \alpha x)$), and a 6-th power oscillator ($U(x) = \frac{1}{6} k x^6$).

The Jupyter notebook 10-ODE_integrators.ipynb contains the lecture notes³. Skeleton code for in-class problem exercises can be found in the notebook 10-ODE_integrators_students.ipynb and in the module integrators.py.⁴

Verlet integrators

The third part introduces a different class of integrators: time reversible (and symplectic) integrators such as the velocity Verlet algorithm, which have worse accuracy than e.g. RK4 but long term stability.

• The notebook 10-ODE-integrators-verlet.ipynb implements the velocity Verlet algorithm and shows how to use it to calculate the orbit of the Earth around the sun.
• 10-ODE-visualize-planets.ipynb shows how to visualize planetary motion in 3D with vpython.

Resources

• Computational Modelling: Chapter 2
• Computational Physics: Chapter 9
• scipy.integrate for high performance integration algorithms for Python, in particular see scipy.integrate.odeint() and scipy.integrate.ode().
• Weisstein, Eric W. "Ordinary Differential Equation." From MathWorld — A Wolfram Web Resource.
• Numerical Recipes in C, WH Press, SA Teukolsky, WT Vetterling, BP Flannery. 2nd ed, 2002. Cambridge University Press. Chapter 16.

---

Footnotes