Differential Equations

Questions to David Rotermund

By Jan Wiersig, modified by Udo Ernst and translated into English by Daniel Harnack. David Rotermund replaced the Matlab code with Python code.

Ordinary Differential Equations

In many fields of science, problems occur where one or several quantities vary dependent on time, e.g. oscillations or decay processes. The description of these processes can be formalized by differential equations. As a simple example, radioactive decay may serve. An amount of radioactive material is give at time $t=0$. We are looking for a function $x(t)$ that gives the amount of the material $x$ which is still present at time $t$. The rate at time $t$ by which the radioactive material decays is proportional to the amount of the material $x(t)$ that exists at time $t$. From these considerations, the following differential equation can be derived:

$\dot{x}=-\alpha x$

with decay constant $\alpha >0$. This differential equation is referred to as \quoting{ordinary}, since it only includes the derivative with respect to one quantity (here time $t$). It is furthermore a differential equation of ‘1. order’, since only the first derivative is taken.

The solution of a differential equation is not a number, but a function. Here, the solution is $x(t)=x_0{e}^{-\alpha t}$, which can be proven by taking the first derivative. Unfortunately, in most practical cases, a differential equation can not be be solved analytically. This holds especially for ordinary systems of differential equations. An example is the Newtonian mechanics of a particle

$m\ddot{\overrightarrow{r}}=\overrightarrow{F}(\overrightarrow{r},\dot{\overrightarrow{r}},t)$

with place $\overrightarrow{r}=(x,y,z)$, velocity $\overrightarrow{v}=\dot{\overrightarrow{r}}$, acceleration $\ddot{\overrightarrow{r}}$, mass $m$, force $\overrightarrow{F}$ and time $t$. This is an ordinary system of differential equations of order 2, since the second derivative with respect to time is taken. Typically one is interested in an initial value problem: Given shall be the place $\overrightarrow{r}(0)$ and the velocity $\overrightarrow{v}(0)$ at the starting time, sought is the place $\overrightarrow{r}(t)$ and velocity $\overrightarrow{v}(t)$ at a later time $t$.

Alternatively to the previous notation, the differential equation can also be interpreted as an ordinary system of differential equations of 1. order with double the number of variables

$\dot{\overrightarrow{r}}=\overrightarrow{v}$

$m\dot{\overrightarrow{v}}=\overrightarrow{F}(\overrightarrow{r},\overrightarrow{v},t).$

Using the ‘phase space variable’ $\overrightarrow{x}=(\overrightarrow{r},\overrightarrow{v})$, it can be transformed to a more compact form

$\dot{\overrightarrow{x}}=\overrightarrow{\tilde{F}}(\overrightarrow{x},t).$

Here, the first three components of $\overrightarrow{\tilde{F}}$ are equal to $\overrightarrow{v}$ and the last three components equal to $\overrightarrow{F}/m$. One can see that it is sufficient to consider ordinary systems of differential equations of 1. order. Note that in the following notation, the tilde $\stackrel{~}{}$ is suppressed.

Euler Method

Figure 8.1.: Flow plan for the Euler integration scheme. The aim is to find the solution of the initial value problem of the differential equation system (DES)

$\dot{\overrightarrow{x}}=\overrightarrow{F}(\overrightarrow{x},t)$

numerically. For this, we first replace the differential quotient by the right sided derivation

$\frac{\overrightarrow{x}(t+\tau )-\overrightarrow{x}(t)}{\tau }+\mathcal{O}(\tau )=\overrightarrow{F}(\overrightarrow{x},t)$

respectively

$\overrightarrow{x}(t+\tau )=\overrightarrow{x}(t)+\overrightarrow{F}(\overrightarrow{x},t)\tau +\mathcal{O}({\tau }^2).$

The truncation error is of 2. order in the step width $\tau$. By following notation we indicate the discretization in time: $t_n=(n-1)\tau$, ${\overrightarrow{x}}_n=\overrightarrow{x}(t_n)$ and ${\overrightarrow{F}}_n=\overrightarrow{F}(\overrightarrow{x}(t_n),t_n)$, where $n=1,2,\dots$. The method of Euler finally follows:

${\overrightarrow{x}}_{n+1}={\overrightarrow{x}}_n+{\overrightarrow{F}}_n\tau .$

This method is illustrated in figure 8.2.

Figure 8.2.: This figure shows a comparison between the true solution (continuous line) and a numerical solution found by the Euler method in one time step (dotted).

Schematically, the algorithm is as follows:

• Choose a step width $\tau$
• Specify an initial condition ${\overrightarrow{x}}_1$
• Set the counter $n$ to $1$
• While $t_n\le t_{\text{max}}$:
  • Evaluate the right side of the DES ${\overrightarrow{F}}_n$
  • Find the new places and velocities via ${\overrightarrow{x}}_{n+1}={\overrightarrow{x}}_n+{\overrightarrow{F}}_n\tau$
  • Increase the counter by 1

The result is a trajectory ${\overrightarrow{x}}_1,{\overrightarrow{x}}_2,\dots$, which approximates the (in the general case unknown) true solution $x(t)$.

The local error is characterized as the truncation error per time step. It is of order $\mathcal{O}({\tau }^n)$, which gives $n=2$ for the method of Euler. The global error is calculated from the time interval of the length needed to calculate the trajectory $t_{\text{max}}=N\tau$. $N$ is the number of time steps. An approximation for the maximal error is

$\text{global error}\propto N\times \text{(local error)}=N\mathcal{O}({\tau }^n)=t_{\text{max}}\mathcal{O}({\tau }^{n-1}).$

The global error of the method of Euler is thus of order $\mathcal{O}(\tau )$.

Example: The Planar Pendulum

Figure 8.3.: Schematic depiction of a planar pendulum.

As an example we will discuss the planar pendulum. The pendulum consists of a point-shaped mass $m$, that is fixed on a massless rod of length $L$ (see fig. 8.2). The only degree of freedom of this system is the angle of deflection $\mathrm{\Theta }(t)\in [-\pi ,\pi ]$. The total energy is the sum of kinetic and potential energy

$E=\frac{1}{2}m({\dot{x}}^2+{\dot{y}}^2)+mgy$

with gravitational acceleration $g$. From fig. 8.2, the following relations can be extracted:

Taking the time derivation while considering $\dot{L}=0$

This yields

${\dot{x}}^2+{\dot{y}}^2=L^2{\dot{\mathrm{\Theta }}}^2(\sin {\mathrm{\Theta }}^2+\cos {\mathrm{\Theta }}^2)=L^2{\dot{\mathrm{\Theta }}}^2.$

Insertion into the energy equation finally yields

$E=\frac{1}{2}m{L}^2{\dot{\mathrm{\Theta }}}^2-mgL\cos \mathrm{\Theta }.$

This is the energy of the pendulum, expressed by the angle of deflection $\mathrm{\Theta }$ and the respective angular velocity $\dot{\mathrm{\Theta }}$. An equation of motion for the angle can be derived from the conservation of energy

$\dot{E}=m{L}^2\dot{\mathrm{\Theta }}\ddot{\mathrm{\Theta }}+mgL\dot{\mathrm{\Theta }}\sin \mathrm{\Theta }\stackrel{!}{=}0.$

This gives the equation of motion for the planar pendulum

$\ddot{\mathrm{\Theta }}=-\frac{g}{L}\sin \mathrm{\Theta }.$

In the approximation of small deflections, the problem can be solved analytically. With $\sin \mathrm{\Theta }\approx \mathrm{\Theta }$, the differential equation of the harmonic oscillator follows:

$\ddot{\mathrm{\Theta }}=-\frac{g}{L}\mathrm{\Theta }.$

We thus choose the ansatz

$\mathrm{\Theta }(t)=C_1\cos (\omega t+C_2)$

with amplitude $C_1$, phase $C_2$ and angular frequency $\omega$. Taking the second derivative of time yields

$\ddot{\mathrm{\Theta }}=-{\omega }^2\mathrm{\Theta }\Rightarrow \omega =\sqrt{\frac{g}{L}}.$

The known result is, that in the approximation of small deflections, the period of oscillation is independent of the deflection.

The general case of arbitrarily big deflections is not as easy to handle. We will solve this problem numerically. The analysis above shows, that it is sensible to choose $\sqrt{L/g}$ as the natural unit of time, i.e.

$t=\tilde{t}\sqrt{\frac{L}{g}}.$

The important point is, $\tilde{t}$ is a dimensionless quantity. This choice of scaling yields

$\frac{d^2\mathrm{\Theta }}{d{t}^2}=\frac{d^2\mathrm{\Theta }}{d{\tilde{t}}^2}\frac{g}{L}=-\frac{g}{L}\sin \mathrm{\Theta }.$

The equation of motion can thus be written with the scaled time as

$\ddot{\mathrm{\Theta }}=-\sin \mathrm{\Theta }.$

This equation of motion is free from both units and parameters! Introducing the natural unit thus not only prevents range errors, but also simplifies the problem.

With the generalized speed $\gamma =\dot{\mathrm{\Theta }}$, the equation of motion can be written as a DES of 1. order

$\dot{\mathrm{\Theta }}=\gamma$

$\dot{\gamma }=-\sin \mathrm{\Theta }.$

The following is a program that solves the equations of motion numerically via the method of Euler.

${\mathrm{\Theta }}_{n+1}={\mathrm{\Theta }}_n+\tau {\gamma }_n$

${\gamma }_{n+1}={\gamma }_n-\tau \sin {\mathrm{\Theta }}_n$

For the sake of clarity, the program is divided into functions.

import numpy as np
import matplotlib.pyplot as plt
import typing


def derive_pendulum(t: float, s: np.ndarray):
    return np.array([s[1], -np.sin(s[0])])


def euler_integration(t: float, x: np.ndarray, tau: float, derivs: typing.Callable):
    return x + tau * derivs(t, x)


# set initial conditions and constants
try:
    theta0: float = float(input("initial angle = "))
except ValueError:
    theta0 = 90
    print(f"theta0 = {theta0}°")

x: np.ndarray = np.array([theta0 * np.pi / 180, 0])

try:
    tmax: float = float(input("time = "))
except ValueError:
    tmax = 100
    print(f"tmax = {tmax}")

try:
    tau: float = float(input("step width = "))
except ValueError:
    tau = 0.001
    print(f"tau = {tau}")

time: np.ndarray = np.linspace(
    start=0,
    stop=tmax,
    num=int(np.ceil(tmax / tau)),
    endpoint=True,
)

t_plot = np.zeros_like(time)
th_plot = np.zeros_like(time)

for i, t in enumerate(time):
    # remember angle and time
    t_plot[i] = t
    th_plot[i] = x[0] * 180 / np.pi

    x = euler_integration(t, x, tau, derive_pendulum)

# figure
plt.plot(t_plot, th_plot)
plt.xlabel("time")
plt.ylabel("theta")
plt.show()

Fig. 8.4(a) shows a solution of the equation of motion by the method of Euler with step width $\tau =0.1$. We observe that the amplitude, and hence the energy, of the oscillation increases strongly with time. This violates energy conservation and has to be considered an artifact of the discretization of time. The non-physical increase of energy can be slowed down by choosing a smaller step width. This is demonstrated in fig.8.4(b) for the case of $\tau =0.05$. Unfortunately, the deviation from the true solution is still considerable, despite the small step width.

Figure 8.4.: Numerical solution of the pendulum dynamics following the method of Euler. The initial angle is ${10}^{\circ }$, and the initial speed is 0. The step width is (a) $\tau =0.1$ and (b) $\tau =0.05$. Note the different scaling of the y axes.

Runge-Kutta Method

Now we will discuss a method that is superior to that of Euler, since the error term is of higher order. The idea is to find the changes in time (derivatives) of the variable at a sampling point that lies at the half of the interval, over which the integration step is supposed to run. The actual integration step is then performed with the derivative at this point, and is one order more precise than the Euler method, as we will show in the following.

We again consider the DES

$\dot{\overrightarrow{x}}=\overrightarrow{F}(\overrightarrow{x},t)$

In the following we will compare the Taylor series

$\overrightarrow{x}(t+\tau )=\overrightarrow{x}(t)+\dot{\overrightarrow{x}}(t)\tau +\ddot{\overrightarrow{x}}(t)\frac{{\tau }^2}{2}+\mathcal{O}({\tau }^3)$

with the method of Euler with half the step width

$\overrightarrow{x}{}^{\ast }(t+\frac{\tau }{2})=\overrightarrow{x}(t)+\overrightarrow{F}(\overrightarrow{x},t)\frac{\tau }{2}$.

It can be seen that

This yields the Runge-Kutta method of 2. order

$\overrightarrow{x}(t+\tau )=\overrightarrow{x}(t)+\overrightarrow{F}(\overrightarrow{x}{}^{\ast }(t+\frac{\tau }{2}),t+\frac{\tau }{2})\tau +\mathcal{O}({\tau }^3)$

$\overrightarrow{x}{}^{\ast }(t+\frac{\tau }{2})=\overrightarrow{x}(t)+\overrightarrow{F}(\overrightarrow{x},t)\frac{\tau }{2}.$

The truncation error is of 3. order in $\tau$, i.e. one order higher as compared to the method of Euler. The calculation method is illustrated in fig.8.5.

Figure 8.5.: Scheme of the Runge-Kutta method, 2. order.

To better understand this method, we investigate a very simple differential equation

$\dot{x}=x.$

With the initial condition $x(0)=1$, the solution is

$x(t)=e^t.$

Let us now consider a translation $\tau$ in time. We make use of the rules of powers and construct the Taylor series up to the fifth order

$x(t+\tau )=e^{t+\tau }=e^t{e}^{\tau }$

$x(t+\tau )=x(t)(1+\tau +\frac{{\tau }^2}{2}+\frac{{\tau }^3}{6}+\frac{{\tau }^4}{24}+\mathcal{O}({\tau }^5)).$

Comparing this with the result obtained by the method of Euler

$x_{n+1}=x_n+\tau x_n=x_n(1+\tau )$

shows that the method of Euler gives the exact solution up to the first order. Now we look at the Runge-Kutta method of 2. order, which is

$x_{n+1}=x_n+\tau x^{\ast }$

$x^{\ast }=x_n+\frac{\tau }{2}{x}_n.$

Inserting the lower equation in the upper one yields

$x_{n+1}=x_n+\tau x_n+\frac{{\tau }^2}{2}{x}_n$

$x_{n+1}=x_n(1+\tau +\frac{{\tau }^2}{2}).$

Thus we see that the Runge-Kutta method of 2. order represents the exact solution up to the 2. order.

Runge-Kutta Method of 4. Order

In complete analogy, the Runge-Kutta method of 4. order can be derived. We will thus skip this procedure and directly look at the result:

$\overrightarrow{x}(t+\tau )=\overrightarrow{x}(t)+({\overrightarrow{F}}_1+2{\overrightarrow{F}}_2+2{\overrightarrow{F}}_3+{\overrightarrow{F}}_4)\frac{\tau }{6}+\mathcal{O}({\tau }^5)$

${\overrightarrow{F}}_1=\overrightarrow{F}(\overrightarrow{x},t)$

${\overrightarrow{F}}_2=\overrightarrow{F}(\overrightarrow{x}+\frac{\tau }{2}{\overrightarrow{F}}_1,t+\frac{\tau }{2})$

${\overrightarrow{F}}_3=\overrightarrow{F}(\overrightarrow{x}+\frac{\tau }{2}{\overrightarrow{F}}_2,t+\frac{\tau }{2})$

${\overrightarrow{F}}_4=\overrightarrow{F}(\overrightarrow{x}+\tau {\overrightarrow{F}}_3,t+\tau )$

The truncation error is of 5. order in the step width $\tau$. Let us first consider again the simple differential equation $\dot{x}=x$. The auxiliary quantities $F_i$ can be found effortlessly,

$F_1=x$

$F_2=x+\frac{\tau }{2}{F}_1=x+\frac{\tau }{2}x$

$F_3=x+\frac{\tau }{2}{F}_2=x+\frac{\tau }{2}x+\frac{{\tau }^2}{4}x$

$F_4=x+\tau F_3=x+\tau x+\frac{{\tau }^2}{2}x+\frac{{\tau }^3}{4}x.$

Inserting and rearranging yields

$x(t+\tau )=x(t)(1+\tau +\frac{{\tau }^2}{2}+\frac{{\tau }^3}{6}+\frac{{\tau }^4}{24}).$

We can see that the Runge-Kutta method of 4. order is the exact solution up to the 4. order of the series. We thus assume that this method needs far less time steps than the method of Euler. Fig. 8.6(a) confirms this. The result for the planar pendulum with a small deflection angle of ${10}^{\circ }$ solved by the Runge-Kutta method of 4. order is plotted. The case of big deflection angles (here ${170}^{\circ }$) is shown in fig. 8.6(b). From the comparison with fig.8.6(a) we can deduct that the period of oscillation depends on the deflection of the pendulum.

Figure 8.6.: Numerical solution of the dynamics of the planar pendulum via Runge-Kutta, 4. order. the step width is $\tau =0.2$. The initial angle is (a) ${10}^{\circ }$ and (b) ${170}^{\circ }$.

Adaptive Methods

Thinking about the previously introduced methods, two questions arise:

• Given a DES, how should the step width $\tau$ be chosen?
• Can the step width $\tau$ be chose dynamically, when the variables of the DES show less or more temporal variations? Such a case is given for example for the pendulum with a large deflection angle near ${180}^{\circ }$ (see the following figure). With a flexible step width, the calculation could be accelerated a lot.

Figure 8.7.: Runge-Kutta method of 4. order: pendulum with initial angle = ${179.5}^{\circ }$ and $\tau =0.2$.

Both questions lead toward the so called adaptive methods, where the step width $\tau$ is adjusted in every time time step. The procedure is the following: Find $\overrightarrow{x}(t+\tau )$ first with a “big” time step $\tau$. This gives a rough approximation ${\overrightarrow{x}}_g$. Then calculate the far better approximation ${\overrightarrow{x}}_k$ with two successive “small” time steps $\tau /2$.

With these two approximations, a relative error $\mathrm{\Delta }$ can be approximated,

with the machine precision $\epsilon$. We want that this error is approximately equal to a given value ${\mathrm{\Delta }}_0$. To achieve this, we recall that the local error of the Runge-Kutta method of 4. order is proportional to ${\tau }^5$. Thus, $\mathrm{\Delta }$ is proportional to ${\tau }^5$. The step width has to be adjusted such that ${\mathrm{\Delta }}_0\propto {\tau }_0^5$. It follows for the approximated step width:

Here, $S_1<1$ is a “security factor”, which assures that the approximated step width is not to big. A second “security factor” $S_2>1$ is introduced to prevent too large steps $\tau \to {\tau }_{new}$,

Feasible values are $S_1=0.9$ and $S_2=1.25$. By the way: As the new initial condition for the next time step ${\overrightarrow{x}}_k$ should be used, since this value is more accurate than ${\overrightarrow{x}}_g$.

Python Solver

import numpy as np
import scipy
import matplotlib.pyplot as plt


def derive_pendulum(t, s):
    return [s[1], -np.sin(s[0])]


# set initial conditions and constants
try:
    theta0: float = float(input("initial angle = "))
except ValueError:
    theta0 = 90
    print(f"theta0 = {theta0}°")

x: np.ndarray = np.array([theta0 * np.pi / 180, 0])

tmin: float = 0.0
try:
    tmax: float = float(input("time = "))
except ValueError:
    tmax = 100
    print(f"tmax = {tmax}")

try:
    tau: float = float(input("step width = "))
except ValueError:
    tau = 0.1
    print(f"tau = {tau}")

time: np.ndarray = np.linspace(
    start=tmin,
    stop=tmax,
    num=int(np.ceil((tmax - tmin) / tau)),
    endpoint=True,
)

t_eval = np.linspace(tmin, tmax, 1000)

sol = scipy.integrate.solve_ivp(
    fun=derive_pendulum, t_span=[tmin, tmax], y0=x, t_eval=time, method="RK45"
)

plt.plot(sol.t, sol.y[0] / np.pi * 180)
plt.xlabel("time")
plt.ylabel("theta")
plt.show()

scipy.integrate.solve_ivp

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options)

Solve an initial value problem for a system of ODEs.

This function numerically integrates a system of ordinary differential equations given an initial value:

Here t is a 1-D independent variable (time), y(t) is an N-D vector-valued function (state), and an N-D vector-valued function f(t, y) determines the differential equations. The goal is to find y(t) approximately satisfying the differential equations, given an initial value y(t0)=y0.

Some of the solvers support integration in the complex domain, but note that for stiff ODE solvers, the right-hand side must be complex-differentiable (satisfy Cauchy-Riemann equations). To solve a problem in the complex domain, pass y0 with a complex data type. Another option always available is to rewrite your problem for real and imaginary parts separately.

method:

The source code is Open Source and can be found on GitHub.