Reduction of Order - Executable Engineering: Interactive Numerical Methods in Python

Up to this point, we have focused on first-order differential equations. Any higher-order ordinary differential equation can be converted into a coupled system of first-order equations via reduction of order by introducing auxiliary state variables for the lower-order derivatives.

For example, consider a general second-order ODE:

\frac{\mathrm{d}^2 y}{\mathrm{d}t^2} = f\!\left(t, y, \frac{\mathrm{d}y}{\mathrm{d}t}\right)

(1)

Defining the auxiliary velocity variable $z = \frac{\mathrm{d}y}{\mathrm{d}t}$ transforms the single second-order equation into a system of two first-order equations:

\begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t} &= z \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= f(t, y, z) \end{aligned}

(2)

which can be integrated directly using any vector Runge-Kutta or multistep solver!

import numpy as np
from scipy.integrate import solve_ivp

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Example: Swinging Pendulum¶

The nonlinear equation of motion for a rigid pendulum of mass $m$ and length $l$ under gravity $g$ is:

m l \frac{\mathrm{d}^2\Theta(t)}{\mathrm{d}t^2} = -m g \sin\!\big(\Theta(t)\big)

(3)

Applying reduction of order, we define the two-component state vector $\mathbf{y}(t)$ :

\mathbf{y}(t) = \begin{bmatrix} y_1(t) \\ y_2(t) \end{bmatrix} = \begin{bmatrix} \Theta(t) \\ \dot{\Theta}(t) \end{bmatrix}

(4)

Differentiating $\mathbf{y}(t)$ with respect to time gives the first-order vector ODE system:

\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \begin{bmatrix} \dot{\Theta}(t) \\ \ddot{\Theta}(t) \end{bmatrix} = \begin{bmatrix} y_2 \\ -\frac{g}{l} \sin(y_1) \end{bmatrix} = \mathbf{f}(t, \mathbf{y})

(5)

def pendulum_rhs(t, y, g=9.81, l=1.0):
    theta, theta_dot = y
    return np.array([theta_dot, -(g / l) * np.sin(theta)])

y0 = [np.pi / 4.0, 0.0]
sol = solve_ivp(pendulum_rhs, (0.0, 10.0), y0, method="RK45", max_step=0.05)

exe.plot_pendulum_system(sol.t, sol.y[0], sol.y[1])