Support material

Support material#

Background#

The problem of solving differential equations is not just a mathematical concern. Many physical equations that describe the evolution and equilibrium of very diverse systems are differential equations. There are two types of differential equations:

Ordinary Differential Equations (ODE) in which the function being sought depends only on one variable, usually time \(t\).
Partial Differential Equations in which it also depends on other variables, usually space.

In this course, we are only interested in ordinary differential equations. PDEs are usually solved by the finite element method which requires a course on its own. The notions discussed and the observations made are also largely applicable to partial differential equations.

Note

This course is not a background paper on the underlying mathematical concepts. It is intended to provide the scientist with the necessary notions to make choices regarding the treatment of ODEs.

Problem formulation#

An ODE can be described by several characteristics. In particular:

Its order.
Its number of equations: depending on the problem, only one here but many ODEs can be presented as systems of equations.
The explicit time dependence. If this one does not appear, the equation is said to be autonomous.
Linearity: if the unknown function and its derivatives appear only multiplied by constants, then the equation is linear.

On the contrary, for example if they appear to a power other than 1 or in a sine, then the equation is non-linear. Here the equation is linear.

Example

As an example, we can consider the case of a free and undamped oscillator. Its motion is described by the following equation:

(1)#\[\ddot x + \omega_0^2 x = 0\]

The equation is linear of the second order. Its exact solution is

\[x(t) = a \cos (\omega_0 t)\]

where \(a \in \mathbb R^+\) is any constant. There is therefore an infinity of solutions.

Depending on these characteristics, an equation will or will not have an exact known solution. Only a minority of them have such a solution. Although they are not the subject of this course, the exact solutions will serve as a reference to test the approximate solutions.

To address ODEs on a systematic basis, a generic formulation is usually used that takes the form:

(2)#\[\dot X = f(X, t)\]

Where:

\(X: t \mapsto X(t) = [X_0(t),\dots, X_{N-1}(t)]^T \) is the unknown function.
\(f\) is the ODE

This formulation is both simple and surprising, it implies that all ODEs can be written as first order equations. Of course, this advantage has a price: the decrease of the order goes through an increase of the number of equations through the creation of intermediate functions.

Example: ODE reformulation

From a practical point of view, formulating a problem in this way can take a bit of practice. Let’s take the example of the equation (1).

The unknown function \(x\) can be replaced by a vector function:

\[\begin{split}X = \begin{bmatrix} X[0] \\ X[1] \end{bmatrix} = \begin{bmatrix} x \\ \dot x \end{bmatrix}\end{split}\]

As a consequence:

\[\begin{split}\dot X = \begin{bmatrix} \dot x \\ \ddot x \end{bmatrix}\end{split}\]

The equation is rewritten as

(3)#\[\begin{split}\underbrace{\begin{bmatrix} \dot x \\ \ddot x \end{bmatrix}}_{\dot X} = \underbrace{\begin{bmatrix} X[1] \\ - \omega_0^2 X[0] \end{bmatrix}}_{f(X, t)}\end{split}\]

Phase space, vector flow and integral curves#

This reformulation is not just a mathematical trick. It also informs us about the physical nature of the problem being modeled. Indeed, the vector \(X\) describes the state variables of the system. Thus, if \(X\) and \(t\) are known, the state of the system is fully determined. This has several essential consequences:

The solution of an ODE is related to the choice of initial conditions of \(X\).
One can therefore plot the solutions of an equation in the \(X\) space called phase space.
The ODE describes the evolution of \(X\) at any point in phase space and at any time. By analogy, we can see \(X\) as a position and \(\dot X = f(X, t)\) as the velocity of a particle in the phase space. The vector field defined by \(f\) thus behaves like the velocity field associated with a fluid flow, so it is called the vector flow.
The solution of an ODE is called the integral curve of the flow and behaves like a stream-line of it.

Note

From a mathematical point of view, the Picard-Lindelöf (a. k. a. Cauchy-Lipschitz) theorem proves the existence and uniqueness of the solution if the equation \(f\) is continuous.

Example: Phase plane and Vector flow

We propose to construct the phase plane and the vector stream associated with the undamped free oscillator (see below). We note that this one is a vortex centered on the origin of the plane. We can therefore immediately guess that the solutions will be closed curves, which implies that they are oscillatory in nature.

Approximate solutions#

Background#

The solution of an ODE in the general case is not known. However, we can calculate numerical solutions. These are fundamentally different from the exact solutions. Suppose that we look for the solution from a starting point \(X_0\) and on a time interval \([t_a, t_b]\):

The numerical solutions are approximate, i.e. they imply a certain amount of error which must be made negligible depending on the problem being addressed.
They are discrete in time, which means that the solution is described only at certain times \(\left\lbrace t_0, \ldots, t_{N-1} \right\rbrace\) where \(t_0 = t_a\) and \(t_{N-1}= t_b\). The distance between these instants is called time step and noted \(dt = t_{k+1} - t_k\). The values of X associated to these time steps are noted \(\left\lbrace X_0, \ldots, X_{N-1} \right\rbrace\)

Problem formulation#

To solve this problem, we can use integration schemes to determine \(X_{k+1}\) from \(X_k\). From a starting point \(X_0\), they allow to determine all the \(X_k\) values step by step on the desired duration. From a mathematical point of view, the problem is therefore to determine:

\[X_{k+1} = X_k + \int_{t=t_k}^{t_{k+1}} f(X(t), t) dt \]

This integral is of course not known and must be approximated efficiently.

Naive solution: Forward Euler#

The simplest approach to solve this problem is to evaluate \(\dot X_k\), the slope of in \((X_k, t_k)\) and follow this direction to \(t_k\). This leads to the following integration scheme:

\[X_{k+1} = X_k + f(X_k, t_k) dt\]

This approach means that:

\[\int_{t=t_k}^{t_{k+1}} f(X(t), t) dt \approx f(X_k,t_k) dt\]

This is equivalent to approximating the integral by the method of left rectangles (or left Riemann sum). This method of integration is called the explicit Euler method or the forward Euler method. It is attributed to Leonhard Euler. This method has the merit of being easy to write:

def Euler(func, X0, t):
    """
    Euler solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt - 1):
        X[i + 1] = X[i] + func(X[i], t[i]) * dt
    return X

The above example shows that the error between the exact solution and the approximate solution is clearly visible and is probably not negligible. We can reduce the time step \(dt\) to limit this error.

Reducing \(dt\) does give a better approximation. We say that the algorithm converges to the solution when \(dt\) tends to 0. This rate of convergence can be shown by measuring the error:

\[e(dt) = \sqrt{\dfrac{1}{k}\sum_k (x_k - x_s(t_k))^2}\]

Where \(x_s\) is the exact solution and \(e\) the error. Let’s plot this error as a function of \(dt\):

We see that Euler’s method converges to the exact solution when \(dt\) decreases. We can therefore chose the time step to limit the error to a value of our choice.

Runge-Kutta 4 method#

There are many improvements of the Euler forward method. The Runge and Kutta methods are an example widely used in the literature, in particular the Runge-Kutta 4 (forward) algorithm which is written as follows:

\[\begin{split}\begin{align} k_1 & = f(X_k, t_k) \\ k_2 & = f(X_k + \dfrac{k_1}{2}dt, t_k + \dfrac{dt}{2}) \\ k_3 & = f(X_k + \dfrac{k_2}{2}dt, t_k + \dfrac{dt}{2}) \\ k_4 & = f(X_k + k_3dt, t_k + dt) \\ X_{k+1} & = X_{k} + \dfrac{1}{6}\left( k_1 + 2 k_2 + 2k_3 + k_4 \right) \end{align}\end{split}\]

def RK4(func, X0, t):
    """
    Runge Kutta 4 solver.
    """
    dt = t[1] - t[0]
    nt = len(t)
    X = np.zeros([nt, len(X0)])
    X[0] = X0
    for i in range(nt - 1):
        k1 = func(X[i], t[i])
        k2 = func(X[i] + dt / 2.0 * k1, t[i] + dt / 2.0)
        k3 = func(X[i] + dt / 2.0 * k2, t[i] + dt / 2.0)
        k4 = func(X[i] + dt * k3, t[i] + dt)
        X[i + 1] = X[i] + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    return X

Built-in methods#

Also built-in integration methods are available. The odeint function of Scipy is an example. We propose to compare it with the two others.

from scipy.integrate import odeint

The graph above shows that the error levels and convergence rates of the 3 methods above are very different. We can first compare the forward Euler and RK4 methods. The latter evaluates the ODE 4 times per iteration and can therefore be considered as four times slower that the Forward Euler method. However, its error level is several orders of magnitude lower and its convergence rate is much better than forward Euler. Concerning odeint, we notice that this last algorithm has a much lower error level than Euler but a comparable convergence rate. In fact it uses an error control system which we will not detail here. However, we must remember that odeint is fast and reliable in many cases.

Support material

Contents

Support material#

Background#

Problem formulation#

Phase space, vector flow and integral curves#

Approximate solutions#

Background#

Problem formulation#

Naive solution: Forward Euler#

Runge-Kutta 4 method#

Built-in methods#