34. Systems of ODEs and Higher Order ODEs

Version of 2022-10-11

References:

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

• Section 5.9 Higher Order Equations and Systems of Differential Equations in Burden&Faires

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.

34.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}\]

34.2. Test Cases

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

Motion of a (Damped) Mass-Spring System in One Dimension

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

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

Note on 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.

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.

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

Aside. 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.

using PyPlot
include("NumericalMethods.jl")
using .NumericalMethods: approx4

WARNING: replacing module NumericalMethods.
WARNING: ignoring conflicting import of NumericalMethods.approx4 into Main

34.3. Example A. The Damped Mass-Spring System

f_mass_spring(t, u) = [ u[2], -(k/M)*u[1] - (D/M)*u[2] ];

E_mass_spring(u) = (K * u[1]^2 + M * u[2]^2)/2;

function y_mass_spring(t; t_0, u_0, k, M, D)
    (y_0, v_0) = u_0
        #println("y_0=$y_0, v_0=$v_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
            #println("Underdamped, A=$A, B=$B")
        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
            #println("Overdamped, A=$A, B=$B")
        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
            #println("Critically damped, A=$A, B=$B")
        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.

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.

Solving the (Damped) Mass-Spring System

First, undamped, so with sinusoidal solutions

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

M = 1.0
K = 1.0
D = 0.0
u_0 = [1.0, 0.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(u_0)
E = E_mass_spring.(U)
title("Energy variation")
plot(t, E .- E_0)
xlabel("t")
grid(true)

BoundsError

Stacktrace:
  [1] getindex
    @ ./number.jl:98 [inlined]
  [2] E_mass_spring(u::Float64)
    @ Main ./In[15]:1
  [3] _broadcast_getindex_evalf
    @ ./broadcast.jl:670 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:643 [inlined]
  [5] getindex
    @ ./broadcast.jl:597 [inlined]
  [6] copy
    @ ./broadcast.jl:899 [inlined]
  [7] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(E_mass_spring), Tuple{Matrix{Float64}}})
    @ Base.Broadcast ./broadcast.jl:860
  [8] top-level scope
    @ In[20]:44
  [9] eval
    @ ./boot.jl:368 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1428

Damped

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

34.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

Example B. The Freely Rotating Pendulum (Under Construction)

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.

20.0

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

34.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.

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 cicular. However, note the discrepancy between the start and end points: these should be the same, as they are (visually) the 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

---

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International