Numerical solution of differential equations - initial value problem

Introduction¶

Writing (a single) differential equation¶

Assume that a first-order differential equation can be written in explicit form:

where $f(t, y)$ is a given function and $y'=dy/dt$.

Additionally, an initial condition is given:

The goal of solving the differential equation is to compute the function $y(t)$ that solves the initial value problem defined above. Under certain conditions on the function $f(t, y)$, the initial value problem has a unique solution on an interval containing $t_0$.

In numerical solving we always compute a table of the function $y(t_i)$ that solves the given initial value problem. The nodes $t_i$ are usually equidistant:

and we call $h$ the (time) step (of integration).

Here we will look at some numerical methods for solving differential equations with an initial condition.

Euler’s method¶

The explicit Euler’s method is based on the Taylor series expansion of the function $y$:

We make a method error of $\mathcal{O}(h^2)$ because we neglect the second and higher-order derivatives; now, given a known value $y(t)$ and the derivative $y'(t)=f(t,y)$, we can estimate the value at the next time step $t+h$. Once we know the values at $t+h$, we repeat the procedure!

Steps of Euler’s method:

1. Set $i=0$, $t_0$, $y_0=y(t_0)$.

2. Compute the function value at $t_{i+1}=t_i+h$: $y_{i+1}= y_i + f(t_i, y_i)\,h.$

3. $i=i+1$ and continue at step 2.

We solve the differential equation on the interval $[t_0,t_n]$ and $h=(t_n-t_0)/n$. $n$ is the number of integration steps (how many times we perform step 2 in the algorithm above).

Numerical solution of the initial value problem:

at the values of the independent variable:

Error of Euler’s method¶

The error of Euler’s method at each step is of order $\mathcal{O}(h^2)$.

Since over the interval from $t_0$ to $t_n$ we make this error $n$ times, the cumulative error is $n\,\mathcal{O}(h^2)=\frac{t_n-t_0}{h}\,\mathcal{O}(h^2)=\mathcal{O}(h)$.

Locally the error is second order, but globally the error is first order, and because Euler’s method is so inaccurate, we rarely use it in practice!

Error estimate¶

The exact solution $y(t_n)$ at step size $h$ is:

where $y_{n,h}$ is the numerical approximation and $E_h$ is the method error. Since the global error is first order, we can write the error as:

Similarly, for step size $2h$ we can write:

where $y_{n,2h}$ is the numerical approximation and $E_{2h}$ is the method error:

Assuming that the constant $k$ is the same for step $h$ and step $2h$, we can determine the error estimate for the better approximation $E_h$. Clearly:

then we determine the error estimate:

A note on the implicit Euler’s method¶

In the explicit Euler’s method we compute the solution at $t_{i+1}$ from the computed value at $t_i$.

If the unknown value of the solution at $t_{i+1}$, that is $y_{i+1}$, appeared on the right-hand side as well, we would be talking about the implicit Euler’s method (or backward Euler’s method):

Because the sought value $y_{i+1}$ appears on both sides of the equation, to determine $y_{i+1}$ we must solve a (nonlinear) equation. The advantage of the implicit Euler’s method is that it is more stable (e.g., for stiff systems, which we will learn about later) than the explicit form, but it is numerically more demanding (because of solving the equation).

Numerical implementation¶

First we import the required libraries:

Next let us define Euler’s method:

Let us prepare a function for the error estimate (in numerical terms it would be better to include the error estimate inside the euler function, but for clarity we present it separately):

Numerical example¶

As an example, let us solve the differential equation that describes a falling body subjected to gravity and air resistance:

According to Newton’s second law, we can write the differential equation:

where $m$ is the mass, $g$ the gravitational acceleration, $c$ the air-resistance coefficient, and $v$ the velocity. We would like to solve the differential equation subject to the initial condition:

The right-hand side / first-derivative function $f(t,y)$ is:

and the initial condition:

Let us define the right-hand-side function:

Let us define the initial condition and the time vector where we want the result:

We call the euler function to compute the values of $y$ (velocity $v$):

We plot the result:

Let us now check the effect of the time step:

We observe that the numerical error decreases substantially when the step changes from 1 to 0.1!

Let us also estimate the error at 100 and 1000 subintervals:

When we reduced the step to a tenth, the error decreased proportionally as well (first-order error).

Let us also look at the case where the air resistance c is an argument of the euler function and is passed via **kwargs to the f_air_resistance() function:

Second-order Runge–Kutta method¶

Euler’s method is first order (the global error $\mathcal{O}(h)$ is indeed first order). If we wanted to derive a second-order-error method, we would need to use the Taylor series expansion of $y(t+h)$, where we neglect the third and higher derivatives:

The local error of the method will thus be third order and the global error second order.

We use the substitutions $y'(t)=f(t,y)$ and $y''(t)=f'(t,y)$:

Because the right-hand side $f(t,y)$ depends on both the independent variable $t$ and the dependent variable $y$, we must use implicit differentiation:

We substitute into the Taylor series expression:

As can be seen from the expression above, we need additional derivatives. This poses a certain difficulty that can be avoided in various ways; in what follows we will look at the Runge–Kutta approach. Because we will need the expression above later, we have named it here $y(t+h)_{\textrm{Taylor}}$.

Idea of the Runge–Kutta approach¶

The dilemma above is resolved by the Runge–Kutta method (developed in 1901) with an idea we have already encountered in Gaussian integration: it tries to find a more accurate solution with a weighted additional value of the function $f$:

where $c_0$, $c_1$, $p$, and $q$ are unknown constants (in principle from 0 up to and including 1). If we used $c_1=0$ in the expression above, we would derive a first-order method; but with an additional function value ($A$) it will turn out that we derive a second-order method.

We continue the search for the unknown constants $c_0$, $c_1$, $p$, $q$ by writing $A$ in the form of a first-order Taylor series:

Now we substitute the derived $B$ back into the expression for $y(t+h)_{\textrm{Runge-Kutta}}$:

We continue the derivation and rearrange the equation so that it resembles the one derived above with the Taylor series $y(t+h)_{\textrm{Taylor}}$:

Now let us compare with the expression derived above:

We find that for equality to hold we must have:

So we have three equations and four unknowns. We can therefore choose one of the constants arbitrarily and compute the other three. If, for example, we choose $c_0=0$, we would call this the modified Euler’s method and the remaining parameters would be: $c_1=1$, $p=q=1/2$. The choice of parameters has no essential effect on the solution. It is worth noting, however, that even the second-order Runge–Kutta method is rarely used, since better methods exist.

We substitute the parameters $c_0$, $c_1$, $p$, and $q$ into the first equation of this section. Once the initial time $t_0$ and the initial condition $y_0$ are defined, we use the second-order Runge–Kutta method:

Fourth-order Runge–Kutta method¶

Similarly to how we derived the second-order Runge–Kutta method, the fourth-order Runge–Kutta method is derived. For the fourth-order method, too, there are several variants, and by the fourth-order Runge–Kutta method we mean the following method:

where:

The steps of the fourth-order Runge–Kutta method are:

1. Set $i=0$ and $t_0$, $y_0=y(t_0)$,

2. Compute the coefficients: $k_0$, $k_1$, $k_2$, $k_3$,

3. Compute the value of the solution of the differential equation at $t_{i+1}=t_i+h$: $\quad y_{i+1}=y_i+\frac{1}{6}(k_0+2\,k_1+2\,k_2+k_3),$

4. $i=i+1$ and continue at step 2.

Error of the fourth-order Runge–Kutta method¶

The fourth-order Runge–Kutta method is named so because it has a local error of fifth order $\mathcal{O}(h^5)$, but we make this error $n$ times, so the global error is fourth order $\mathcal{O}(h^4)$.

Error estimate¶

The exact result $y(t_n)$ at step size $h$ is:

where $y_{n,h}$ is the numerical approximation of the solution and $E_h$ is the method error. Since the global error is fourth order, we can write the error as:

Similarly, for step size $2h$ we can write:

where $y_{n,2h}$ is the numerical approximation of the solution and $E_{2h}$ is the method error:

Assuming that the constant $k$ is the same for step $h$ and step $2h$, we can compute the error estimate for the better approximation $E_h$.

First, it is true that:

it follows:

and then we determine the error estimate of the more accurate solution:

Numerical implementation¶

Function for the error estimate:

Numerical example¶

Let us now look at the example of computing the velocity of a falling mass:

Let us provide the initial condition and the time vector where we want the result:

For comparison, let us compute the solution with the euler and runge_kutta_4 functions:

We plot the result:

Let us also look at the numerical error:

When the step is reduced to a tenth, the error decreased by about 10⁴ times (which matches the expectation for a fourth-order method).

Using scipy to solve ordinary differential equations¶

The scipy package has many numerical methods implemented for solving initial value problems of ordinary differential equations. Here we will look at the function scipy.integrate.solve_ivp, which is suitable for most initial value problems.

scipy.integrate.solve_ivp¶

Syntax for use (IVP - Initial Value Problem):

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

An explanation of all the arguments is in the documentation; here we will highlight a few:

• fun is the right-hand side (func(t, y...)),

• t_span a tuple (t0, tf) defining the initial time t0 and final time tf,

• y0 a list of initial value(s),

• method defines the numerical method (the default is the explicit ‘RK45’; for stiff systems the following are better: ‘Radau’, ‘BDF’, ‘LSODA’),

• t_eval is used if we want the solution at specific values of the independent variable,

• dense_output whether a continuous solution is also prepared, default False,

• events for tracking events that stop solving the initial value problem (e.g., when the mass touches the ground and the relative distance to the ground is zero).

• args for passing arguments into the function fun (see the example below).

The result of calling solve_ivp is an object with attributes (selected):

• t the values of the independent variable at which the result is computed,

• y the result,

• sol the continuous-solution function (only in the case of dense_output=True)

• success is True if the solution was computed successfully or terminated due to an event (events).

Numerical example¶

Let us define the initial condition and the time vector where we want the result (we show only the first ten elements):

We pass the arguments (g=9.81, m=1., c=0.4) into the function via the args parameter:

We could also use a lambda expression (documentation):

y_ivp = solve_ivp(lambda t,y: f_air_resistance(t, y, c=0.4), t_span=(t[0], t[-1]), y0=[v0])

The lambda expression used above, lambda t, y: f_air_resistance(t, y, c=0.4), is equivalent to:

def function_name(t, y):
    return f_air_resistance(t, y, c=0.4)

Let us look at the result:

Because solve_ivp is already set up for a system of first-order ordinary differential equations, the result is given in two dimensions. The sought result is therefore: y_ivp.y[0].

For comparison, here is also the result with our own fourth-order Runge–Kutta implementation:

We plot the result:

System of ordinary differential equations¶

Above we looked at solving the initial value problem of a single ordinary differential equation; now we will generalize the solving to an initial value problem for a system of $m$ first-order ordinary differential equations.

Such a system is written as:

where additional (initial) conditions are given:

By $t$ we denote the independent variable (it need not always be time) and by $\mathbf{f}$ the vector of right-hand sides.

We compute the solution of a system of $m$ differential equations $\mathbf{y}'=\mathbf{f}(t, \mathbf{y})$ at $a=t_0, t_1,\dots,t_{n-1}=b$, that is $m$ function values at each step $y_k(t_i)$, $k=0,1,\dots,m-1$ and $i=0,1,\dots,n-1$. It can be proved that any method we derived for first-order differential equations can also be used for a system of first-order ordinary differential equations if we replace the scalar quantities with the corresponding vector ones.

Numerical implementation¶

The Euler and fourth-order Runge–Kutta methods need only minor adjustments (y0 is a numerical array, f returns a list of derivative values):

Numerical example¶

We extend the falling mass to planar motion; let the magnitude of the air-resistance force be defined as:

The air-resistance force acts in the opposite direction to the velocity vector.

With the help of the figure, we define the force in the $x$ direction:

Similarly, the force in the $y$ direction is:

According to Newton’s second law, we write a system of two (coupled) first-order differential equations:

$m$ is the mass, $g$ the gravitational acceleration, $c$ the air-resistance coefficient, and $v_x$ and $v_y$ the velocities in the $x$ and $y$ directions, respectively. We would like to solve the differential equation subject to the initial condition:

Let us define the list of right-hand sides:

Let us define the initial condition and the time vector where we want the result:

For comparison, let us compute the solution with the runge_kutta_4 and solve_ivp functions:

Let us look at the result:

We plot the result:

Reformulating a higher-order differential equation into a system of first-order differential equations¶

We will look at how an ordinary differential equation of arbitrary order $n$:

with initial conditions:

is reformulated into a system of first-order differential equations.

First, instead of derivatives we introduce new variables:

The $n$th-order differential equation written with the new variables is:

We differentiate the new functions with respect to the independent variable:

and

We have obtained a system of first-order ordinary differential equations:

with initial conditions:

Numerical example¶

We return to the falling mass:

But this time let us write Newton’s second law in terms of displacement (a second-order differential equation):

$m$ is the mass, $g$ the gravitational acceleration, $c$ the air-resistance coefficient, and $x''$ and $y''$ the acceleration in the chosen coordinate system. Initial conditions:

We have a system of two second-order differential equations. By introducing new variables $y_i$:

We prepare the system of first-order differential equations:

Let us define the Python function of right-hand sides / first derivatives:

Let us define the initial condition and the time vector where we want the result:

We compute the solution:

Let us look at the result ($[x, x', y, y']$):

We plot the results; first the velocity, then the position (the $y$ coordinate is positive downward)!

Stability of solving differential equations*¶

The numerical solution of differential equations is subject to method error and round-off error. These errors differ for different numerical methods.

Solving a differential equation is stable if a small change in the initial condition leads to a small change in the computed solution; otherwise we speak of instability of the solution.

Stability depends on the differential equation, on the numerical method used, and on the integration step $h$.

Example of a simple oscillator¶

Let us first look at an example of solving the differential equation of a simple oscillator.

The figure (source: Slavič, Dinamika, mehanska nihanja in mehanika tekočin, 2017) shows a dynamic system (mass $m$, stiffness $k$) whose differential equation is

We reformulate this differential equation into the standard form of free vibration:

where the natural angular frequency is:

and we expect a response of the form:

If the initial conditions are:

the solution of the initial value problem is:

Numerical example¶

First let us define the vector of initial conditions and the function of right-hand sides / first derivatives (we convert the second-order differential equation into a system of first-order differential equations $\mathbf{y}'=\mathbf{f}(t, \mathbf{y})$), where $y_0=x, y_1=x'$:

Let us define the data and the analytical solution:

Solution using the Euler and fourth-order Runge–Kutta methods:

Let us plot the results:

We observe that Euler’s method is unstable, and if we increased the step, the fourth-order Runge–Kutta method would also become unstable.

Why is Euler’s method so unstable?¶

Let us recall Euler’s method:

which tells us that we determine the displacement $x(t+h)$ from the position $x(t)$ and the velocity $x'(t)$. The initial conditions stem from the maximum deflection $x(t=0)=1$, and at that point the velocity is $x'(t=0)=0$, which means that $x(t+h)=1$. Already in the first step we therefore make a relatively large error. But why does the value then start to grow in an alternating fashion?

Recall that the analytical solution is $x(t)=x_0\,\cos(\omega_0\,t)$, and therefore $x'(t)=-\omega_0\,x_0\,\sin(\omega_0\,t)$.

We substitute the prepared expressions into Euler’s method and rearrange:

Assume we look at the state at a time $t=\pi/(2\omega_0)$, when $\cos(\omega_0\,t)=0$ and $\sin(\omega_0\,t)=1$:

If the absolute value of the expression $A$ is greater than 1, then at time $t+h$ the value will be greater than in the previous step, and probably again in the next one. It follows that instability can occur. To avoid it, the following must hold:

Note: in some books you will also see the value $h<2/\omega_0$; a unique bound for all differential equations cannot be defined; in general, however, the step is defined relative to the shortest period $T$ in the differential equation (e.g.: $h<2/\omega_0$ is essentially the same as $h<2/(2\pi/T)$, i.e., $h<T/\pi$). The period $T$ is defined relative to the highest natural frequency of the system $T=1/f_{\textrm{max}}$, which we compute from the eigenvalue of the system.

Additional: Example of the Van der Pol oscillator¶

The purpose of this example is to show how the choice of integrator can affect the speed of solving a problem! The Van der Pol oscillator is described here.

Let us define the list of derivatives:

Solution using the RK45 method (this is an explicit scheme, which is not suitable for stiff systems of differential equations; solving is very slow, so we compute the solution only up to t1/100):

The implicit BDF scheme (Backward Differentiation Formulas) turns out here to be much more suitable. Because of its stability, the steps can be much larger and therefore the solving is much faster:

Additional: symbolic solving of a second-order differential equation¶

We will look at the example shown in the figure, where a mass $m$ is on an incline of slope $\alpha$. The coefficient of friction is $\mu$, and the gravitational acceleration is $g$. The initial velocity is $\dot y(0\,\textrm{s})=v_0$, the displacement $y(0\,\textrm{s})=0\,\textrm{m}$. The equation of motion (only for the $y$ direction) is defined according to Newton’s second law (see the free-body force diagram).

Derivation of the equation of motion¶

The solution of the equation is:

To determine $C_1$ and $C_2$, we substitute $t=0\,\textrm{s}$:

Then we differentiate with respect to time and again substitute $t=0\,\textrm{s}$:

Based on the initial conditions, we have therefore determined the constants:

The solution follows:

Let us prepare functions for the numerical call:

Let us prepare the plot:

Symbolic reformulation of the differential equation into a system of first-order differential equations¶

Let us recall the original differential equation:

Let us define new variables and prepare the function $f$:

Now let us connect the new variables.

Let $d y_0/dt$ be equal to $y_1$:

Let the derivative $d y_1/dt$ (which is essentially $y''$) be equal to the function $f$:

We write the expressions above in vector form:

Let us now recall f_vec:

If we solve numerically, then we write the function $\mathbf{f}(t, \mathbf{y})$ as:

Let us check the function at the initial time $t=0\,$s and at the initial conditions $[y_0, y_1]=[0, v_0]$:

Now let us use Euler’s method:

We plot and compare with the analytical solution: