6.6. Systems of ODEs and Higher Order ODEs#
References:
[Sauer, 2022] Section 6.3, Systems of Ordinary Differential Equations, to Sub-section 6.3.1, Higher order equations.
[Burden et al., 2016] Section 5.9, Higher Order Equations and Systems of Differential Equations.
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.6.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.6.2. Test Cases#
In this and subsequent sections, numerical methods for higher order equations and systems will be compared using several equations seen in the first section of this chapter:
Example 6.5 \(\displaystyle M \frac{d^2 u}{d t^2} = -K u - D \frac{d u}{d t}\) (a Mass-Spring System).
Example 6.6 \(M L \theta'' = -M g\sin\theta - D L \theta'\) (The Freely Rotating Pendulum).
using PyPlot
include("NumericalMethods.jl")
using .NumericalMethods: approx4
The Euler’s method code from before does not quite work, but only slight modification is needed; that “scalar” version
function eulermethod(f, a, b, u_0, n)
h = (b-a)/n
t = range(a, b, n+1)
U = zeros(n+1)
U[1] = u_0
for i in 1:n
U[i+1] = U[i] + f(t[i], U[i])*h
end
return (t, U)
end;
becomes
function eulermethod_system(f, a, b, u_0, n)
# TO DO: one could use multiple dispatch to keep the name "eulermethod".
# The conversion for the system version is mainly "U[i] -> U[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
# The following three lines and the one in the for loop below change for the system version
n_unknowns = length(u_0)
U = zeros(n+1, n_unknowns)
U[1,:] = u_0 # Only for system version
for i in 1:n
U[i+1,:] = U[i,:] + f(t[i], U[i,:])*h # For the system version
end
return (t, U)
end;
Note. Here and below, these notes follow the convention of using lowercase letters for exact solutions; uppercase for numerical approximations.
6.6.3. Solving the Damped Mass-Spring System#
f_mass_spring(t, u) = [ u[2], -(K/M)*u[1] - (D/M)*u[2] ];
E_mass_spring(y, Dy) = (K * y^2 + M * Dy^2)/2;
function y_mass_spring(t; t_0, u_0, K, M, D)
(y_0, v_0) = u_0
discriminant = D^2 - 4K*M
if discriminant < 0 # underdamped
omega = sqrt(4K*M - D^2)/(2M)
A = y_0
B = (v_0 + y_0*D/(2M))/omega
return exp(-D/(2M)*(t-t_0)) * ( A*cos(omega*(t-t_0)) + B*sin(omega*(t-t_0)))
elseif discriminant > 0 # overdamped
Delta = sqrt(discriminant)
lambda_plus = (-D + Delta)/(2M)
lambda_minus = (-D - Delta)/(2M)
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
lambda = -D/(2M)
A = y_0
B = v_0 - A * lambda
return (A + B*t)*exp(lambda*(t-t_0))
end
end;
function damping(K, M, D)
if D == 0
println("Undamped")
else
discriminant = D^2 - 4K*M
if discriminant < 0
println("Underdamped")
elseif discriminant > 0
println("Overdamped")
else
println("Critically damped")
end
end
end;
The above functions are available in module NumericalMethods;
they will be used in later sections.
6.6.3.1. 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 = 2pi * periods
stepsperperiod = 500
n = Int(stepsperperiod * periods)
(t, U) = eulermethod_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
figure(figsize=[10,4])
title("y for K/M=$(K/M), D=$D by Euler's method with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
# Phase plane diagram; for D=0 the exact solutions are ellipses (circles if M = k)
figure(figsize=[6,6]) # Make axes equal length; orbits should be circular or "circular spirals"
title("The orbit")
plot(Y, DY)
xlabel("y")
ylabel("dy/dt")
plot(Y[1], DY[1], "g*", label="start")
plot(Y[end], DY[end], "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)
6.6.3.2. Next solve with damping#
D = 0.5 # Underdamped: decaying oscillations
#D = 2 # Critically damped
#D = 2.1 # Overdamped: exponential decay
periods = 4
b = 2pi * periods
stepsperperiod = 500
n = Int(stepsperperiod * periods)
(t, U) = eulermethod_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
damping(K, M, D)
figure(figsize=[10,4])
title("y for K/M=$(K/M), D=$D by Euler's method with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
# Phase plane diagram; for D=0 the exact solutions are ellipses (circles if M = K)
figure(figsize=[6,6]) # Make axes equal length; orbits should be circular or "circular spirals"
title("The orbit")
plot(Y, DY)
xlabel("y")
ylabel("dy/dt")
plot(Y[1], DY[1], "g*", label="start")
plot(Y[end], DY[end], "r*", label="end")
legend()
grid(true)
Underdamped
6.6.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:
function rungekutta(f, a, b, u_0, n)
# Use the (classical) Runge-Kutta Method to solve
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
h = (b-a)/n
t = range(a, b, n+1)
u = zeros(length(t))
u[1] = u_0
for i in 1:n
K_1 = f(t[i], u[i])*h
K_2 = f(t[i]+h/2, u[i]+K_1/2)*h
K_3 = f(t[i]+h/2, u[i]+K_2/2)*h
K_4 = f(t[i]+h, u[i]+K_3)*h
u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
end
return (t, u)
end;
After:
function rungekutta_system(f, a, b, u_0, n)
# Use the (classical) Runge-Kutta Method to solve
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h
K_3 = f(t[i]+h/2, u[i,:]+K_2/2)*h
K_4 = f(t[i]+h, u[i,:]+K_3)*h
u[i+1,:] = u[i,:] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
end
return (t, u)
end;
M = 1.0
k = 1.0
D = 0.0
u_0 = [1.0, 0.0]
a = 0.0
periods = 4
b = 2pi * periods
stepsperperiod = 25
n = stepsperperiod * periods
(t, U) = rungekutta_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
figure(figsize=[10,4])
title("y for k/M=$(k/M), D=$D by Runge-Kutta with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
# Phase plane diagram; for D=0 the exact solutions are ellipses (circles if M = k)
figure(figsize=[6,6]) # Make axes equal length; orbits should be circular or "circular spirals"
title("The orbit")
plot(Y, DY)
xlabel("y")
ylabel("dy/dt")
plot(Y[1], DY[1], "g*", label="start")
plot(Y[end], DY[end], "r*", label="end")
legend()
grid(true)
D = 0.5 # Underdamped: decaying oscillations
#D = 2 # Critically damped
#D = 2.1 # Overdamped: exponential decay
periods = 4
b = 2pi * periods
stepsperperiod = 25
n = Int(stepsperperiod * periods)
(t, U) = rungekutta_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
damping(k, M, D)
figure(figsize=[10,4])
title("y for k/M=$(k/M), D=$D by Runge-Kutta with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
Underdamped
6.6.5. Solving the Freely Rotating Pendulum Equation#
For now, this is just briefly explored as a cautionary tail of what can happen when slight changes in the system lead to qualitatively very different solution behavior. So we will look at a few examples for the conservative case \(D=0\), close to the separatrix solutions noted above.
Parameters can all be scaled away to \(M = L = g = 1\) so the critical energy is \(Mg = 1\).
f_pendulum(t, u) = [ u[2], -(g/L)*sin(u[1]) ];
M = g = L = 1.0;
E_0 = 1.0 # Separatrix
#E_0 = 0.999
#E_0 = 1.001
theta_0 = 0.0
omega_0 = sqrt(2(E_0 + M*g*cos(theta_0))/(M*L));
u_0 = [theta_0, omega_0]
a = 0.0
#periods = 8 # periods of the linear approximation, "sin(theta) = theta"
#b = 2pi * sqrt(L/g) * periods
b = 80.0
b = 20.;
#stepsperperiod = 1_000
#n = Int(stepsperperiod * periods)
#h = (b-a)/n
stepsperunittime = 10_000
h = 1/stepsperunittime
n = Int(round((b-a)/h))
(t, U) = eulermethod_system(f_pendulum, a, b, u_0, n)
theta = U[:,1]
omega = U[:,2]
figure(figsize=[10,4])
title("By Euler's method with E = $E_0, step size h = $(approx4(h))")
plot(t, theta/pi, label="theta")
xlabel("t")
ylabel(L"\theta/\pi")
grid(true)
# Phase plane diagram
figure(figsize=[10,4])
title("The orbit")
plot(theta/pi, omega)
xlabel(L"\theta/\pi")
ylabel(L"\omega = d\theta/dt")
plot(theta[1]/pi, omega[1], "g*", label="start")
plot(theta[end]/pi, omega[end], "r*", label="end")
legend()
grid(true)
# Error in the (conserved) energy E
figure(figsize=[10,4])
E = (M*L/2) * omega.^2 - M*g*cos.(theta)
E_error = E .- E_0
title("Error in E(t)")
plot(t, E_error, label="theta")
xlabel("t")
#ylabel(L"\theta/\pi")
grid(true)
#stepsperperiod = 10_000
#n = Int(stepsperperiod * periods)
#h = (b-a)/n
stepsperunittime = 25
stepsperunittime = 10_000
h = 1/stepsperunittime
n = Int(round((b-a)/h))
(t, U) = rungekutta_system(f_pendulum, a, b, u_0, n)
theta = U[:,1]
omega = U[:,2]
figure(figsize=[10,4])
title("By the Runge-Kutta method with E = $E_0, step size h = $(approx4(h))")
plot(t, theta/pi, label="theta")
xlabel("t")
ylabel(L"\theta/\pi")
grid(true)
# Phase plane diagram
figure(figsize=[10,4])
title("The orbit")
plot(theta/pi, omega)
xlabel(L"\theta/\pi")
ylabel(L"\omega = d\theta/dt")
plot(theta[1]/pi, omega[1], "g*", label="start")
plot(theta[end]/pi, omega[end], "r*", label="end")
legend()
grid(true)
# Error in the (conserved) energy E
E = (M*L/2) * omega.^2 - M*g*cos.(theta)
E_error = E .- E_0
figure(figsize=[10,4])
title("Error in E(t)")
plot(t, E_error, label="theta")
xlabel("t")
#ylabel(L"\theta/\pi")
grid(true);
6.6.6. 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.
function explicittrapezoid_system(f, a, b, u_0, n)
# Use the Explict Trapezoid Method (a.k.a Improved Euler) to solve the system
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h, u[i,:]+K_1)*h
u[i+1,:] = u[i,:] + (K_1 + K_2)/2.0
end
return (t, u)
end;
D = 0.5 # Underdamped: decaying oscillations
#D = 2 # Critically damped
#D = 2.1 # Overdamped: exponential decay
periods = 4
b = 2pi * periods
stepsperperiod = 50
n = Int(stepsperperiod * periods)
damping(k, M, D)
(t, U) = explicittrapezoid_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
damping(k, M, D)
figure(figsize=[10,4])
title("y for k/M=$(k/M), D=$D by explicit trapezoid with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
Underdamped
Underdamped
At first glance this is doing 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.
function explicitmidpoint_system(f, a, b, u_0, n)
# Use the Explict Midpoint Method (a.k.a Modified Euler) to solve the system
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h
u[i+1,:] = u[i,:] + K_2
end
return (t, u)
end;
D = 0.5 # Underdamped: decaying oscillations
#D = 2 # Critically damped
#D = 2.1 # Overdamped: exponential decay
periods = 4
b = 2pi * periods
stepsperperiod = 50
n = Int(stepsperperiod * periods)
damping(k, M, D)
(t, U) = explicitmidpoint_system(f_mass_spring, a, b, u_0, n)
Y = U[:,1]
DY = U[:,2]
y = y_mass_spring.(t; t_0=a, u_0=u_0, K=K, M=M, D=D) # Exact solution
damping(k, M, D)
figure(figsize=[10,4])
title("y for k/M=$(k/M), D=$D by explicit midpoint with $periods periods, $stepsperperiod steps per period")
plot(t, Y, label="y computed")
plot(t, y, label="exact solution")
xlabel("t")
legend()
grid(true)
figure(figsize=[10,4])
title("Error in Y")
plot(t, y-Y)
xlabel("t")
grid(true)
Underdamped
Underdamped
