Computational Dynamics: Numerical Methods

An Earth climate model, a system of 10⁸ variables, integrated 100 years forward. Molecular dynamics of a protein, 10⁵ atoms, 10⁹ steps. A Pluto flyby, 9 years of flight, accuracy ±1 km. All of this relies on numerical ODE methods maintaining accuracy where analysis is impossible.

• **Meteorology:** ECMWF (European Centre for Medium-Range Weather Forecasts) uses spectral methods (Fourier-space analog of RK4) to integrate Navier-Stokes equations 10 days forward
• **Neural modeling:** NEURON/Brian2 use adaptive Runge-Kutta to simulate ~10⁴ neurons in real time. CVODE (LLNL), the standard for stiff neuron models
• **Pharmacology:** PK/PD ODE systems with AUTO continuation, finding dose-response bifurcations to predict side effects

Предварительные знания

• Continuous Dynamical Systems
• Bifurcations

Numerical ODE Integrators

SpaceX Dragon's landing trajectory is computed numerically at 100Hz , RK4 integration achieves 10⁻⁸ m accuracy per step, making real-time autonomous landing possible. Almost all real-world ODEs have no analytical solutions. Numerical methods are the workhorse of computational dynamics. Core task: given ẋ = f(x, t), x(t₀) = x₀, find x(T) to a specified accuracy ε.

For smooth non-stiff systems: RK45 (Dormand-Prince), the gold standard. For Hamiltonian systems: symplectic methods (Verlet, leapfrog), structurally preserve energy. For stiff systems: implicit methods (Radau, BDF), explicit methods diverge.

Stiff Systems and Adaptive Methods

**Stiff ODEs** are systems with very disparate time scales: some variables change in nanoseconds, others in seconds. For explicit methods an extremely small step h ≪ min(1/|λ_i|) is required, this is computationally wasteful.

Stiffness indicator: if explicit RK45 constantly reduces step size (rejected steps >> accepted steps), the system is likely stiff. **LSODA** (FORTRAN classic, 1983) automatically switches between explicit Adams methods and implicit BDF, the de-facto standard in scipy.

Numerical Continuation and AUTO

**Continuation** is a powerful technique: instead of finding equilibria at one parameter value r, we 'follow' them as r changes gradually. This allows automatically constructing **bifurcation diagrams**.

AUTO-07p is the gold standard for bifurcation analysis. Alternatives: MATCONT (MATLAB, GUI), pde2path (PDEs), PyDSTool (Python). For neural modeling: XPPAUT (Ermentrout), ODE solver + AUTO with a graphical interface.

Numerical Lyapunov Exponents and Diagrams

Numerical computation of **Lyapunov exponents** is the standard tool for diagnosing chaos. The main challenge: variational vectors grow exponentially, overflowing floating-point range. Solution, periodic orthonormalization (Benettin's method).

Key Ideas

• **RK4:** gold standard for smooth non-stiff systems. Error O(h⁴), halving h gives 16× accuracy gain
• **Stiffness:** |λ_max|/|λ_min| ≫ 1 → explicit methods are inefficient. Radau/BDF, implicit methods with unconditional stability
• **Pseudo-arclength continuation (AUTO):** follow equilibrium branches in parameter space, round fold points, auto-detect bifurcations
• **Lyapunov exponents (Benettin):** QR orthonormalization prevents overflow. Λ₁ > 0 → chaos

Related Topics

Computational dynamics is the practical foundation for the entire course:

• Control Theory — Numerical solution of the Riccati equation (ARE), discretization for digital control, ODE solvers in MPC
• Chaos and Strange Attractors — Numerical Lyapunov exponents are the primary tool for diagnosing chaos. Benettin's method has been standard since 1980
• Bifurcations — AUTO/MATCONT are the computational implementation of bifurcation theory. Continuation produces complete bifurcation diagrams for nonlinear systems

Вопросы для размышления

• Adaptive ODE solvers automatically choose step size h. When might this be dangerous? For instance, in a real-time control system with hard timing constraints?
• Continuation methods like AUTO find all branches of equilibria. But how do one detect chaotic attractors, which are neither fixed points nor limit cycles?
• Numerical methods introduce discretization error. For chaotic systems this means the numerical trajectory diverges from the true one after O(log(1/ε)/λ₁) time. How does this limit climate forecasting?

Связанные уроки

• nm-01