Systems of Ordinary Differential Equations and Stability Theory

Systems of Ordinary Differential Equations (ODEs) are fundamental in modeling complex phenomena where multiple interdependent quantities evolve simultaneously. In this lesson, we transition from scalar equations to vector-valued differential equations, focusing on linear systems, the matrix exponential, and the qualitative behavior of solutions.

1. Linear Systems of First-Order ODEs

A general system of $n$ first-order linear ODEs can be expressed in vector-matrix form: ${\mathbf{x}}^{\prime }(t)=\mathbf{A}(t)\mathbf{x}(t)+\mathbf{g}(t)$ where $\mathbf{x}(t)\in {\mathbb{R}}^n$ is the state vector, $\mathbf{A}(t)\in {\mathbb{R}}^{n\times n}$ is the coefficient matrix, and $\mathbf{g}(t)$ is the non-homogeneous term.

For a system with constant coefficients ($\mathbf{A}$ is independent of $t$), the homogeneous system is: ${\mathbf{x}}^{\prime }=\mathbf{Ax}$

Existence and Uniqueness

By the Picard-Lindelöf Theorem, if $\mathbf{A}(t)$ and $\mathbf{g}(t)$ are continuous on an interval $I$, then for any $t_0\in I$ and ${\mathbf{x}}_0\in {\mathbb{R}}^n$, there exists a unique solution $\mathbf{x}(t)$ existing for all $t\in I$ satisfying $\mathbf{x}(t_0)={\mathbf{x}}_0$.

2. The Matrix Exponential

The solution to the scalar equation $x^{\prime }=ax$ is $x(t)=e^{at}x(0)$. By analogy, the solution to the vector equation ${\mathbf{x}}^{\prime }=\mathbf{Ax}$ is: $\mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0)$

Definition

The matrix exponential $e^{\mathbf{A}t}$ is defined by the power series: $e^{\mathbf{A}t}={\sum }_{k=0}^{\infty }\frac{(\mathbf{A}t{)}^k}{k!}=\mathbf{I}+\mathbf{A}t+\frac{{\mathbf{A}}^2{t}^2}{2!}+\dots$ This series converges absolutely for all $t$ and all square matrices $\mathbf{A}$.

Properties

1. $\frac{d}{dt}{e}^{\mathbf{A}t}=\mathbf{A}{e}^{\mathbf{A}t}=e^{\mathbf{A}t}\mathbf{A}$.
2. $e^{\mathbf{A}(t+s)}=e^{\mathbf{A}t}{e}^{\mathbf{A}s}$ if and only if any matrices involved commute.
3. $(e^{\mathbf{A}t}{)}^{-1}=e^{-\mathbf{A}t}$.
4. If $\mathbf{A}$ is diagonalizable, $\mathbf{A}=\mathbf{PD}{\mathbf{P}}^{-1}$, then $e^{\mathbf{A}t}=\mathbf{P}{e}^{\mathbf{D}t}{\mathbf{P}}^{-1}$, where $e^{\mathbf{D}t}=\text{diag}(e^{{\lambda }_{1}t},\dots ,e^{{\lambda }_{n}t})$.

3. Solving Homogeneous Systems

The general solution is a linear combination of linearly independent solutions $\mathbf{x}(t)={\sum }_{i=1}^n{c}_i{\mathbf{x}}_i(t)$. These solutions are derived from the eigenvalues ${\lambda }_i$ and eigenvectors ${\mathbf{v}}_i$ of $\mathbf{A}$.

Case 1: Distinct Real Eigenvalues

If $\mathbf{A}$ has $n$ linearly independent eigenvectors ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_n$ corresponding to real ${\lambda }_1,\dots ,{\lambda }_n$: $\mathbf{x}(t)=c_1{e}^{{\lambda }_{1}t}{\mathbf{v}}_1+\cdots +c_n{e}^{{\lambda }_{n}t}{\mathbf{v}}_n$

Case 2: Complex Eigenvalues

If $\lambda =\alpha \pm i\beta$ with eigenvectors $\mathbf{v}=\mathbf{u}\pm i\mathbf{w}$, the real-valued solutions are: ${\mathbf{x}}_1(t)=e^{\alpha t}(\cos (\beta t)\mathbf{u}-\sin (\beta t)\mathbf{w})$ ${\mathbf{x}}_2(t)=e^{\alpha t}(\sin (\beta t)\mathbf{u}+\cos (\beta t)\mathbf{w})$

Case 3: Repeated Eigenvalues

If an eigenvalue $\lambda$ has algebraic multiplicity greater than its geometric multiplicity, generalized eigenvectors are used. For a matrix with a single Jordan block of size 2: ${\mathbf{x}}_1(t)=e^{\lambda t}\mathbf{u},{\mathbf{x}}_2(t)=e^{\lambda t}(t\mathbf{u}+\mathbf{v})$ where $(\mathbf{A}-\lambda \mathbf{I})\mathbf{v}=\mathbf{u}$.

4. Phase Plane Analysis (2D Systems)

For ${\mathbf{x}}^{\prime }=\mathbf{Ax}$ where $\mathbf{x}\in {\mathbb{R}}^2$, the equilibrium point at the origin is classified by the eigenvalues:

• Nodes: ${\lambda }_1,{\lambda }_2$ real and same sign.
  • Stable Node: ${\lambda }_1,{\lambda }_2<0$.
  • Unstable Node: ${\lambda }_1,{\lambda }_2>0$.
• Saddle Points: ${\lambda }_1,{\lambda }_2$ real and opposite sign. Always unstable.
• Spiral Points: $\lambda =\alpha \pm i\beta$ with $\alpha \ne 0$.
  • Stable Spiral: $\alpha <0$.
  • Unstable Spiral: $\alpha >0$.
• Centers: $\lambda =\pm i\beta$ (purely imaginary). Marginally stable.

5. Stability Theory and Linearization

Stability Definitions

An equilibrium point ${\mathbf{x}}^{\ast }$ is:

• Stable (Lyapunov): If for every $ϵ>0$, there exists $\delta >0$ such that if $\Vert \mathbf{x}(0)-{\mathbf{x}}^{\ast }\Vert <\delta$, then $\Vert \mathbf{x}(t)-{\mathbf{x}}^{\ast }\Vert <ϵ$ for all $t\ge 0$.
• Asymptotically Stable: If it is stable and ${\lim }_{t\to \infty }\mathbf{x}(t)={\mathbf{x}}^{\ast }$.
• Unstable: If it is not stable.

Linearization and Hartman-Grobman

For a nonlinear system ${\mathbf{x}}^{\prime }=\mathbf{f}(\mathbf{x})$, let ${\mathbf{x}}^{\ast }$ be an equilibrium point. The Hartman-Grobman Theorem states that if ${\mathbf{x}}^{\ast }$ is hyperbolic (no eigenvalues of $Df({\mathbf{x}}^{\ast })$ have zero real part), then the nonlinear flow is topologically conjugate to the linear flow ${\mathbf{y}}^{\prime }=Df({\mathbf{x}}^{\ast })\mathbf{y}$ near the equilibrium.

Lyapunov’s Direct Method

Define a scalar function $V(\mathbf{x})$ such that $V({\mathbf{x}}^{\ast })=0$ and $V(\mathbf{x})>0$ for $\mathbf{x}\ne {\mathbf{x}}^{\ast }$. If $\dot{V}(\mathbf{x})=\nabla V\cdot \mathbf{f}(\mathbf{x})\le 0$, the equilibrium is stable. If $\dot{V}(\mathbf{x})<0$, it is asymptotically stable.

6. Non-Homogeneous Systems

The general solution to ${\mathbf{x}}^{\prime }=\mathbf{Ax}+\mathbf{g}(t)$ is given by the Variation of Parameters formula: $\mathbf{x}(t)=e^{\mathbf{A}t}{\mathbf{x}}_0+{\int }_0^t{e}^{\mathbf{A}(t-s)}\mathbf{g}(s)ds$

7. Application: Lotka-Volterra Predator-Prey Model

The interaction between prey $x$ and predators $y$ is modeled by: $\frac{dx}{dt}=\alpha x-\beta xy,\frac{dy}{dt}=\delta xy-\gamma y$ Linearizing around the interior equilibrium point $(\gamma /\delta ,\alpha /\beta )$ yields: $\mathbf{J}=\left(\begin{array}{cc}0& -\beta \gamma /\delta \\ \delta \alpha /\beta & 0\end{array}\right)$ The eigenvalues are purely imaginary ($\lambda =\pm i\sqrt{\alpha \gamma }$), indicating a center in the linearized system and periodic orbits in the nonlinear system.

Computational Example: Phase Portraits in Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

# Define the system matrix A for a stable spiral
A = np.array([[-0.5, 1.0], 
              [-1.0, -0.5]])

# 1. Compute Matrix Exponential for t=1.0
t = 1.0
Phi = expm(A * t)
print(f"Matrix Exponential e^(At) at t={t}:\n{Phi}")

# 2. Generate Phase Portrait using streamplot
w = 2.0
x, y = np.mgrid[-w:w:20j, -w:w:20j]
u = A[0,0]*x + A[0,1]*y
v = A[1,0]*x + A[1,1]*y

plt.figure(figsize=(8, 8))
plt.streamplot(x, y, u, v, color='cornflowerblue')
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.title("Phase Portrait of a Stable Spiral")
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)
plt.show()