1.4. Higher-order ODEs

The Euler method and the numerical methods developed later in this module are formulated for first-order initial value problems of the form

\[ y' = f(t, y). \]

However, many important mathematical methods involve second-order or higher-order differential equations. Fortunately, any higher-order ODE can be rewritten as a system of first-order ODEs, allowing the same numerical methods to be applied.

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\) 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-1}' &= y_N, \\ 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 solve the system to give an equivalent solution to the \(N\)-th order ODE.

Example 1.3

A mass-spring-damper model consists of objects connected via springs and dampers. A simple example of the application of a model is a single object connected to a surface is shown in Fig. 1.8.

Fig. 1.8 The mass-spring-damper model [Wikipedia, 2008].

The displacement of the object over time can be modelled by the second-order ODE

\[ m \ddot{y} + c \dot{y} + k y = 0. \]

The three terms in the equation represent

• inertial forces: \(m\ddot{y}\)

• damping forces: \(c \dot{y}\)

• spring restoring forces: \(ky\)

Together these determine how the object oscillates and gradually returns to equilibrium.

An object of mass 1 kg is connected to a dampened spring with \(c = 2\) and \(k = 4\). The object is displaced by 1m and then released. Use the Euler method to compute the displacement of the object over the first 5 seconds after it was released.

---

Solution

Rearranging the governing equation to make \(\ddot{y}\) the subject gives

\[ \ddot{y} = \frac{1}{m}(-c \dot{y} - k y). \]

We need to rewrite this as a system of two first-order ODEs. Let \(y_1 = y\) and \(y_2 = \dot{y}\) then we have the system of two first-order ODEs

\[\begin{split} \begin{align*} \dot{y}_1 &= y_2, \\ \dot{y}_2 &= \frac{1}{m}(- c y_2 - k y_1). \end{align*} \end{split}\]

This is now a system of two first-order ODEs of exactly the form considered in the previous section. Since \(m = 1\), \(c = 2\) and \(k = 4\) then writing the system in vector form \(\dot{\mathbf{y}} = f(t, \mathbf{y})\) gives

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

The initial conditions are \(y_1 = 1\) (displacement) and \(y_2 = 0\) (velocity of the object). Using the Euler method with a step length of \(h = 0.1\)

\[\begin{split} \begin{align*} \mathbf{y}_1 &= \mathbf{y}_0 + h \mathbf{f}(t_0, \mathbf{y}_0), \\ &= \begin{pmatrix} 1 \\ 0 \end{pmatrix} + 0.1 \begin{pmatrix} 0 \\ -2(0) - 4(1) \end{pmatrix} = \begin{pmatrix} 1 \\ -0.4 \end{pmatrix}, \\ t_1 &= t_0 + h = 0 + 0.1 = 0.1, \\ \\ \mathbf{y}_2 &= \mathbf{y}_1 + h \mathbf{f}(t_1, \mathbf{y}_1), \\ &= \begin{pmatrix} 1 \\ -0.4 \end{pmatrix} + 0.1 \begin{pmatrix} -0.4 \\ -2(-0.4) - 4(1) \end{pmatrix} = \begin{pmatrix} 0.96 \\ -0.72 \end{pmatrix}, \\ t_2 &= t_1 + h = 0.1 + 0.1 = 0.2, \\ \\ \mathbf{y}_3 &= \mathbf{y}_2 + h \mathbf{f}(t_2, \mathbf{y}_2), \\ &= \begin{pmatrix} 0.96 \\ -0.72 \end{pmatrix} + 0.1 \begin{pmatrix} -0.72 \\ -2(-0.72) - 4(0.96) \end{pmatrix} = \begin{pmatrix} 0.888 \\ -0.96 \end{pmatrix}, \\ t_3 &= t_2 + h = 0.2 + 0.1 = 0.3, \\ & \vdots \end{align*} \end{split}\]

Continuing this process produces approximations at equally spaced time points throughout the interval \(0 \le t \le 5\). Fig. 1.9 shows that the displacement oscillates about the equilibrium position \(y = 0\), while the amplitude gradually decreases due to the damping term \(c\dot{y}\).

Fig. 1.9 Euler approximation of the displacement of the mass over the first 4 seconds.

1.4.1. Code

The code below sets up and solves the mass-spring-damper model example from Example 1.3 using the Euler method.

Python

# Define Mass-Spring model
def spring(t, y):
    return np.array([ y[1], (-c*y[1] - k*y[0]) / mass ])
    

# Define IVP parameters
tspan = [0, 5]  # boundaries of the t domain
y0 = [1, 0]     # initial values
h = 0.1         # step length
mass = 1;       # mass of object
c = 2;          # damping coefficient
k = 4;          # spring constant

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

# Plot solution
fig, ax = plt.subplots()
plt.plot(t, y[:,0], "b-")
plt.xlabel("time (s)")
plt.ylabel("displacement (m)")
plt.show()

MATLAB

% Define mass-spring-damper model
spring = @(t, y, m, c, k) [ y(2) ; (-c*y(2) - k*y(1)) / m];

% Define IVP parameters
tspan = [0, 5];  % boundaries of the t domain
y0 = [1, 0];     % initial values [S, I, R]
h = 0.1;         % step length
m = 1;           % mass of object
c = 2;           % damping coefficient
k = 4;           % spring constant

% Solve the IVP using the Euler method
[t, y] = euler(@(t, y)spring(t, y, m, c, k), tspan, y0, h);

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