# 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

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

with initial conditions

to a first order system, introduction the two functions

Then

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

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

and initial data vector

puts the equation into the form

## 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

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

#### 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

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

and it can be verified that the “energy”

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

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

(Stiffness)

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

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

has first order system form

and the general solution

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

or in system form

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

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\).

(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)
```