Systems of ODEs

⚡ Key Takeaways

A system of ordinary differential equations (ODEs) is a collection of coupled equations relating several unknown functions of a single variable and their derivatives. In words: instead of one quantity evolving in time, several quantities evolve simultaneously, each influencing the others. They appear everywhere — predator–prey dynamics, electrical circuits, chemical kinetics, mechanics, control engineering. Theory crystallized in the 18th–19th centuries (Euler, Lagrange, Cauchy, Poincaré) and matured into the modern dynamical-systems viewpoint in the 20th.

💡 The Core Insight

Every ODE problem is secretly a first-order vector field. Once you bundle all unknowns and their derivatives into a single state vector $\mathbf{x}$ , the entire dynamics is encoded by one map $\mathbf{f}: \mathbb{R}^{n+1} \to \mathbb{R}^n$ that assigns a velocity to every point of state space. Solving the system means finding the flow — the family of curves that are everywhere tangent to this vector field. Geometry replaces formula-hunting: even when no closed-form solution exists, the structure of equilibria, eigenvalues, and trajectories tells you how the system behaves.

DefinitionSystem of First-Order ODEs

A first-order system of $n$ ordinary differential equations is

1

\quad \mathbf{x}'(t) = \mathbf{f}(t, \mathbf{x}(t)), \qquad \mathbf{x}(t)=\begin{pmatrix} x_1(t)\\ \vdots \\ x_n(t)\end{pmatrix},\quad \mathbf{f}=\begin{pmatrix} f_1\\ \vdots \\ f_n\end{pmatrix}.

Derivation hint: Stack the $n$ scalar equations $x_i'=f_i(t,x_1,\dots,x_n)$ into a single vector equation; this is purely notational bundling. Here $t$ is the independent variable (often time), $\mathbf{x}(t)$ is the state vector living in state space (or phase space) $\mathbb{R}^n$ , and $\mathbf{f}$ is the vector field. An initial value problem (IVP) adjoins the condition $\mathbf{x}(t_0)=\mathbf{x}_0$ . The system is autonomous if $\mathbf{f}$ does not depend explicitly on $t$ , i.e. $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ . It is linear if $\mathbf{f}(t,\mathbf{x})=A(t)\mathbf{x}+\mathbf{g}(t)$ for a matrix $A(t)$ and vector $\mathbf{g}(t)$ ; homogeneous if additionally $\mathbf{g}\equiv\mathbf{0}$ .

💡 The Core Insight

Any single $n$ -th order ODE $y^{(n)}=g(t,y,y',\dots,y^{(n-1)})$ becomes a first-order system. Set $x_1=y,\ x_2=y',\ \dots,\ x_n=y^{(n-1)}$ . Then

2

\quad x_1'=x_2,\quad x_2'=x_3,\quad\dots,\quad x_{n-1}'=x_n,\quad x_n'=g(t,x_1,\dots,x_n).

Derivation hint: Each new variable is the next derivative of $y$ ; the chain of definitions makes the first $n-1$ equations trivial, and the original equation supplies the last.

DefinitionLinear Homogeneous System & Matrix Exponential

The system $\mathbf{x}'=A\mathbf{x}$ with constant $n\times n$ matrix $A$ has the unique solution

3

\quad \mathbf{x}(t)=e^{A(t-t_0)}\,\mathbf{x}_0, \qquad e^{At}:=\sum_{k=0}^{\infty}\frac{(At)^k}{k!}.

Derivation hint: Differentiate the power series term by term, $\frac{d}{dt}e^{At}=A e^{At}$ , and verify it satisfies the equation with $e^{A\cdot 0}=I$ ; uniqueness follows from the existence–uniqueness theorem.

TheoremExistence and Uniqueness (Picard–Lindelöf)

Let $\mathbf{f}(t,\mathbf{x})$ be continuous on an open region $D\subseteq \mathbb{R}\times\mathbb{R}^n$ containing $(t_0,\mathbf{x}_0)$ , and suppose $\mathbf{f}$ is Lipschitz in $\mathbf{x}$ : there exists $L>0$ with

4

\quad \|\mathbf{f}(t,\mathbf{x})-\mathbf{f}(t,\mathbf{y})\|\le L\,\|\mathbf{x}-\mathbf{y}\| \quad \text{for all } (t,\mathbf{x}),(t,\mathbf{y})\in D.

Derivation hint: The bound controls how fast nearby trajectories can separate; it is exactly what makes the Picard iteration a contraction. Then the IVP $\mathbf{x}'=\mathbf{f}(t,\mathbf{x}),\ \mathbf{x}(t_0)=\mathbf{x}_0$ has a unique solution on some interval $(t_0-h,t_0+h)$ .

TheoremEigenvalue Solution of Linear Systems

If $A$ is diagonalizable with eigenpairs $(\lambda_i,\mathbf{v}_i)$ , $i=1,\dots,n$ , the general solution of $\mathbf{x}'=A\mathbf{x}$ is

5

\quad \mathbf{x}(t)=\sum_{i=1}^{n} c_i\, e^{\lambda_i t}\,\mathbf{v}_i, \qquad c_i\in\mathbb{C}.

Derivation hint: Substitute the trial solution $\mathbf{x}=e^{\lambda t}\mathbf{v}$ into $\mathbf{x}'=A\mathbf{x}$ to get $(A-\lambda I)\mathbf{v}=\mathbf{0}$ ; superposition of these independent modes gives the general solution. For complex conjugate eigenvalues $\lambda=\alpha\pm i\beta$ , real solutions are extracted via Euler's formula, yielding $e^{\alpha t}(\cos\beta t,\ \sin\beta t)$ combinations. For a repeated eigenvalue lacking enough eigenvectors, generalized eigenvectors introduce polynomial-in- $t$ terms.

ProofPicard Iteration Proves Existence and Uniqueness

We sketch the standard proof of Theorem [4]. 1. Reformulate as an integral equation. A continuous $\mathbf{x}(t)$ solves the IVP iff

6

\quad \mathbf{x}(t)=\mathbf{x}_0+\int_{t_0}^{t}\mathbf{f}(s,\mathbf{x}(s))\,ds.

Derivation hint: Integrate $\mathbf{x}'=\mathbf{f}$ from $t_0$ to $t$ and apply the Fundamental Theorem of Calculus; the constant of integration is fixed by $\mathbf{x}(t_0)=\mathbf{x}_0$ . 2. Define the Picard operator $(\mathcal{T}\mathbf{x})(t)=\mathbf{x}_0+\int_{t_0}^t \mathbf{f}(s,\mathbf{x}(s))\,ds$ acting on continuous functions on $[t_0-h,t_0+h]$ valued in a closed ball about $\mathbf{x}_0$ . 3. Show $\mathcal{T}$ is a contraction for small $h$ . Using the Lipschitz bound,

7

\quad \|(\mathcal{T}\mathbf{x})(t)-(\mathcal{T}\mathbf{y})(t)\|\le L\int_{t_0}^{t}\|\mathbf{x}(s)-\mathbf{y}(s)\|\,ds \le Lh\,\sup_{s}\|\mathbf{x}(s)-\mathbf{y}(s)\|.

Derivation hint: Subtract the two integral forms, pull the norm inside the integral, and apply the Lipschitz inequality [4]; choose $h<1/L$ so $Lh<1$ . 4. Apply the Banach fixed-point theorem. Since $\mathcal{T}$ contracts a complete metric space (continuous functions with sup norm), it has a unique fixed point $\mathbf{x}^*$ , which is exactly the unique solution of the IVP. $\blacksquare$

∎

ProofWhy the Matrix Exponential Solves the Linear System

1. Define $\Phi(t)=e^{At}$ via the absolutely convergent series in [3]. The series converges for every $t$ because $\|(At)^k/k!\|\le (\|A\|\,|t|)^k/k!$ and $\sum$ gives $e^{\|A\||t|}$ . 2. Differentiate term by term (legal due to uniform convergence on compacts):

8

\quad \frac{d}{dt}e^{At}=\sum_{k=1}^{\infty}\frac{A^k t^{k-1}}{(k-1)!}=A\sum_{j=0}^{\infty}\frac{(At)^j}{j!}=A\,e^{At}.

Derivation hint: Shift the index $k\to j=k-1$ and factor one $A$ out front. 3. Thus $\mathbf{x}(t)=e^{At}\mathbf{x}_0$ satisfies $\mathbf{x}'=A\mathbf{x}$ and $\mathbf{x}(0)=\mathbf{x}_0$ . By Picard–Lindelöf (the field $A\mathbf{x}$ is globally Lipschitz with $L=\|A\|$ ) it is the unique solution. $\blacksquare$

∎

mathematician · Henri Poincaré | 1854–1912

Poincaré transformed the subject from a quest for explicit formulas into the qualitative theory of dynamical systems. Frustrated that most systems (notably the three-body problem) cannot be solved in closed form, he asked instead: what does the family of all solutions look like? He invented phase-space portraits, the notion of recurrence, Poincaré sections, index theory for vector fields, and the first glimpse of chaos in his 1890 study of celestial mechanics. His maxim — that geometry and topology of trajectories matter more than formulas — underpins everything in the modern treatment of systems of ODEs.

🧠 Intuition

Picture state space $\mathbb{R}^n$ with an arrow attached to every point: that arrow is $\mathbf{f}(\mathbf{x})$ , the instantaneous velocity. A solution is a particle released into this wind, drifting along so its velocity always matches the local arrow. The collection of all such streamlines is the phase portrait.

Equilibria ( $\mathbf{f}(\mathbf{x}^*)=\mathbf{0}$ ) are points of dead calm: a particle placed exactly there never moves.
Eigenvalues of the linearization at an equilibrium tell you the local weather: real negative ⇒ trajectories converge (stable node), real positive ⇒ they flee (unstable node), opposite signs ⇒ saddle, complex with negative real part ⇒ inward spiral, pure imaginary ⇒ closed loops (center).
For a 2D linear system, the trace and determinant of $A$ alone classify the portrait: $\det A<0$ gives a saddle; $\det A>0$ with $\operatorname{tr}A$ deciding spiral/node and stability.

The big mental shift: you don't need the formula to understand the motion. You read the flow off the field.

Figure

Phase portrait of a stable spiral: trajectory winding into the origin under a 2×2 linear system with complex eigenvalues (negative real part).

ExampleCoupled Decay / Mixing Tanks (Real Eigenvalues)

Given. Two brine tanks. Let $x_1,x_2$ be salt (kg) in tanks 1 and 2. Flow gives

x_1'=-2x_1+x_2,\qquad x_2'=2x_1-3x_2,\qquad x_1(0)=10,\ x_2(0)=0.

So $A=\begin{pmatrix}-2 & 1\\ 2 & -3\end{pmatrix}$ .

Step 1 — Eigenvalues. $\det(A-\lambda I)=(-2-\lambda)(-3-\lambda)-2=\lambda^2+5\lambda+4=(\lambda+1)(\lambda+4)$ , so $\lambda_1=-1,\ \lambda_2=-4$ .

Step 2 — Eigenvectors. For $\lambda_1=-1$ : $(A+I)\mathbf v=\begin{pmatrix}-1&1\\2&-2\end{pmatrix}\mathbf v=0\Rightarrow \mathbf v_1=(1,1)^T$ . For $\lambda_2=-4$ : $(A+4I)\mathbf v=\begin{pmatrix}2&1\\2&1\end{pmatrix}\mathbf v=0\Rightarrow \mathbf v_2=(1,-2)^T$ .

Step 3 — General solution (eq. [5]).

\mathbf x(t)=c_1 e^{-t}\binom{1}{1}+c_2 e^{-4t}\binom{1}{-2}.

Step 4 — Apply initial conditions. At $t=0$ : $c_1+c_2=10,\ c_1-2c_2=0\Rightarrow c_1=\tfrac{20}{3},\ c_2=\tfrac{10}{3}.$

Boxed answer.

\boxed{\,x_1(t)=\tfrac{20}{3}e^{-t}+\tfrac{10}{3}e^{-4t}\ \text{kg},\qquad x_2(t)=\tfrac{20}{3}e^{-t}-\tfrac{20}{3}e^{-4t}\ \text{kg}\,}

Both eigenvalues negative ⇒ the origin is a stable node; all salt eventually drains, $\mathbf x(t)\to(0,0)$ as $t\to\infty$ .

ExampleOscillator as a System (Complex Eigenvalues)

Given. Mass–spring with light damping: $y''+2y'+5y=0$ , $y(0)=1,\ y'(0)=0$ (units: $y$ in m, $t$ in s).

Step 1 — Convert to a system (eq. [2]). Let $x_1=y,\ x_2=y'$ :

x_1'=x_2,\qquad x_2'=-5x_1-2x_2,\qquad A=\begin{pmatrix}0&1\\-5&-2\end{pmatrix}.

Step 2 — Eigenvalues. $\lambda^2+2\lambda+5=0\Rightarrow \lambda=-1\pm 2i$ . So $\alpha=-1$ , $\beta=2$ .

Step 3 — Real solution form. Complex eigenvalues with negative real part ⇒ a decaying spiral:

y(t)=e^{-t}\big(C_1\cos 2t + C_2\sin 2t\big).

Step 4 — Initial conditions. $y(0)=C_1=1$ . Then $y'(t)=e^{-t}(-C_1\cos2t-C_2\sin2t-2C_1\sin2t+2C_2\cos2t)$ ; at $t=0$ : $-C_1+2C_2=0\Rightarrow C_2=\tfrac12$ .

Boxed answer.

\boxed{\,y(t)=e^{-t}\!\left(\cos 2t+\tfrac12\sin 2t\right)\ \text{m}\,}

The phase portrait is a stable focus (cf. the SVG above): oscillations of angular frequency $2\ \text{rad/s}$ with amplitude decaying like $e^{-t}$ .

⚠️ Common mistake

Eigenvalues are not just a trick for linear systems — but they only describe the linearization of nonlinear ones. For $\mathbf{x}'=\mathbf{f}(\mathbf{x})$ , you compute eigenvalues of the Jacobian $D\mathbf{f}(\mathbf{x}^*)$ at an equilibrium. The Hartman–Grobman theorem guarantees the local portrait matches the linear one only when the equilibrium is hyperbolic</strong> (no eigenvalue has zero real part). For a center (pure imaginary eigenvalues), linearization is silent: the true nonlinear flow could spiral in, spiral out, or form closed orbits. Other frequent errors:

Forgetting that a repeated eigenvalue may be defective (too few eigenvectors), requiring generalized eigenvectors and $t\,e^{\lambda t}$ terms — not simply two copies of $e^{\lambda t}\mathbf v$ .
Assuming solutions exist for all time. The Lipschitz condition gives only a local interval; nonlinear systems can blow up in finite time.
Treating $e^{A}e^{B}=e^{A+B}$ as always valid — it holds only if $A$ and $B$ commute ( $AB=BA$ ).

CounterexampleA Center That Lies, and a Finite-Time Blowup

Center failure. Consider $x'=-y+x(x^2+y^2),\ y'=x+y(x^2+y^2)$ . The Jacobian at the origin is $\begin{pmatrix}0&-1\\1&0\end{pmatrix}$ with eigenvalues $\pm i$ — a linear center. But in polar coordinates $r'=r^3>0$ : every trajectory spirals outward. The linearization predicts closed loops; the true flow is an unstable spiral. Hyperbolicity fails, so linearization gives the wrong answer.

Finite-time blowup. The scalar (hence 1D system) IVP $x'=x^2,\ x(0)=1$ is perfectly smooth and Lipschitz locally, yet its solution $x(t)=\dfrac{1}{1-t}$ escapes to $+\infty$ as $t\to 1^-$ . No global solution exists despite the field being analytic everywhere.

connection

Systems of ODEs sit at a crossroads of analysis, algebra, and geometry.

They are powered by Linear Algebra: eigenvalues, eigenvectors, diagonalization, and Jordan form are the entire toolkit for constant-coefficient linear systems.
The matrix exponential [3] is a special case of the Matrix Exponential and connects to operator theory and Lie groups.
Qualitative analysis flows into Dynamical Systems and Stability Theory (Lyapunov functions, limit cycles, bifurcations).
The continuous-time picture has a discrete cousin in Difference Equations, where $e^{At}$ is replaced by powers of a matrix.
When the independent variable becomes space and time, you generalize to Partial Differential Equations; method-of-lines discretization turns a PDE back into a large ODE system.
Numerical solution links to Numerical Methods for ODEs (Euler, Runge–Kutta, stiff solvers).

🔬 Where this appears

⏳ History

- 1736–1755 — Leonhard Euler systematizes solution methods for linear ODEs and the exponential ansatz

e^{\lambda t}

, foreshadowing eigenvalue methods.

- 1762–1788 — Joseph-Louis Lagrange develops variation of parameters and treats coupled equations arising in mechanics; his Mécanique analytique (1788) frames dynamics as systems of equations.

- 1788–1833 — Hamilton & Jacobi recast mechanics as first-order systems in phase space (Hamilton's equations), the prototype of state-space thinking.

- 1820s–1840s — Augustin-Louis Cauchy establishes the first rigorous existence proofs (the Cauchy–Lipschitz method), giving the IVP a sound foundation.

- 1877 — Rudolf Lipschitz sharpens the existence–uniqueness condition; 1890s — Émile Picard & Ernst Lindelöf give the iteration/fixed-point proof (Theorem [4]).

- 1881–1886 — Henri Poincaré founds the qualitative theory: phase portraits, recurrence, index theory, and the discovery of chaotic behavior in the three-body problem.

- 1892 — Aleksandr Lyapunov publishes The General Problem of the Stability of Motion, introducing Lyapunov functions and the modern definition of stability.

- 1927–1928 — van der Pol & van der Mark study self-sustained electronic oscillations, popularizing limit cycles and relaxation oscillations in coupled circuits.

- 1963 — Edward Lorenz numerically discovers sensitive dependence on initial conditions in a three-variable ODE system, launching chaos theory.

open-problem

The deepest difficulties live in the nonlinear, non-hyperbolic, and chaotic regimes.

No general closed-form solutions. Most nonlinear systems (the three-body problem being the historic example) have no explicit solution in elementary functions; Poincaré proved no general analytic first integrals exist for the restricted three-body problem.
Chaos and predictability. Even fully deterministic systems like Lorenz's exhibit exponential sensitivity (positive Lyapunov exponents), so long-term prediction is fundamentally limited despite the equations being known exactly.
Bifurcation theory classifies how phase portraits change as parameters vary (saddle-node, Hopf, period-doubling cascades), but a complete classification in dimensions $\ge 3$ is open.
Global existence vs. blowup. Determining whether a given nonlinear system has globally defined solutions can be subtle; for PDE analogues (Navier–Stokes) it is a Millennium Prize Problem, and finite-dimensional analogues feed intuition.
Structural stability and the Poincaré–Bendixson boundary. In the plane, the Poincaré–Bendixson theorem forbids chaos (trajectories must approach equilibria or limit cycles); in three or more dimensions this guardrail vanishes, and the full taxonomy of attractors remains an active research frontier.

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also