Ordinary Differential Equations (17): Physics and Engineering Applications

Differential equations are not a pure mathematical game — they are the language for understanding the physical world. From celestial motion to circuit response, from spring vibration to chemical reactions, almost all dynamical system behaviors can be described by differential equations. In this chapter, we explore the core applications of ordinary differential equations in physics and engineering, building a bridge from mathematics to practice.

Differential Equations in Classical Mechanics

Mathematical Expression of Newton's Laws of Motion

Newton's second lawis the cornerstone of mechanics, but it is essentially a differential equation:Hereis the position vector, and the forcemay depend on position, velocity, and time.

Important observation: Given the expression for force, an object's motion is completely determined by its initial position and initial velocity. This embodies the deterministic nature of Newtonian mechanics.

Free Fall with Air Resistance

Consider an object falling from height, subject to both gravity and air resistance:whereis the drag coefficient andis the velocity (downward positive).

Terminal velocity: When, acceleration is zero, reaching terminal velocity:

Solution process: This is a separable first-order ODE:Using partial fraction decomposition:After integration:

(assuming zero initial velocity)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

def falling_with_drag():
    """Free fall with air resistance"""
    m = 70  # Mass kg (human body)
    g = 9.8  # Gravitational acceleration
    k = 0.25  # Drag coefficient (depends on body size and posture)
    
    v_terminal = np.sqrt(m * g / k)
    print(f"Terminal velocity: {v_terminal:.1f} m/s ({v_terminal * 3.6:.1f} km/h)")
    
    def drag_ode(v, t):
        return g - (k/m) * v**2
    
    t = np.linspace(0, 30, 500)
    v0 = 0
    
    # Numerical solution
    v_numerical = odeint(drag_ode, v0, t).flatten()
    
    # Analytical solution
    v_analytical = v_terminal * np.tanh(g * t / v_terminal)
    
    # Position (integrate velocity)
    x_numerical = np.cumsum(v_numerical) * (t[1] - t[0])
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Velocity-time
    ax1 = axes[0]
    ax1.plot(t, v_numerical, 'b-', linewidth=2, label='Numerical')
    ax1.plot(t, v_analytical, 'r--', linewidth=2, label='Analytical')
    ax1.axhline(y=v_terminal, color='k', linestyle=':', label=f'Terminal: {v_terminal:.1f} m/s')
    ax1.set_xlabel('Time (s)', fontsize=12)
    ax1.set_ylabel('Velocity (m/s)', fontsize=12)
    ax1.set_title('Velocity vs Time', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Position-time
    ax2 = axes[1]
    ax2.plot(t, x_numerical, 'g-', linewidth=2)
    ax2.set_xlabel('Time (s)', fontsize=12)
    ax2.set_ylabel('Distance (m)', fontsize=12)
    ax2.set_title('Distance vs Time', fontsize=14, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # Comparison of different drag coefficients
    ax3 = axes[2]
    k_values = [0.1, 0.25, 0.5, 1.0]
    for k_val in k_values:
        v_term = np.sqrt(m * g / k_val)
        v = v_term * np.tanh(g * t / v_term)
        ax3.plot(t, v, linewidth=2, label=f'k = {k_val}')
    
    ax3.set_xlabel('Time (s)', fontsize=12)
    ax3.set_ylabel('Velocity (m/s)', fontsize=12)
    ax3.set_title('Effect of Drag Coefficient', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('falling_with_drag.png', dpi=150, bbox_inches='tight')
    plt.show()

falling_with_drag()

Planetary Motion and the Kepler Problem

The two-body problem is one of the most important problems in classical mechanics. Consider a planet of massorbiting a star of mass():Introducing polar coordinatesand using conservation of angular momentum, we obtain the Binet equation:where.

This is a simple harmonic oscillator equation! The general solution is:This is the polar equation of a conic section, whereis the eccentricity: -: Ellipse (bound orbit) -: Parabola (critical escape) -: Hyperbola (escape orbit)

def kepler_orbits():
    """Visualization of Kepler orbits"""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # Orbits with different eccentricities
    ax1 = axes[0]
    theta = np.linspace(0, 2*np.pi, 1000)
    
    eccentricities = [0, 0.3, 0.6, 0.9]
    p = 1  # Semi-latus rectum
    
    for e in eccentricities:
        r = p / (1 + e * np.cos(theta))
        # Only plot bound portion
        valid = r > 0
        x = r[valid] * np.cos(theta[valid])
        y = r[valid] * np.sin(theta[valid])
        ax1.plot(x, y, linewidth=2, label=f'e = {e}')
    
    # Parabola and hyperbola
    for e, style in [(1.0, '--'), (1.5, ':')]:
        theta_range = np.linspace(-np.pi*0.8, np.pi*0.8, 500)
        r = p / (1 + e * np.cos(theta_range))
        valid = r > 0
        x = r[valid] * np.cos(theta_range[valid])
        y = r[valid] * np.sin(theta_range[valid])
        ax1.plot(x, y, linestyle=style, linewidth=2, label=f'e = {e}')
    
    ax1.plot(0, 0, 'yo', markersize=15, label='Sun')
    ax1.set_xlim(-4, 2)
    ax1.set_ylim(-3, 3)
    ax1.set_aspect('equal')
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title('Kepler Orbits: Different Eccentricities', fontsize=14, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    
    # Numerical orbit simulation
    ax2 = axes[1]
    
    def kepler_ode(t, state, GM):
        x, y, vx, vy = state
        r = np.sqrt(x**2 + y**2)
        ax = -GM * x / r**3
        ay = -GM * y / r**3
        return [vx, vy, ax, ay]
    
    GM = 1.0
    
    # Different initial conditions
    initial_conditions = [
        ([1, 0, 0, 1.0], 'Circle'),
        ([1, 0, 0, 0.8], 'Ellipse'),
        ([1, 0, 0, 1.414], 'Parabola-like'),
    ]
    
    for (state0, name) in initial_conditions:
        t_span = (0, 20)
        t_eval = np.linspace(0, 20, 2000)
        
        sol = solve_ivp(kepler_ode, t_span, state0, args=(GM,), 
                       t_eval=t_eval, method='DOP853', rtol=1e-10)
        
        ax2.plot(sol.y[0], sol.y[1], linewidth=2, label=name)
    
    ax2.plot(0, 0, 'yo', markersize=15)
    ax2.set_xlim(-3, 2)
    ax2.set_ylim(-2, 2)
    ax2.set_aspect('equal')
    ax2.set_xlabel('x', fontsize=12)
    ax2.set_ylabel('y', fontsize=12)
    ax2.set_title('Numerical Orbit Simulation', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('kepler_orbits.png', dpi=150, bbox_inches='tight')
    plt.show()

kepler_orbits()

In-Depth Analysis of Vibration Systems

Multi-Degree-of-Freedom Vibration Systems

Real vibration systems often have multiple degrees of freedom. Considermasses connected by springs:where,, andare mass, damping, and stiffness matrices respectively.

Modal analysis: For undamped free vibrationAssuming solution form, we obtain the eigenvalue problem:Eigenvaluescorrespond to natural frequencies, and eigenvectorscorrespond to mode shapes.

Three-Degree-of-Freedom Spring System

def multi_dof_vibration():
    """Modal analysis of three-degree-of-freedom vibration system"""
    from scipy.linalg import eigh
    
    # Mass matrix (assuming unit mass)
    M = np.array([[1, 0, 0],
                  [0, 1, 0],
                  [0, 0, 1]])
    
    # Stiffness matrix (three masses connected by springs, fixed at both ends)
    k = 1.0
    K = np.array([[2*k, -k, 0],
                  [-k, 2*k, -k],
                  [0, -k, 2*k]])
    
    # Solve generalized eigenvalue problem
    eigenvalues, eigenvectors = eigh(K, M)
    frequencies = np.sqrt(eigenvalues)
    
    print("Natural frequencies:", frequencies)
    print("\nMode shapes (column vectors):")
    print(eigenvectors)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Plot mode shapes
    ax1 = axes[0, 0]
    x_pos = [0, 1, 2, 3]  # Equilibrium positions
    
    for i, mode in enumerate(eigenvectors.T):
        # Normalize mode shape
        mode_normalized = mode / np.max(np.abs(mode))
        displacement = np.concatenate([[0], mode_normalized, [0]])  # Fixed boundaries
        ax1.plot(x_pos + [4], displacement, 'o-', linewidth=2, markersize=10,
                label=f'Mode {i+1}: ω = {frequencies[i]:.3f}')
    
    ax1.axhline(y=0, color='k', linewidth=0.5)
    ax1.set_xlabel('Mass position', fontsize=12)
    ax1.set_ylabel('Displacement (normalized)', fontsize=12)
    ax1.set_title('Mode Shapes', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Time domain response (excite first mass)
    ax2 = axes[0, 1]
    
    def three_dof_system(t, state):
        x = state[:3]
        v = state[3:]
        a = -np.linalg.solve(M, K @ x)
        return np.concatenate([v, a])
    
    # Initial condition: displace only first mass
    state0 = [1, 0, 0, 0, 0, 0]
    t_span = (0, 50)
    t_eval = np.linspace(0, 50, 2000)
    
    sol = solve_ivp(three_dof_system, t_span, state0, t_eval=t_eval, method='RK45')
    
    for i in range(3):
        ax2.plot(sol.t, sol.y[i], linewidth=1.5, label=f'Mass {i+1}')
    
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Displacement', fontsize=12)
    ax2.set_title('Free Vibration Response', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Energy transfer
    ax3 = axes[1, 0]
    
    # Calculate kinetic energy of each mass
    KE = np.zeros((3, len(sol.t)))
    for i in range(3):
        KE[i] = 0.5 * M[i, i] * sol.y[i+3]**2
    
    for i in range(3):
        ax3.plot(sol.t, KE[i], linewidth=1.5, label=f'Mass {i+1}')
    
    ax3.plot(sol.t, np.sum(KE, axis=0), 'k--', linewidth=2, label='Total KE')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Kinetic Energy', fontsize=12)
    ax3.set_title('Energy Distribution', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Phase space trajectory (first mass)
    ax4 = axes[1, 1]
    ax4.plot(sol.y[0], sol.y[3], 'b-', linewidth=1, alpha=0.7)
    ax4.plot(sol.y[0][0], sol.y[3][0], 'go', markersize=10, label='Start')
    ax4.set_xlabel('$x_1$', fontsize=12)
    ax4.set_ylabel('$v_1$', fontsize=12)
    ax4.set_title('Phase Portrait (Mass 1)', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('multi_dof_vibration.png', dpi=150, bbox_inches='tight')
    plt.show()

multi_dof_vibration()

Nonlinear Vibration: The Duffing Oscillator

When displacement is large, springs may exhibit nonlinearity:This is the famous Duffing equation.

-: Hardening spring (stiffens with extension) -: Softening spring (softens with extension)

The Duffing system exhibits rich dynamical behavior: - Multiple steady-state solutions (jump phenomenon) - Subharmonic and superharmonic resonance - Chaotic motion

def duffing_oscillator():
    """Dynamics of the Duffing oscillator"""
    
    def duffing(t, state, delta, alpha, beta, gamma, omega):
        x, v = state
        return [v, -delta*v - alpha*x - beta*x**3 + gamma*np.cos(omega*t)]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Parameters
    delta = 0.3  # Damping
    alpha = -1.0  # Linear term (double-well potential)
    beta = 1.0   # Nonlinear term
    omega = 1.2  # Driving frequency
    
    # Different driving amplitudes
    gamma_values = [0.2, 0.3, 0.37, 0.5]
    
    for ax, gamma in zip(axes.flatten(), gamma_values):
        t_span = (0, 500)
        t_eval = np.linspace(400, 500, 2000)  # Only steady state
        
        sol = solve_ivp(duffing, t_span, [0.1, 0], args=(delta, alpha, beta, gamma, omega),
                       t_eval=t_eval, method='RK45')
        
        ax.plot(sol.y[0], sol.y[1], 'b-', linewidth=0.5, alpha=0.7)
        ax.set_xlabel('x', fontsize=12)
        ax.set_ylabel('v', fontsize=12)
        ax.set_title(f'γ = {gamma}', fontsize=14, fontweight='bold')
        ax.grid(True, alpha=0.3)
    
    plt.suptitle('Duffing Oscillator: Transition to Chaos', fontsize=16, fontweight='bold', y=1.02)
    plt.tight_layout()
    plt.savefig('duffing_oscillator.png', dpi=150, bbox_inches='tight')
    plt.show()

duffing_oscillator()

Circuits and Electromagnetism

General RLC Networks

Complex circuit networks can be described by differential equation systems. Using the state-space approach:

Choose capacitor voltagesand inductor currentsas state variables:$

$ Using Kirchhoff's laws to establish relationships between,and state variables:whereandis the input (voltage or current source).

Coupled Circuits: Transformers

Transformers connect two circuits through magnetic coupling:$

$ whereis the mutual inductance, satisfying.

def coupled_circuits():
    """Analysis of coupled circuits (transformer)"""
    
    # Parameters
    L1 = 1e-3   # Primary inductance 1mH
    L2 = 4e-3   # Secondary inductance 4mH
    M = 1.5e-3  # Mutual inductance
    R1 = 10     # Primary resistance
    R2 = 100    # Secondary resistance (load)
    
    # Coupling coefficient
    k = M / np.sqrt(L1 * L2)
    print(f"Coupling coefficient k = {k:.2f}")
    
    # Turns ratio
    n = np.sqrt(L2 / L1)
    print(f"Turns ratio n = {n:.1f}")
    
    def transformer(t, state, V0, omega):
        i1, i2 = state
        # v1 = V0*sin(ω t) = L1*di1/dt + M*di2/dt + R1*i1
        # 0 = M*di1/dt + L2*di2/dt + R2*i2
        
        det = L1*L2 - M**2
        v1 = V0 * np.sin(omega * t)
        
        di1_dt = (L2*(v1 - R1*i1) + M*R2*i2) / det
        di2_dt = (-M*(v1 - R1*i1) - L1*R2*i2) / det
        
        return [di1_dt, di2_dt]
    
    V0 = 10  # Input voltage amplitude
    omega = 2*np.pi*1000  # 1kHz
    
    t_span = (0, 0.02)  # 20ms
    t_eval = np.linspace(0, 0.02, 2000)
    
    sol = solve_ivp(transformer, t_span, [0, 0], args=(V0, omega),
                   t_eval=t_eval, method='RK45')
    
    # Calculate voltages
    v1 = V0 * np.sin(omega * sol.t)
    v2 = R2 * sol.y[1]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Currents
    ax1 = axes[0, 0]
    ax1.plot(sol.t*1000, sol.y[0]*1000, 'b-', linewidth=2, label='$i_1$(primary)')
    ax1.plot(sol.t*1000, sol.y[1]*1000, 'r-', linewidth=2, label='$i_2$(secondary)')
    ax1.set_xlabel('Time (ms)', fontsize=12)
    ax1.set_ylabel('Current (mA)', fontsize=12)
    ax1.set_title('Transformer Currents', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Voltages
    ax2 = axes[0, 1]
    ax2.plot(sol.t*1000, v1, 'b-', linewidth=2, label='$v_1$(input)')
    ax2.plot(sol.t*1000, v2, 'r-', linewidth=2, label='$v_2$(output)')
    ax2.set_xlabel('Time (ms)', fontsize=12)
    ax2.set_ylabel('Voltage (V)', fontsize=12)
    ax2.set_title('Transformer Voltages', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Power
    ax3 = axes[1, 0]
    p1 = v1 * sol.y[0]  # Input power
    p2 = v2 * sol.y[1]  # Output power
    ax3.plot(sol.t*1000, p1*1000, 'b-', linewidth=2, label='Input power')
    ax3.plot(sol.t*1000, p2*1000, 'r-', linewidth=2, label='Output power')
    ax3.set_xlabel('Time (ms)', fontsize=12)
    ax3.set_ylabel('Power (mW)', fontsize=12)
    ax3.set_title('Power Transfer', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Frequency response
    ax4 = axes[1, 1]
    freqs = np.logspace(2, 5, 200)
    gain = []
    
    for f in freqs:
        omega_i = 2*np.pi*f
        # Steady-state analysis (using phasor method)
        # Simplified: assume steady state, calculate gain
        s = 1j*omega_i
        Z11 = R1 + s*L1
        Z12 = s*M
        Z21 = s*M
        Z22 = R2 + s*L2
        
        # Output voltage / input voltage
        H = -Z21*R2 / (Z11*Z22 - Z12*Z21)
        gain.append(np.abs(H))
    
    ax4.semilogx(freqs, 20*np.log10(gain), 'g-', linewidth=2)
    ax4.set_xlabel('Frequency (Hz)', fontsize=12)
    ax4.set_ylabel('Gain (dB)', fontsize=12)
    ax4.set_title('Frequency Response', fontsize=14, fontweight='bold')
    ax4.grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('coupled_circuits.png', dpi=150, bbox_inches='tight')
    plt.show()

coupled_circuits()

Nonlinear Circuits: Diodes and Transistors

Real circuit elements are often nonlinear. For example, a diode's I-V characteristic:whereis the reverse saturation current andis the thermal voltage.

Circuit equation with a diode:Such nonlinear ODEs usually can only be solved numerically.

Heat Conduction and Diffusion

Newton's Law of Cooling

Heat exchange between an object and its environment follows:The solution is:

Application: Estimating time of death in forensic science.

Multi-Region Heat Conduction

Consider two objects in contact through an interface:$

$ This is a coupled linear system.

def thermal_conduction():
    """Heat conduction problems"""
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Simple Newton cooling
    ax1 = axes[0]
    T_env = 20  # Environmental temperature
    T0 = 90     # Initial temperature
    k_values = [0.05, 0.1, 0.2, 0.5]
    t = np.linspace(0, 60, 500)
    
    for k in k_values:
        T = T_env + (T0 - T_env) * np.exp(-k * t)
        ax1.plot(t, T, linewidth=2, label=f'k = {k}')
    
    ax1.axhline(y=T_env, color='k', linestyle='--', label='Environment')
    ax1.set_xlabel('Time (min)', fontsize=12)
    ax1.set_ylabel('Temperature (° C)', fontsize=12)
    ax1.set_title('Newton\'s Law of Cooling', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Two-region heat conduction
    ax2 = axes[1]
    
    def two_region(t, state):
        T1, T2 = state
        C1, C2 = 100, 50  # Heat capacities
        h = 5    # Interface heat transfer coefficient
        k = 1    # Environment heat transfer coefficient
        T_env = 20
        
        dT1 = -h/C1 * (T1 - T2)
        dT2 = h/C2 * (T1 - T2) - k/C2 * (T2 - T_env)
        return [dT1, dT2]
    
    sol = solve_ivp(two_region, [0, 120], [100, 50], t_eval=np.linspace(0, 120, 500))
    
    ax2.plot(sol.t, sol.y[0], 'r-', linewidth=2, label='Region 1')
    ax2.plot(sol.t, sol.y[1], 'b-', linewidth=2, label='Region 2')
    ax2.axhline(y=20, color='k', linestyle='--', label='Environment')
    ax2.set_xlabel('Time (min)', fontsize=12)
    ax2.set_ylabel('Temperature (° C)', fontsize=12)
    ax2.set_title('Two-Region Heat Transfer', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Forensic application: time of death estimation
    ax3 = axes[2]
    
    # Inverse calculation
    T_body = 37  # Normal body temperature
    T_env = 20
    T_measured = 30  # Measured temperature
    k = 0.1  # Assumed cooling constant
    
    # t = -ln((T_measured - T_env)/(T_body - T_env))/k
    t_death = -np.log((T_measured - T_env)/(T_body - T_env))/k
    print(f"Estimated time of death: {t_death:.1f} hours ago")
    
    t = np.linspace(-t_death-2, 10, 500)
    T = T_env + (T_body - T_env) * np.exp(-k * (t + t_death))
    
    ax3.plot(t, T, 'g-', linewidth=2)
    ax3.axvline(x=0, color='r', linestyle='--', label='Measurement time')
    ax3.axvline(x=-t_death, color='b', linestyle='--', label=f'Death time (~{t_death:.1f}h ago)')
    ax3.axhline(y=T_measured, color='gray', linestyle=':', alpha=0.5)
    ax3.plot(0, T_measured, 'ro', markersize=10)
    ax3.set_xlabel('Time (hours)', fontsize=12)
    ax3.set_ylabel('Body Temperature (° C)', fontsize=12)
    ax3.set_title('Forensic Application', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('thermal_conduction.png', dpi=150, bbox_inches='tight')
    plt.show()

thermal_conduction()

Chemical Reaction Kinetics

First-Order Reactions

Radioactive decay and many chemical reactions follow first-order kinetics:Solution: Half-life:

Consecutive Reactions

Consider:$

$$

Solutions: - - -

Reversible Reactions and Chemical Equilibrium:Equilibrium constant.

def chemical_kinetics():
    """Chemical reaction kinetics"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # First-order reaction (half-life)
    ax1 = axes[0, 0]
    t = np.linspace(0, 10, 500)
    
    half_lives = [1, 2, 4]
    for t_half in half_lives:
        k = np.log(2) / t_half
        N = np.exp(-k * t)
        ax1.plot(t, N, linewidth=2, label=f'$t_{{ 1/2 }}$= {t_half}')
    
    ax1.axhline(y=0.5, color='k', linestyle='--', alpha=0.5)
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('$N/N_0$', fontsize=12)
    ax1.set_title('First-Order Kinetics', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Consecutive reactions
    ax2 = axes[0, 1]
    
    k1, k2 = 1.0, 0.5
    A0 = 1.0
    
    A = A0 * np.exp(-k1 * t)
    B = k1 * A0 / (k2 - k1) * (np.exp(-k1 * t) - np.exp(-k2 * t))
    C = A0 - A - B
    
    ax2.plot(t, A, 'b-', linewidth=2, label='[A]')
    ax2.plot(t, B, 'g-', linewidth=2, label='[B]')
    ax2.plot(t, C, 'r-', linewidth=2, label='[C]')
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Concentration', fontsize=12)
    ax2.set_title('Consecutive Reactions: A → B → C', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Reversible reaction
    ax3 = axes[1, 0]
    
    def reversible(t, state, kf, kr):
        A, B = state
        return [-kf*A + kr*B, kf*A - kr*B]
    
    kf, kr = 1.0, 0.3
    K_eq = kf / kr
    A_eq = 1 / (1 + K_eq)
    B_eq = K_eq / (1 + K_eq)
    
    sol = solve_ivp(reversible, [0, 10], [1, 0], args=(kf, kr),
                   t_eval=np.linspace(0, 10, 500))
    
    ax3.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='[A]')
    ax3.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='[B]')
    ax3.axhline(y=A_eq, color='b', linestyle='--', alpha=0.5)
    ax3.axhline(y=B_eq, color='r', linestyle='--', alpha=0.5)
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Concentration', fontsize=12)
    ax3.set_title(f'Reversible Reaction:$K_{{ eq }}$= {K_eq:.2f}', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Michaelis-Menten enzyme kinetics
    ax4 = axes[1, 1]
    
    def michaelis_menten(t, state, Vmax, Km):
        S = state[0]
        return [-Vmax * S / (Km + S)]
    
    Vmax = 1.0
    Km_values = [0.1, 0.5, 2.0]
    
    for Km in Km_values:
        sol = solve_ivp(michaelis_menten, [0, 10], [1.0], args=(Vmax, Km),
                       t_eval=np.linspace(0, 10, 500))
        ax4.plot(sol.t, sol.y[0], linewidth=2, label=f'$K_m$= {Km}')
    
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('[S]', fontsize=12)
    ax4.set_title('Michaelis-Menten Kinetics', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('chemical_kinetics.png', dpi=150, bbox_inches='tight')
    plt.show()

chemical_kinetics()

Biological Systems and Population Dynamics

Lotka-Volterra Predator-Prey Model$

-: Prey population -: Predator population -: Prey growth rate -: Predation rate -: Predator death rate -: Predator conversion efficiency

Equilibrium analysis:

1.: Extinction (unstable saddle) 2.: Coexistence (center, periodic solutions)

Competitive Exclusion Principle

Two species competing for the same resource:$

Competition coefficientsdetermine coexistence or exclusion.

def population_dynamics():
    """Population dynamics models"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Lotka-Volterra
    ax1 = axes[0, 0]
    
    def lotka_volterra(t, state, alpha, beta, delta, gamma):
        x, y = state
        return [alpha*x - beta*x*y, delta*x*y - gamma*y]
    
    alpha, beta, delta, gamma = 1.0, 0.1, 0.05, 0.5
    
    t_span = (0, 100)
    t_eval = np.linspace(0, 100, 2000)
    
    for x0, y0 in [(10, 2), (20, 5), (5, 10)]:
        sol = solve_ivp(lotka_volterra, t_span, [x0, y0], 
                       args=(alpha, beta, delta, gamma), t_eval=t_eval)
        ax1.plot(sol.y[0], sol.y[1], linewidth=1.5)
    
    # Equilibrium point
    x_eq, y_eq = gamma/delta, alpha/beta
    ax1.plot(x_eq, y_eq, 'r*', markersize=15, label='Equilibrium')
    
    ax1.set_xlabel('Prey (x)', fontsize=12)
    ax1.set_ylabel('Predator (y)', fontsize=12)
    ax1.set_title('Lotka-Volterra Phase Portrait', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Time evolution
    ax2 = axes[0, 1]
    sol = solve_ivp(lotka_volterra, t_span, [10, 2], 
                   args=(alpha, beta, delta, gamma), t_eval=t_eval)
    
    ax2.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Prey')
    ax2.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Predator')
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Population', fontsize=12)
    ax2.set_title('Population Dynamics', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Competition model
    ax3 = axes[1, 0]
    
    def competition(t, state, r1, r2, K1, K2, a12, a21):
        x, y = state
        return [r1*x*(1 - (x + a12*y)/K1),
                r2*y*(1 - (y + a21*x)/K2)]
    
    # Coexistence case
    r1, r2 = 1.0, 1.0
    K1, K2 = 100, 100
    a12, a21 = 0.5, 0.5
    
    sol = solve_ivp(competition, [0, 50], [10, 10], 
                   args=(r1, r2, K1, K2, a12, a21),
                   t_eval=np.linspace(0, 50, 500))
    
    ax3.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Species 1')
    ax3.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Species 2')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Population', fontsize=12)
    ax3.set_title('Competition: Coexistence', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Exclusion case
    ax4 = axes[1, 1]
    a12, a21 = 1.5, 1.5
    
    sol = solve_ivp(competition, [0, 50], [10, 9], 
                   args=(r1, r2, K1, K2, a12, a21),
                   t_eval=np.linspace(0, 50, 500))
    
    ax4.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Species 1')
    ax4.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Species 2')
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('Population', fontsize=12)
    ax4.set_title('Competition: Exclusion', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('population_dynamics.png', dpi=150, bbox_inches='tight')
    plt.show()

population_dynamics()

Epidemic Model: SIR$

$$

The basic reproduction numberdetermines whether an epidemic occurs.

def sir_model():
    """SIR epidemic model"""
    
    def sir(t, state, beta, gamma):
        S, I, R = state
        N = S + I + R
        return [-beta*S*I/N, beta*S*I/N - gamma*I, gamma*I]
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Basic SIR dynamics
    ax1 = axes[0]
    
    beta = 0.3   # Transmission rate
    gamma = 0.1  # Recovery rate
    N = 1000
    I0 = 1
    R0 = 0
    S0 = N - I0 - R0
    
    R_0 = beta * S0 / (gamma * N)
    print(f"Basic reproduction number R ₀ = {R_0:.2f}")
    
    sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                   t_eval=np.linspace(0, 200, 1000))
    
    ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Susceptible')
    ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Infected')
    ax1.plot(sol.t, sol.y[2], 'g-', linewidth=2, label='Recovered')
    ax1.set_xlabel('Time (days)', fontsize=12)
    ax1.set_ylabel('Population', fontsize=12)
    ax1.set_title(f'SIR Model ($R_0$= {R_0:.1f})', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Effect of different R0
    ax2 = axes[1]
    
    R0_values = [1.5, 2.5, 4.0]
    
    for R0_val in R0_values:
        beta = R0_val * gamma * N / S0
        sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                       t_eval=np.linspace(0, 200, 1000))
        ax2.plot(sol.t, sol.y[1], linewidth=2, label=f'$R_0$= {R0_val}')
    
    ax2.set_xlabel('Time (days)', fontsize=12)
    ax2.set_ylabel('Infected', fontsize=12)
    ax2.set_title('Effect of$R_0$on Peak', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Phase plane
    ax3 = axes[2]
    
    beta = 0.3
    sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                   t_eval=np.linspace(0, 200, 1000))
    
    ax3.plot(sol.y[0], sol.y[1], 'b-', linewidth=2)
    ax3.plot(S0, I0, 'go', markersize=10, label='Start')
    ax3.plot(sol.y[0][-1], sol.y[1][-1], 'ro', markersize=10, label='End')
    
    # Threshold line
    S_threshold = gamma * N / beta
    ax3.axvline(x=S_threshold, color='k', linestyle='--', label=f'Threshold S = {S_threshold:.0f}')
    
    ax3.set_xlabel('Susceptible', fontsize=12)
    ax3.set_ylabel('Infected', fontsize=12)
    ax3.set_title('SIR Phase Plane', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('sir_model.png', dpi=150, bbox_inches='tight')
    plt.show()

sir_model()

Fluid Mechanics ODEs

Simplification of Navier-Stokes Equations

The full Navier-Stokes equations are partial differential equations, but in certain symmetric cases they simplify to ODEs.

Poiseuille flow (pipe flow): Steady laminar flow satisfiesThe solution is a parabolic velocity profile:

Boundary Layer Equation

The Blasius equation describes the flat plate boundary layer:Boundary conditions:,This is a third-order nonlinear BVP (boundary value problem).

def fluid_mechanics():
    """ODEs in fluid mechanics"""
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Poiseuille flow
    ax1 = axes[0]
    
    R = 1  # Pipe radius
    r = np.linspace(0, R, 100)
    
    # Normalized velocity profile
    u = 1 - (r/R)**2
    
    ax1.plot(u, r, 'b-', linewidth=2)
    ax1.plot(u, -r, 'b-', linewidth=2)
    ax1.fill_betweenx(r, 0, u, alpha=0.3)
    ax1.fill_betweenx(-r, 0, u, alpha=0.3)
    
    ax1.axvline(x=0, color='k', linewidth=2)
    ax1.axhline(y=R, color='k', linewidth=2)
    ax1.axhline(y=-R, color='k', linewidth=2)
    
    ax1.set_xlabel('Velocity u/u_max', fontsize=12)
    ax1.set_ylabel('Radial position r', fontsize=12)
    ax1.set_title('Poiseuille Flow Profile', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # Blasius boundary layer
    ax2 = axes[1]
    
    from scipy.integrate import solve_bvp
    
    def blasius_ode(eta, y):
        f, f_prime, f_double_prime = y
        return [f_prime, f_double_prime, -0.5 * f * f_double_prime]
    
    def blasius_bc(ya, yb):
        return [ya[0], ya[1], yb[1] - 1]
    
    eta = np.linspace(0, 10, 100)
    y_init = np.zeros((3, eta.size))
    y_init[1] = eta / 10  # Initial guess
    
    sol = solve_bvp(blasius_ode, blasius_bc, eta, y_init)
    
    ax2.plot(sol.y[1], sol.x, 'b-', linewidth=2, label="$f'(\\eta)$= u/U")
    ax2.axhline(y=5, color='k', linestyle='--', alpha=0.5, label='δ₉₉')
    ax2.set_xlabel("$f'$(normalized velocity)", fontsize=12)
    ax2.set_ylabel("$\\eta$(similarity variable)", fontsize=12)
    ax2.set_title('Blasius Boundary Layer', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 1.2)
    
    # Boundary layer thickness vs x
    ax3 = axes[2]
    
    x = np.linspace(0.1, 10, 100)
    Re_x = 1000 * x  # Assume some Re number
    delta = 5 * x / np.sqrt(Re_x)  # Blasius result
    
    ax3.plot(x, delta, 'g-', linewidth=2)
    ax3.set_xlabel('x (distance from leading edge)', fontsize=12)
    ax3.set_ylabel('δ (boundary layer thickness)', fontsize=12)
    ax3.set_title('Boundary Layer Growth', fontsize=14, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fluid_mechanics.png', dpi=150, bbox_inches='tight')
    plt.show()

fluid_mechanics()

Gyroscope and Rigid Body Dynamics

Euler Equations

Euler's equations for rigid body rotation about a fixed point:$

$ This is a set of nonlinear coupled ODEs.

Torque-Free Symmetric Top

Whenand no external torque: where.

The solution is circular motion:This is precession.

def gyroscope_dynamics():
    """Gyroscope dynamics"""
    
    def euler_equations(t, state, I1, I2, I3):
        omega1, omega2, omega3 = state
        domega1 = (I2 - I3) * omega2 * omega3 / I1
        domega2 = (I3 - I1) * omega3 * omega1 / I2
        domega3 = (I1 - I2) * omega1 * omega2 / I3
        return [domega1, domega2, domega3]
    
    fig = plt.figure(figsize=(14, 10))
    
    # Symmetric top (I1 = I2)
    ax1 = fig.add_subplot(2, 2, 1)
    
    I1, I2, I3 = 1.0, 1.0, 0.5  # Oblate top
    omega0 = [0.1, 0, 10]  # Mainly rotating about z-axis
    
    sol = solve_ivp(euler_equations, [0, 10], omega0, args=(I1, I2, I3),
                   t_eval=np.linspace(0, 10, 1000))
    
    ax1.plot(sol.t, sol.y[0], 'r-', linewidth=2, label='$\\omega_1$')
    ax1.plot(sol.t, sol.y[1], 'g-', linewidth=2, label='$\\omega_2$')
    ax1.plot(sol.t, sol.y[2], 'b-', linewidth=2, label='$\\omega_3$')
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('Angular velocity', fontsize=12)
    ax1.set_title('Symmetric Top (I ₁ = I ₂)', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 3D trajectory
    ax2 = fig.add_subplot(2, 2, 2, projection='3d')
    ax2.plot(sol.y[0], sol.y[1], sol.y[2], 'b-', linewidth=1)
    ax2.set_xlabel('$\\omega_1$')
    ax2.set_ylabel('$\\omega_2$')
    ax2.set_zlabel('$\\omega_3$')
    ax2.set_title('Angular Velocity Trajectory', fontsize=14, fontweight='bold')
    
    # Asymmetric case (Euler top)
    ax3 = fig.add_subplot(2, 2, 3)
    
    I1, I2, I3 = 1.0, 2.0, 3.0
    omega0 = [1, 0.1, 0.1]
    
    sol = solve_ivp(euler_equations, [0, 50], omega0, args=(I1, I2, I3),
                   t_eval=np.linspace(0, 50, 2000))
    
    ax3.plot(sol.t, sol.y[0], 'r-', linewidth=1, label='$\\omega_1$')
    ax3.plot(sol.t, sol.y[1], 'g-', linewidth=1, label='$\\omega_2$')
    ax3.plot(sol.t, sol.y[2], 'b-', linewidth=1, label='$\\omega_3$')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Angular velocity', fontsize=12)
    ax3.set_title('Asymmetric Top (I ₁ ≠ I ₂ ≠ I ₃)', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Energy and angular momentum conservation
    ax4 = fig.add_subplot(2, 2, 4)
    
    # Calculate energy and angular momentum
    E = 0.5 * (I1*sol.y[0]**2 + I2*sol.y[1]**2 + I3*sol.y[2]**2)
    L_sq = (I1*sol.y[0])**2 + (I2*sol.y[1])**2 + (I3*sol.y[2])**2
    
    ax4.plot(sol.t, E/E[0], 'b-', linewidth=2, label='Energy (normalized)')
    ax4.plot(sol.t, np.sqrt(L_sq)/np.sqrt(L_sq[0]), 'r--', linewidth=2, label='|L| (normalized)')
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('Conserved quantity', fontsize=12)
    ax4.set_title('Conservation Laws', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_ylim(0.99, 1.01)
    
    plt.tight_layout()
    plt.savefig('gyroscope_dynamics.png', dpi=150, bbox_inches='tight')
    plt.show()

gyroscope_dynamics()

Summary

In this chapter, we explored the broad applications of ordinary differential equations in physics and engineering:

Mechanical Systems

• Point mass motion: From simple free fall to complex planetary motion
• Vibration systems: Single and multi-degree-of-freedom, linear and nonlinear
• Rigid body dynamics: Euler equations and gyroscopic motion

Electromagnetism

• Circuit analysis: RLC circuits, coupled circuits, nonlinear elements
• Frequency response: Transfer functions and Bode plots

Thermal and Chemical

• Heat conduction: Newton's law of cooling and its applications
• Reaction kinetics: First-order, consecutive, and reversible reactions

Biological and Ecological

• Population dynamics: Predator-prey models, competitive exclusion
• Epidemic models: SIR model and basic reproduction number

Control Systems

• Feedback control: Stability and pole placement
• PID control: Most common control strategy in industry

Fluid Mechanics

• Simplified models: Poiseuille flow, boundary layer theory

Exercises

Basic Problems

1. Akg object falls from high altitude with air resistance. Find the terminal velocity and time to reach 99% of terminal velocity.

2. Solve the damped harmonic oscillatorwith initial conditions,.

3. In an RLC circuit withH,μ F,Ω, determine if the circuit is underdamped or overdamped and find the natural frequency.

4. A first-order reaction has a half-life of 2 hours. Find the time for the reactant concentration to drop to 10% of its initial value.

5. In the SIR model, if/day,/day, calculateand explain its physical meaning.

Advanced Problems

6. Double pendulum: Derive the equations of motion and solve numerically. Observe sensitivity to initial conditions.

7. Van der Pol oscillator:Analyze the limit cycle behavior for different values of.

8. Chemical oscillation (simplified Belousov-Zhabotinsky reaction model): Explore oscillatory behavior in parameter space.

9. Spacecraft attitude control: Design a simple feedback control law to stabilize satellite attitude angle to zero.

10. Epidemic control: Introduce vaccination ratein the SIR model and analyze how to chooseto prevent outbreak.

Programming Problems

11. Implement a general second-order system simulator that can analyze step response, frequency response, and phase plane trajectories for different damping ratios.

12. Write a program to solve the restricted three-body problem and visualize orbits.

13. Implement a PID controller for temperature control system simulation, including auto-tuning capability.

14. Simulate multi-degree-of-freedom building structure response under earthquake excitation, analyzing effects of different interstory stiffness.

Discussion Questions

15. Why can a gyroscope maintain stable orientation? Explain from the perspective of angular momentum conservation.

16. In control systems, why does pure integral control lead to instability or oscillation?

17. In ecosystems, why are the periodic oscillations predicted by the Lotka-Volterra model rarely observed in reality?

18. From physical intuition, explain why critical damping is "optimal"— neither oscillating nor taking longest to reach equilibrium.

References

1. Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos. CRC Press.
2. Ogata, K. (2010). Modern Control Engineering. Prentice Hall.
3. Murray, J. D. (2002). Mathematical Biology. Springer.
4. Thornton, S. T., & Marion, J. B. (2004). Classical Dynamics of Particles and Systems. Brooks/Cole.
5. Dorf, R. C., & Bishop, R. H. (2017). Modern Control Systems. Pearson.
6. Batchelor, G. K. (2000). An Introduction to Fluid Dynamics. Cambridge University Press.

---

Next Chapter: Chapter 18: Advanced Topics and Summary →

---

This article is Chapter 17 of the "Erta Ordinary Differential Equations" series.