1.4. Higher-order ODEs

The Euler method solves a first-order ODE of the form \(y' = f(t, y)\) so what happens when we want to solve a higher order ODE? The solution is we write a higher order ODE as system of first-order ODEs and apply the Euler method to each of these equations. For example consider the \(N\)-th order ODE

\[ y^{(N)} = f(t, y, y', y'' ,\ldots ,y^{(N-1)}). \]

If we introduce functions \(y_1, y_2, \ldots, y_N\)[1] where \(y_1=y\), \(y_2 =y'\), \(y_3 =y''\) and so on then

\[\begin{split} \begin{align*} y_1' &= y_2,\\ y_2' &= y_3,\\ &\vdots \\ y_N' &= f(t, y_1 , y_2 , y_3 , \ldots, y_{N-1}). \end{align*} \end{split}\]

This is a system of \(N\) first-order ODEs so we can apply the Euler method to each equation in the system to give an equivalent solution to the \(n\)th order ODE.

Example 1.2

Rewrite the following ODE as a system of first-order ODEs

\[ y''' + yy'' -2y' + ty'' - 10 = 0. \]

Solution (click to show)

First rearrange the ODE so that the highest derivative is the subject

\[ y''' = -yy'' + 2y' - ty'' + 10, \]

then let \(y_1 = y\), \(y_2 = y'\), \(y_3 = y''\) then we can rewrite this ODE as

\[\begin{split} \begin{align*} y_1' &= y_2, \\ y_2' &= y_3, \\ y_3' &= - y_1 y_3 + 2 y_2 - t y_3 + 10. \end{align*} \end{split}\]

1.4.1. Solving systems of ordinary differential equations

We have seen that we can rewrite an \(N\)-th order ODE into a system of \(N\) first-order ODEs and we can apply the Euler method to solve each one. Since we applying the same method to multiple equations it makes sense to group them for convenience. Consider a system of \(N\) first-order ODEs

\[\begin{split} \begin{align*} y_1' &= f_1(t, y_1, \ldots, y_N), \\ y_2' &= f_2(t, y_1, \ldots, y_N), \\ & \vdots \\ y_N' &= f_N(t, y_1, \ldots, y_N). \end{align*} \end{split}\]

Let \(\vec{y}\) and \(\vec{f}(t, \vec{y})\) be the vectors

\[\begin{split} \begin{align*} \vec{y} &= \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix}, & \vec{f}(t, \vec{y}) &= \begin{pmatrix} f_1(t, \vec{y}) \\ f_2(t, \vec{y}) \\ \vdots \\ f_N(t, \vec{y}) \end{pmatrix}, \end{align*} \end{split}\]

then we can write the system in vector form

\[ \begin{align*} \vec{y}' = \vec{f}(t, \vec{y}). \end{align*} \]

The Euler method so solving a system of first-order ODEs is

\[ \begin{align*} \vec{y}_{n+1} &= \vec{y}_n + h \vec{f}(t_n, \vec{y}_n). \end{align*} \]

Where we have an initial value problem that is defined using a system of ODEs we need to know the initial value for each equation in the system, i.e., \(y_1(t_0), y_2(t_0), \ldots\). The solution is an \(n \times N\) matrix \(Y\) where the rows are formed from the solutions for \(\vec{y}_n\), i.e.,

\[\begin{split} \begin{align*} Y = \begin{pmatrix} \leftarrow & \textbf{y}_0 & \rightarrow \\ \leftarrow & \textbf{y}_1 & \rightarrow \\ & \vdots \\ \leftarrow & \textbf{y}_n & \rightarrow \end{pmatrix} \end{align*} \end{split}\]

Example 1.3

The van der Pol equation is a model of a non-conservative oscillator with non-linear damping which is

\[ \ddot{y} - \mu (1 - y^2) \dot{y} + y = 0, \]

where \(y\) is the position of the oscillator, \(\mu\) a parameter for the strength of the dampening and \(\dot{y}\) and \(\ddot{y}\) are the first and second-order derivatives with respect to time \(t\).

Use the Euler method to calculate the solution to the van der Pol equation with \(\mu=1\) over the domain \(t \in [0, 20]\) with initial conditions \(y(0)=2\) and \(\dot{y}(0)=0\) using a step length of \(h=0.1\).

Solution (click to show)

First we need to write the ODE as a system of two first-order ODEs. Let \(y_1=y\) and \(y_2=\dot{y}\) then

\[\begin{split}\begin{align*} \dot{y}_1 &= y_2, \\ \dot{y}_2 &= \mu (1 - y_1^2) y_2 - y_1, \end{align*} \end{split}\]

which can be written in vector form as

\[\begin{split}\begin{align*} \vec{y} &= \begin{pmatrix} y_1 \\ y_2 \end{pmatrix}, & \vec{f}(t, \vec{y}) &= \begin{pmatrix} y_2 \\ \mu (1 - y_1^2) y_2 - y_1 \end{pmatrix}. \end{align*} \end{split}\]

The initial conditions are \(y(0) = 2\) and \(\dot{y}(0) = 0\) so \(\vec{y}_0 = \begin{pmatrix} 2 , 0 \end{pmatrix}\). The values of \(t\) are

\[\begin{split} \vec{t} = \begin{pmatrix} 0.0 \\ 0.1 \\ 0.2 \\ \vdots \\ 19.9 \\ 20.0 \end{pmatrix}.\end{split}\]

Calculating the first few steps of the Euler method with \(\mu=1\) and \(h=0.1\)

\[\begin{split} \begin{align*} \vec{y}_1 &= \vec{y}_0 + h \vec{f}(t_0, \vec{y}_0) \\ &= \begin{pmatrix} 2 \\ 0 \end{pmatrix} + 0.1 \begin{pmatrix} 0 \\ 1(1 - 2^2)(0) - 2 \end{pmatrix} = \begin{pmatrix} 2 \\ -0.2 \end{pmatrix}, \\ \vec{y}_2 &= \vec{y}_1 + h \vec{f}(t_1, \vec{y}_1) \\ &= \begin{pmatrix} 2 \\ -0.2 \end{pmatrix} + 0.1 \begin{pmatrix} -0.2 \\ 1(1 - 2^2)(-0.2) - 2 \end{pmatrix} = \begin{pmatrix} 1.98 \\ -0.34 \end{pmatrix}, \\ \vec{y}_3 &= \vec{y}_2 + h \vec{f}(t_2, \vec{y}_2) \\ &= \begin{pmatrix} 1.98 \\ -0.34 \end{pmatrix} + 0.1 \begin{pmatrix} -0.34 \\ 1(1 - 1.98^2)(-0.34) - 1.98 \end{pmatrix} = \begin{pmatrix} 1.946 \\ -0.348706 \end{pmatrix}, \\ &\vdots \end{align*} \end{split}\]

Continuing for the whole domain we get

\[\begin{split} Y = \begin{pmatrix} 2 & -0.2 \\ 1.98 & -0.34 \\ \vdots & \vdots \\ -0.828052 & 1.314658 \\ -0.696586 & 1.438787 \end{pmatrix}. \end{split}\]

The first column of \(Y\) is the solution \(y\) for the van der Pol equation which has been plotted in Fig. 1.9

Fig. 1.9 Solution to the van der Pol equation with \(y(0)=2\) and \(\mu=1\) using the Euler method.

1.4.2. Code

The code used to compute the solutions to Example 1.3 is shown below (this assumes the functions solveIVP() and euler() have already been defined). Note that the y0 array and the values returned by the f() function are \(1 \times 2\) arrays where the number of columns is the same as the number of equations in the system of ODEs.

Python

# Define ODE function (van der Pol equation)
def f(t,y):
    return np.array([ y[1], mu* (1 - y[0] ** 2) * y[1] - y[0] ])

# Define IVP parameters
tspan = [0, 20]  # boundaries of the t domain
y0 = [2, 0]      # initial values
h = 0.1          # step length
mu = 1           # dampening parameter

# Solve the IVP using the Euler method
t, y = solveIVP(f, tspan, y0, h, euler)

# Plot solution
fig, ax = plt.subplots()
plt.plot(t, y[:,0], "b-")
plt.xlabel("$t$", fontsize=12)
plt.ylabel("$y$", fontsize=12)
plt.show()

MATLAB

% Define ODE function (van der Pol equation)
f = @(t, y, mu) [ y(2), mu * (1 - y(1) ^ 2) * y(2) - y(1) ];

% Define IVP parameters
tspan = [0, 20];    % boundaries of the t domain
y0 = [2, 0];        % initial values
h = 0.1;            % step length
mu = 1;             % dampening parameter

% Solve the IVP using the Euler method
[t, y] = solveIVP(@(t,y)f(t, y, mu), tspan, y0, h, @euler);

% Plot solution
plot(t, y(:, 1), 'b', LineWidth=2)
xlabel('$t$', Fontsize=14, Interpreter='latex')
ylabel('$y$', Fontsize=14, Interpreter='latex')
axis padded

---

[1]

Here \(y_1, y_2, \ldots\) are functions and should not be confused for the solution values \(y_1=y(t_1), y_2=y(t_2), \ldots\).