6.4. Systems of ODEs and Higher Order ODEs

References:

• Section 6.3 Systems of Ordinary Differential Equations in [Sauer, 2022], to Sub-section 6.3.1 Higher order equations.

• Section 5.9 Higher Order Equations and Systems of Differential Equations in [Burden et al., 2016].

The short version of this section is that the numerical methods and algorithms developed so for for the initial value problem

\[\begin{split} \begin{split} \frac{d u}{d t} &= f(t, u(t)), \quad a \leq t \leq b \\ u(a) &= u_0 \end{split} \end{split}\]

all also work for system of first order ODEs by simply letting \(u\) and \(f\) be vector-valued, and for that, the Python code requires only one small change.

Also, higher order ODE’s (and systems of them) can be converted into systems of first order ODEs.

6.4.1. Converting a second order ODE to a first order system

To convert

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

with initial conditions

\[ y(a) = y_0, \; y'(a) = v_0 \]

to a first order system, introduction the two functions

\[\begin{split} \begin{split} u_1(t) &= y(t) \\ u_2(t) &= \frac{d y}{d t}, = u_1'(t) \end{split} \end{split}\]

Then

\[ y'' = u_1' = f(t, u_0, u_1) \]

and combining with the definition of \(u_1\) gives the system

\[\begin{split} \begin{split} u_0' &= u_1 \\ u_1' &= f(t, u_0, u_1) \\ &\text{with initial conditions} \\ u_0(a) &= y_0 \\ u_1(a) &= v_0 \end{split} \end{split}\]

Next this can be put into vector form. Defining the vector-valued functions

\[\begin{split} \begin{split} \tilde{u}(t) &= \langle u_1(t), u_2(t) \rangle \\ \tilde{f}(t, \tilde{u}(t)) &= \left\langle u_1(t), f(t, u_2(t), u_2(t)) \right\rangle \end{split} \end{split}\]

and initial data vector

\[\tilde{u}_0 = \langle u_{0,1}, u_{0,2} \rangle = \langle y_0, v_0 \rangle\]

puts the equation into the form

\[\begin{split} \begin{split} \frac{d \tilde{u}}{d t} &= \tilde{f}(t, \tilde{u}(t)), \quad a \leq t \leq b \\ \tilde{u}(a) &= \tilde{u}_0 \end{split} \end{split}\]

\[\begin{split} \begin{split} \frac{d \tilde{u}}{d t} &= \tilde{f}(t, \tilde{u}(t)), \quad a \leq t \leq b \\ \tilde{u}(a) &= \tilde{u}_0 \end{split} \end{split}\]

6.4.2. Test Cases

In this and subsequent sections, numerical methods for higher order equations and systems will be compared using several test cases:

Test Case A: Motion of a (Damped) Mass-Spring System in One Dimension

A simple mathematical model of a damped mass-spring system is

(6.1)\[\begin{split} \begin{split} M \frac{d^2 y}{d t^2} &= -K y - D \frac{d y}{d t} \\ & \text{with initial conditions} \\ y(a) &= y_0 \\ \left. \frac{dy}{dt} \right|_{t=a} &= v_0 \end{split} \end{split}\]

where \(K\) is the spring constant and \(D\) is the coefficient of friction, or drag.

The first order system form can be left in terms of \(y\) and \(y'\) as

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} y \\ y' \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -K & -D \end{array}\right] \left[\begin{array}{c} y \\ y' \end{array}\right] \end{split}\]

Exact solutions

For testing of numerical methods in this and subsequent sections, here are the exact solutions.

They depend on whether

• \(D < D_0 := 2 \sqrt{K M}\): underdamped,

• \(D > D_0\) : overdamped, or

• \(D = D_0\) : critically damped.

We will mostly explore the first two more “generic” cases.

For the underdamped case, the general solution is

\[ y(t) = e^{-(D/(2M)) (t-a)} [ A \cos(\omega(t-a)) + B \sin (\omega(t-a))], \quad \omega = \frac{\sqrt{4KM - D^2}}{2M} \]

For the above initial conditions, \(A = y_0\) and \(B = (v_0 + y_0 D/(2M)/\omega\).

An important special case of this is the undamped system \(M \frac{d^2 y}{d t^2} = -K y\) for which the solutions become

\[ y(t) = A \cos(\omega(t-a)) + B \sin (\omega(t-a)), \quad \omega =\sqrt{K/M} \]

and it can be verified that the “energy”

\[ E(t) = \frac{M}{2}(y'(t))^2 + \frac{K}{2}(y(t))^2 = \frac{1}{2}(K u_1^2 + {M}u_2^2) \]

is conserved: \(dE/dt = 0\). Conserved quantities can provide a useful check of the acccuracy of numerical method, so we will look at this below.

For the overdamped case, the general solution is

\[ y(t) = A e^{\lambda_+ (t-a)} + B e^{\lambda_- (t-a)}, \quad \lambda_\pm = \frac{-D \pm \Delta}{2M},\; \Delta = \sqrt{D^2 - 4KM} \]

For the above initial conditions, \(A = M(v_0 - \lambda_- y_0)/\Delta\) and \(B = y_0 - A\).

Remark 6.2 (Stiffness)

Fixing \(M\) and scaling \(K = D \to \infty\), \(\Delta = D \sqrt{1 - 4M/D} \approx D - 2M\) so

\[ \lambda_- \approx -\frac{D}{M} + 1 \to -\infty, \quad \lambda_+ \approx -1. \]

Thus the time scales of the two exponential decays become hugely different, with the fast term \(B e^{\lambda_- (t-a)}\) becoming negligible compared to the slower decaying \(A e^{\lambda_+ (t-a)}\).

This is a simple example of stiffness, and influences the choice of a good numerical method for solving such equations.

The variable can be rescaled to the case \(K = M = 1\), so that will be done from now on, but of course you can easily experiment with other parameter values by editing copies of the Jupyter notebooks.

Test Case B: A “Fast-Slow” Equation

The equation

\[ y'' + (K+1) y' + K y = 0, \quad y(0) = y_0, y'(0) = v_0 \]

has first order system form

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} y \\ y' \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ -K & -(K+1) \end{array}\right] \left[\begin{array}{c} y \\ y' \end{array}\right] \end{split}\]

and the general solution

\[ y(t) = A e^{-t} + B e^{-Kt} \]

so for large \(K\), it has two very disparate time scales, with only the slower scale of much significance after an initial transient.

This is a convenient “toy” example for testing two refinements to algorithms:

• Variable time step sizes, so that they can be short during the initial transient and longer later, when only the \(e^{-t}\) behavior is significant.

• Implicit methods that can effectively suppress the fast but extremely small \(e^{-kt}\) terms while hanling the larger, slower terms accurately.

The examples below will use \(K=100\), but as usual, you are free to experiment with other values.

Test Case C: The Freely Rotating Pendulum

Both the above equations are constant coefficient linear, which is convenient for the sake of having exact solution to compare with, but one famous nonlinear example is worth exporing too.

A pendulum with mass \(m\) concentrated at a distnace \(L\) from the axis of rotation and that can rotate freely in a vertical plane about that axis and with possible friction proportional to \(D\), can be modeled in terms of its angular position \(\theta\) and angular velocity \(\omega = \theta'\) by

\[ M L \theta'' = -M g\sin\theta - D L \theta', \quad \theta(0) = \theta_0, \theta'(0) = \omega_0 \]

or in system form

\[\begin{split} \frac{d}{dt} \left[\begin{array}{c} \theta \\ \omega \end{array}\right] = \left[\begin{array}{c}\omega \\ -\frac{g}{L}\sin\theta -\frac{D}{M} \omega \end{array}\right] \end{split}\]

These notes will mostly look at the frictionelss case \(D=0\), which has conserved energy

\[ E(\theta, \omega) = \frac{ML}{2} \omega^2 - M g\cos\theta \]

For this, the solution fall into three qualitatively different cases depending on whether the energy is less than, equal to, or greater than the “critical energy” \(Mg\), which is the energy of the unstable stationary solutions \(\theta(t) = \pi (\mod 2\pi)\), \(\omega(t) = 0\): “balancing at the top”:

• For \(E < Mg\), a solution can never reach the top, so the pendulum rocks back and forth, reach maximum height at \(\theta = \pm \arccos(-E/(Mg))\)

• For \(E > Mg\), solutions have angular speed \(|\omega| \geq \sqrt{E -Mg} > 0\) so it never drops to zero, and so the direction of rotation can never reverse: solutions rotate in one direction for ever.

• For \(E = Mg\), one special type of solution is those up-side down stationary ones. Any other solution always has \(|\omega| = \sqrt{E - Mg\cos\theta} > 0\), and so never stops or reverses direction but instead approaches the above stationary point asymptotically both as \(t \to \infty\) and \(t \to \infty\). To visualize concretely, the solution starting at the bottom with \(\theta(0) = 0\), \(\omega(0) = \sqrt{2g/L}\) has \(\theta(t) \to \pm \pi\) and \(\omega(t) \to 0\) as \(t \to \pm\infty\).

Remark 6.3 (Separatrices)

This last kind of special solution is known as a separatrix due to separating the other two qualitatively different sorts of solution. They are also known as heteroclinic orbits, for “asymptotically” starting and ending at different stationary solutions in each time direction — or homoclinic if you consider the angle as a “mod \(2\pi\)” value describing a position, so that \(\theta = \pm\pi\) are the same location and the solutions start and end at the same stationary point.

import numpy as np
# Shortcuts for some favorite mathemtcial functions and numbers:
from numpy import sqrt, sin, cos, pi
from matplotlib.pyplot import figure, plot, grid, title, xlabel, ylabel, legend, show

The Euler’s method code from before does not quite work, but only slight modification is needed; that “scalar” version

def eulerMethod(f, a, b, u_0, n):
    h = (b-a)/n
    t = np.linspace(a, b, n+1)
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        u[i+1] = u[i] + f(t[i], u[i])*h
    return (t, u)

becomes

def eulerMethodSystem(f, a, b, u_0, n):
    """Use Euler's Method to solve du/dt = f_mass_spring(t, u) for t in [a, b], with initial value u(a) = u_0
    Modified from function eulermethod to handle systems.
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)
    
    # Only the following three lines change for the system version
    n_unknowns = len(u_0)
    u = np.zeros([n+1, n_unknowns])
    u[0] = np.array(u_0)  # In case u_0 is a single number (the scalar case)

    for i in range(n):
        u[i+1] = u[i] + f(t[i], u[i])*h
    return (t, u)

6.4.3. Solving the Mass-Spring System

def f_mass_spring(t, u):
    return np.array([ u[1], -(K/M)*u[0] - (D/M)*u[1]])

def E_mass_spring(y, Dy):
    return (K * y**2 + M * Dy**2)/2

def y_mass_spring(t, t_0, u_0, K, M, D):
    (y_0, v_0) = u_0
    discriminant = D**2 - 4*K*M
    if discriminant < 0:  # underdamped
        omega = sqrt(4*K*M - D**2)/(2*M)
        A = y_0
        B = (v_0 + y_0*D/(2*M))/omega
        return exp(-D/(2*M)*(t-t_0)) * ( A*cos(omega*(t-t_0)) + B*sin(omega*(t-t_0)))
    elif discriminant > 0:  # overdamped
        Delta = sqrt(discriminant)
        lambda_plus = (-D + Delta)/(2*M)
        lambda_minus = (-D - Delta)/(2*M)
        A = M*(v_0 - lambda_minus * y_0)/Delta
        B = y_0 - A
        return A*exp(lambda_plus*(t-t_0)) + B*exp(lambda_minus*(t-t_0))
    else:
        q = -D/(2*M)
        A = y_0
        B = v_0 - A * q
        return (A + B*t)*exp(q*(t-t_0))

def damping(K, M, D):
    if D == 0:
        print("Undamped")
    else:
        discriminant = D**2 - 4*K*M
        if discriminant < 0:
            print("Underdamped")
        elif discriminant > 0:
            print("Overdamped")
        else:
            print("Critically damped")

The above functions are available in module numericalMethods; they will be used in later sections.

First solve without damping, so the solutions have sinusoidal solutions

Note: the orbits go clockwise for undamped (and underdamped) systems.

M = 1.0
K = 1.0
D = 0.0
y_0 = 1.0
Dy_0 = 0.0
u_0 = [y_0, Dy_0]
a = 0.0
periods = 4
b = 2 * pi * periods

stepsperperiod = 1000
n = stepsperperiod * periods

(t, U) = eulerMethodSystem(f_mass_spring, a, b, u_0, n)
Y = U[:,0]
DY = U[:,1]

figure(figsize=[14,7])
title(f"y and dy/dt with {K/M=}, {D=} by Euler's method with {stepsperperiod} steps per period")
plot(t, Y, label="y")
plot(t, DY, label="dy/dt")
legend()
xlabel("t")
grid(True)

# Phase plane diagram; for D=0 the exact solutions are ellipses (circles if M = k)

figure(figsize=[8,8])  # Make axes equal length; orbits should be circular or "circular spirals" 
title(f"The orbits of the mass-spring system, {K/M=}, {D=} by Euler's method with {stepsperperiod} steps per period")
plot(Y, DY)
xlabel("y")
ylabel("dy/dt")
plot(Y[0], DY[0], "g*", label="start")
plot(Y[-1], DY[-1], "r*", label="end")
legend()
grid(True)

figure(figsize=[10,4])
E_0 = E_mass_spring(y_0, Dy_0)
E = E_mass_spring(Y, DY)
title("Energy variation")
plot(t, E - E_0)
xlabel("t")
grid(True)

Damped

D = 0.05 # Underdamped: decaying oscillations
#D = 1.1 # Overdamped: exponential decay
(t, U) = eulerMethodSystem(f_mass_spring, a, b, u_0, n)
Y = U[:,0]
DY = U[:,1]

figure(figsize=[14,7])
title(f"y and dy/dt with {K/M=}, {D=} by Euler's method with {stepsperperiod} steps per period")
plot(t, Y, label="y")
plot(t, DY, label="dy/dt")
legend()
xlabel("t")
grid(True)

figure(figsize=[8,8])  # Make axes equal length
title(f"The orbits of the mass-spring system, {K/M=}, {D=} by Euler's method with {stepsperperiod} steps per period")
plot(Y, DY)
xlabel("y")
ylabel("dy/dt")
plot(Y[0], DY[0], "g*", label="start")
plot(Y[-1], DY[-1], "r*", label="end")
legend()
grid(True)

6.4.4. The “Classical” Runge-Kutta Method, Extended to Systems of Equations

As above, the previous “scalar” function for this method needs just three lines of code modified.

Before:

def rungeKutta(f, a, b, u_0, n):
    """Use the (classical) Runge-Kutta Method
    to solve du/dt = f_mass_spring(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.
    u = np.empty_like(t)
    u[0] = u_0
    for i in range(n):
        K_1 = f_mass_spring(t[i], u[i])*h
        K_2 = f_mass_spring(t[i]+h/2, u[i]+K_1/2)*h
        K_3 = f_mass_spring(t[i]+h/2, u[i]+K_2/2)*h
        K_4 = f_mass_spring(t[i]+h, u[i]+K_3)*h
        u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
    return (t, u)

After:

def rungeKuttaSystem(f, a, b, u_0, n):
    """Use the (classical) Runge-Kutta Method
    to solve system du/dt = f_mass_spring(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)  # Note: "n" counts steps, so there are n+1 values for t.

    # Only the following three lines change for the system version.
    n_unknowns = len(u_0)
    u = np.zeros([n+1, n_unknowns])
    u[0] = np.array(u_0)

    for i in range(n):
        K_1 = f_mass_spring(t[i], u[i])*h
        K_2 = f_mass_spring(t[i]+h/2, u[i]+K_1/2)*h
        K_3 = f_mass_spring(t[i]+h/2, u[i]+K_2/2)*h
        K_4 = f_mass_spring(t[i]+h, u[i]+K_3)*h
        u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
    return (t, u)

M = 1.0
k = 1.0
D = 0.0
y_0 = 1.0
Dy_0 = 0.0
u_0 = [y_0, Dy_0]
a = 0.0
periods = 4
b = 2 * pi * periods

stepsperperiod = 25
n = stepsperperiod * periods

(t, U) = rungeKuttaSystem(f_mass_spring, a, b, u_0, n)
y = U[:,0]
Dy = U[:,1]

figure(figsize=[14,7])
title(f"y and dy/dt with {k/M=}, {D=} by Runge-Kutta with {stepsperperiod} steps per period")
plot(t, y, ".:", label="y")
plot(t, Dy, ".:", label="dy/dt")
legend()
xlabel("t")
grid(True)

figure(figsize=[8,8])  # Make axes equal length
title(f"The orbits of the mass-spring system, {k/M=}, {D=} by Runge-Kutta with {stepsperperiod} steps per period")
plot(y, Dy, ".:")
xlabel("y")
ylabel("dy/dt")
plot(y[0], Dy[0], "g*", label="start")
plot(y[-1], Dy[-1], "r*", label="end")
legend()
grid(True)

6.4.5. Appendix: the Explicit Trapezoid and Midpoint Methods for systems

Yet again, the previous functions for these methods need just three lines of code modified.

The demos are just for the non-dissipative case, where the solution is known to be \(y = \cos t\), \(dt/dt = -\sin t\).

For a fairer comparison of “accuracy vs computational effort” to the Runge-Kutta method, twice as many time steps are used so that the same number of function evaluations are used for these three methods.

def explicitTrapezoidSystem(f, a, b, u_0, n):
    """Use the Explict Trapezoid Method (a.k.a Improved Euler)
    to solve system du/dt = f_mass_spring(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)

    # Only the following three lines change for the systems version
    n_unknowns = len(u_0)
    u = np.zeros([n+1, n_unknowns])
    u[0] = np.array(u_0)

    for i in range(n):
        K_1 = f_mass_spring(t[i], u[i])*h
        K_2 = f_mass_spring(t[i]+h, u[i]+K_1)*h
        u[i+1] = u[i] + (K_1 + K_2)/2.
    return (t, u)

M = 1.0
k = 1.0
D = 0.0
y_0 = 1.0
Dy_0 = 0.0
u_0 = [y_0, Dy_0]
a = 0.0
periods = 4
b = 2 * pi * periods

stepsperperiod = 50
n = stepsperperiod * periods

(t, U) = explicitTrapezoidSystem(f_mass_spring, a, b, u_0, n)
y = U[:,0]
Dy = U[:,1]

figure(figsize=[14,7])
title(f"y and dy/dt with {k/M=}, {D=} by explicit trapezoid with {stepsperperiod} steps per period")
plot(t, y, ".:", label="y")
plot(t, Dy, ".:", label="dy/dt")
legend()
xlabel("t")
grid(True)

figure(figsize=[8,8])  # Make axes equal length
title(f"The orbits of the mass-spring system, {k/M=}, {D=} by explicit trapezoid with {stepsperperiod} steps per period")
plot(y, Dy, ":")
xlabel("y")
ylabel("dy/dt")
plot(y[0], Dy[0], "g*", label="start")
plot(y[-1], Dy[-1], "r*", label="end")
legend()
grid(True)

At first glance this is foing well, keeping the orbits circular. However note the discrepancy between the start and end points: these should be the same, as they are (visually) with the Runge-Kutta method.

def explicitMidpointSystem(f, a, b, u_0, n):
    """Use the Explicit Midpoint Method (a.k.a Modified Euler)
    to solve system du/dt = f_mass_spring(t, u) for t in [a, b], with initial value u(a) = u_0
    """
    h = (b-a)/n
    t = np.linspace(a, b, n+1)

    # Only the following three lines change for the systems version.
    n_unknowns = len(u_0)
    u = np.zeros([n+1, n_unknowns])
    u[0] = np.array(u_0)

    for i in range(n):
        K_1 = f_mass_spring(t[i], u[i])*h
        K_2 = f_mass_spring(t[i]+h/2, u[i]+K_1/2)*h
        u[i+1] = u[i] + K_2
    return (t, u)

M = 1.0
k = 1.0
D = 0.0
y_0 = 1.0
Dy_0 = 0.0
u_0 = [y_0, Dy_0]
a = 0.0
periods = 4
b = 2 * pi * periods

stepsperperiod = 50
n = stepsperperiod * periods

(t, U) = explicitMidpointSystem(f_mass_spring, a, b, u_0, n)
y = U[:,0]
Dy = U[:,1]

figure(figsize=[14,7])
title(f"y and dy/dt with {k/M=}, {D=} by explicit midpoint with {stepsperperiod} steps per period")
plot(t, y, ".:", label="y")
plot(t, Dy, ".:", label="dy/dt")
legend()
xlabel("t")
grid(True)

figure(figsize=[8,8])  # Make axes equal length
title(f"The orbits of the mass-spring system, {k/M=}, {D=} by explicit midpoint with {stepsperperiod} steps per period")
plot(y, Dy, ":")
xlabel("y")
ylabel("dy/dt")
plot(y[0], Dy[0], "g*", label="start")
plot(y[-1], Dy[-1], "r*", label="end")
legend()
grid(True)