Learning #32 – Differential Equations in Physics: Oscillators, Chaos, PDEs and Numerical Methods

A differential equation relates a function to its derivatives. The solutions describe how physical systems evolve: how a population grows, how an electrical signal propagates, how a planet orbits, how heat diffuses, how a pendulum swings. Every branch of physics is ultimately a collection of differential equations and systematic techniques for solving them. This post builds those techniques from first principles, culminating in a discussion of chaotic systems where long-term prediction becomes impossible despite fully deterministic equations.

1. First-Order ODEs: Decay, Growth and Phase Lines

A first-order ODE has the form dy/dt = f(t, y). The fundamental techniques for solving it depend on the structure of f.

Separation of Variables

If f(t, y) = g(t) h(y), the equation is separable and can be integrated directly: ∫dy/h(y) = ∫g(t) dt. Classic examples: radioactive decay (h(y)=y, g(t)=−λ), giving y(t)=y₀e^−λt; and the logistic equation.

dN/dt = r N (1 − N/K)

N = population,  r = intrinsic growth rate,  K = carrying capacity

Exact solution:  N(t) = K / (1 + ((K−N₀)/N₀) e^{−rt})

Fixed points: N* = 0 (unstable),  N* = K (stable)
  → determined from phase line of dN/dt = f(N)

Linear First-Order ODEs and the Integrating Factor

For dy/dt + p(t)y = q(t), multiply both sides by the integrating factor μ(t) = exp(∫p(t) dt) to make the left side an exact derivative: d(μy)/dt = μ q(t). This technique appears in series RC circuit discharge, Newton’s law of cooling, and the Langevin equation for Brownian motion.

2. Second-Order Linear ODEs: Oscillators and Resonance

The prototype is the damped harmonic oscillator: m&ddot;x + b˙x + kx = F(t). The homogeneous solution (F=0) depends on the discriminant of the characteristic polynomial mλ² + bλ + k = 0.

m x'' + b x' + k x = F cos(ωt)

ω₀ = √(k/m)     (natural frequency)
ζ  = b/(2mω₀)  (damping ratio)

Characteristic roots: λ = −ζω₀ ± ω₀√(ζ²−1)

Underdamped (ζ < 1): x(t) = A e^{−ζω₀t} cos(ω_d t + ϕ)
  ω_d = ω₀ √(1−ζ²)   (damped natural frequency)

Critically damped (ζ = 1): x(t) = (A + Bt) e^{−ω₀t}

Overdamped (ζ > 1): x(t) = A e^{λ₁t} + B e^{λ₂t},  λ₁,λ₂ < 0

Quality factor:  Q = ω₀/2ζω₀ = 1/(2ζ)
  Q ≫ 1  → resonant system (bell, LC circuit, optical cavity)
  Q ≈ 0.5  → critically damped door-closer

Driven Oscillator and Resonance

When F(t) = F₀cos(ωt), the steady-state particular solution has amplitude:

A(ω) = F₀/m / √((ω₀² − ω²)² + (2ζω₀ω)²)

Maximum amplitude at:  ω_res = ω₀ √(1 − 2ζ²)  (for ζ < 1/√2)

Peak amplitude:  A_max = F₀/(2mζω₀²√(1−ζ²)) ≈ F₀Q/(mω₀²)  (small ζ)

Phase lag:  φ(ω) = arctan(2ζωω₀ / (ω₀² − ω²))
  φ → 0    as ω → 0  (in phase with force)
  φ = π/2 at ω = ω₀  (quadrature at resonance)
  φ → π    as ω → ∞  (180° out of phase)

3. Systems of ODEs and Phase Plane Analysis

Many physical systems are described by autonomous systems of ODEs: dx/dt = f(x,y), dy/dt = g(x,y). The phase plane (x,y) visualises trajectories without solving explicitly. The key tool is linearisation around fixed points.

Fixed points (x*, y*): f(x*,y*) = g(x*,y*) = 0

Jacobian at (x*,y*):
  J = | ∂f/∂x   ∂f/∂y |
      | ∂g/∂x   ∂g/∂y |

Eigenvalues λ₁, λ₂ of J determine local stability:
  Re(λ) < 0 both: stable node/spiral (attractor)
  Re(λ) > 0 both: unstable node/spiral
  Re(λ) opposite signs: saddle point (unstable)
  Re(λ) = 0: centre (neutral, true nonlinear analysis needed)

Trace τ = λ₁+λ₂ = tr(J),   Determinant Δ = λ₁λ₂ = det(J)
Stable node:   τ < 0, Δ > 0, τ² > 4Δ
Stable spiral: τ < 0, Δ > 0, τ² < 4Δ

Lotka–Volterra Predator–Prey System

dx/dt = αx − βxy, dy/dt = δxy − γy (x = prey, y = predator). Two fixed points: (0,0) (saddle) and (γ/δ, α/β) (centre in the nonlinear model — closed orbits). The conserved quantity is V(x,y) = δx − γln x + βy − αln y.

4. Partial Differential Equations: Wave, Heat and Laplace

When a function depends on both space and time (or multiple space coordinates), its evolution is governed by a PDE. The three canonical second-order linear PDEs each represent a different physical archetype.

The Wave Equation

∂²u/∂t² = c² ∂²u/∂x²

D'Alembert solution:  u(x,t) = f(x − ct) + g(x + ct)
  f: right-travelling wave, g: left-travelling wave
  given u(x,0) = φ(x), u_t(x,0) = ψ(x):
  f(ξ) = [φ(ξ) − Κ(ξ)]/2,  g(ξ) = [φ(ξ) + Κ(ξ)]/2
  where Κ(x) = (1/c) ∫ ψ dx

Standing wave on string [0,L] with u(0,t)=u(L,t)=0:
  u_n(x,t) = sin(nπx/L) [A_n cos(ω_n t) + B_n sin(ω_n t)]
  ω_n = nπc/L  (normal mode frequencies)

The Heat Equation

∂u/∂t = α ∂²u/∂x²       α = thermal diffusivity [m²/s]

Fourier series solution on [0,L] with u(0,t)=u(L,t)=0:
  u(x,t) = ∑_n B_n sin(nπx/L) exp(−α(nπ/L)² t)

B_n = (2/L) ∫₀ᴸ u(x,0) sin(nπx/L) dx

Characteristic diffusion length:  ℓ ≈ √(αt)
  Thermal diffusivity: Cu ≈ 1.17×10⁻⁴ m²/s, Si ≈ 8×10⁻⁵ m²/s

Laplace’s Equation

∇²u = 0 governs electrostatic potential in free space, steady-state temperature distributions, incompressible irrotational fluid flow, and complex-analytic functions. Solutions (harmonic functions) satisfy the mean-value property and have no local extrema in the interior of the domain. Separation of variables in Cartesian, polar, and spherical coordinates gives families of solutions (Legendre polynomials, spherical harmonics) that are the basis of multipole expansions in electrostatics and gravity.

5. Numerical Methods: Euler and Runge–Kutta

Most ODEs encountered in practice cannot be solved analytically. Numerical integration methods discretise time and advance the solution step by step. The key trade-off is between accuracy (local truncation error) and computational cost.

Given dy/dt = f(t,y),  step size h,  current (t_n, y_n):

k₁ = f(t_n,          y_n)
k₂ = f(t_n + h/2,   y_n + h k₁/2)
k₃ = f(t_n + h/2,   y_n + h k₂/2)
k₄ = f(t_n + h,     y_n + h k₃)

y_{n+1} = y_n + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)

Local truncation error: O(h⁵)   Global error: O(h⁴)

Comparison:
  Euler:        LTE O(h²), GE O(h)   [1 f-eval/step]
  Trapezoidal:  LTE O(h³), GE O(h²)  [2 f-evals/step]
  RK4:          LTE O(h⁵), GE O(h⁴)  [4 f-evals/step]
  Dormand-Prince RK45: adaptive, used by scipy.integrate.solve_ivp

Stiff ODEs (where some time scales are much faster than others) require implicit methods (e.g., backward Euler, Crank–Nicolson) or implicit Runge–Kutta to avoid requiring extremely small step sizes. A classic stiff example is the Van der Pol oscillator with large damping parameter μ.

6. Chaos and Nonlinear Dynamics

Deterministic chaos is not randomness — the equations are exact and perfectly deterministic. What makes a system chaotic is sensitive dependence on initial conditions: nearby trajectories diverge exponentially fast, making long-term prediction impossible in practice despite being possible in principle.

dx/dt = σ(y − x)         σ = 10
dy/dt = x(ρ − z) − y   ρ = 28
dz/dt = xy − βz          β = 8/3

Max Lyapunov exponent λ₁ ≈ +0.906 (positive ⇒ chaotic)
  ⇒ separation: |Δ(t)| ~ |Δ(0)| e^{λ₁ t}
  Predictability horizon: T ≈ ln(Δ_max/|Δ(0)|) / λ₁
  With |Δ(0)|=10⁻10, Δ_max=1: T ≈ 23/0.906 ≈ 25 Lorenz time-units

Attractor dimension (Kaplan-Yorke): D_KY ≈ 2.06

Bifurcation Diagrams and the Feigenbaum Constant

As a control parameter is varied, a nonlinear system can undergo sudden qualitative changes in behaviour (bifurcations). The logistic map x_n+1 = r x_n(1−x_n) shows a period-doubling cascade as r increases from 3 toward 4: stable fixed point → period-2 cycle → period-4 → … → chaos at r≈3.57. The ratio of successive bifurcation intervals converges to the Feigenbaum constant δ ≈ 4.6692… — universal for all one-dimensional maps with a quadratic maximum (period-doubling universality).

Differential Equations as Physical Intuition

The progression in this post — from separation of variables to RK4 integration to the Lorenz attractor — mirrors the historical development of dynamics over three centuries. Newton’s second law launched the programme of writing physics as differential equations. Fourier’s heat equation and the development of spectral methods expanded it to distributed systems. Poincaré’s geometric approach to differential equations (1880s) introduced phase planes, fixed points, and the first glimpses of chaotic behaviour.

The interactive simulations let you adjust damping and driving frequency to trace the resonance curve of a harmonic oscillator in real time, watch a predator–prey limit cycle in phase space, and observe the Lorenz attractor’s sensitive dependence on initial conditions by launching two nearly identical trajectories and watching them diverge. These are not just computational exercises; they are the same physical phenomena that engineers, biologists, climatologists, and physicists encounter daily.