Ordinary Differential Equations (6): Linear Systems of Differential Equations

When multiple variables influence each other, a single differential equation is no longer sufficient. The populations of predators and prey, currents and voltages in circuits, concentrations of substances in chemical reactions — the relationships between them need to be described by systems of differential equations. In this chapter, we'll see that the combination of linear algebra and differential equations produces an extremely elegant theory: matrix exponentials, fundamental matrices, and phase space analysis — these tools make the behavior of complex systems clearly visible.

Starting with an Ecology Problem

Imagine a simplified ecosystem: rabbits and wolves. Rabbits eat grass and reproduce, wolves eat rabbits and reproduce. Letbe the number of rabbits andbe the number of wolves.

Intuitive modeling: - The growth rate of rabbits depends on their own number (more means faster reproduction) minus the number eaten by wolves - The growth rate of wolves depends on how many rabbits are available to eat minus natural death

A simplified linear model: This is a linear system of differential equations. How do we solve it? How do we understand the long-term behavior ofand?

Vector Form and Matrix Notation

Vector Representation

Write the system in vector form:Abbreviated as:whereis the state vector andis the coefficient matrix.

General Form

Homogeneous linear system:

Non-homogeneous linear system:This chapter mainly focuses on the constant coefficient case:does not depend on.

Review: The One-Dimensional Case

For the one-dimensional equation, the solution is.

Natural conjecture: For, is the solution?

Yes! But first we need to define the matrix exponential.

Matrix Exponential

Definition

Recall the Taylor expansion of the scalar exponential function:

Definition of matrix exponential:This series converges for any square matrix

Basic Properties

Property 1: (exponential of zero matrix is identity matrix)

Property 2:This is exactly the differentiation property we need!

Property 3: If (andcommute), then

Warning: If, then generallyThis is an important difference between matrices and scalars.

Property 4:is always invertible,

Property 5:

Solution of Linear Systems

Theorem: The initial value problemhas unique solution:

Verification: - ✓ - ✓

How to Compute the Matrix Exponential?

Direct computation using the definition is tedious. There are several more practical methods.

Method 1: Diagonalization

Ifis diagonalizable:, where.

Then:$

$ And the exponential of a diagonal matrix is simple:

Method 2: Eigenvalues and Eigenvectors

Lethave eigenvaluesand corresponding eigenvectors.

If the eigenvectors are linearly independent, the general solution is:

Why? Because, sois a solution of:

Example: 2×2 Matrix

Consider the ecology model from the beginning:

Step 1: Find eigenvalues

Characteristic polynomial:$

Complex eigenvalues! Step 2: Find eigenvectors

For:Take Step 3: Write the general solution

For complex eigenvaluesand complex eigenvector:

Real part solution:Imaginary part solution:This gives spiral trajectories: oscillating and growing (because).

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

# Define the system
A = np.array([[2, -1], 
              [1, 0.5]])

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

# Define differential equation
def system(x, t, A):
    return A @ x

# Solve for multiple initial values
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

t = np.linspace(0, 5, 500)

# Left plot: phase space trajectories
ax1 = axes[0]
initial_conditions = [
    [1, 0], [0, 1], [-1, 0], [0, -1],
    [1, 1], [-1, 1], [1, -1], [-1, -1],
    [2, 0], [0, 2]
]

for x0 in initial_conditions:
    sol = odeint(system, x0, t, args=(A,))
    ax1.plot(sol[:, 0], sol[:, 1], linewidth=1.5, alpha=0.7)
    ax1.scatter(x0[0], x0[1], s=50, zorder=3)

# Draw eigenvector directions
for i, (ev, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
    if np.isreal(ev):
        scale = 2
        ax1.arrow(0, 0, scale*vec[0].real, scale*vec[1].real, 
                 head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
                 linewidth=2, label=f'$\\lambda={ev:.2f}$')

ax1.set_xlabel('x (rabbits)', fontsize=12)
ax1.set_ylabel('y (wolves)', fontsize=12)
ax1.set_title('Phase Portrait: Unstable Spiral', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(-3, 3)
ax1.set_ylim(-3, 3)
ax1.set_aspect('equal')

# Right plot: time evolution
ax2 = axes[1]
x0 = [1, 0.5]
sol = odeint(system, x0, t, args=(A,))

ax2.plot(t, sol[:, 0], 'b-', linewidth=2, label='x(t) - Rabbits')
ax2.plot(t, sol[:, 1], 'r-', linewidth=2, label='y(t) - Wolves')
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Population', fontsize=12)
ax2.set_title('Time Evolution: Oscillating Growth', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('linear_system_spiral.png', dpi=150, bbox_inches='tight')
plt.show()

# Verify matrix exponential
t_test = 1.0
x0 = np.array([1, 0.5])

# Method 1: using scipy
exp_At = expm(A * t_test)
x_expm = exp_At @ x0

# Method 2: using ODE solver
sol_ode = odeint(system, x0, [0, t_test], args=(A,))[-1]

print(f"\nSolution at t={t_test}:")
print(f"  Matrix exponential method: {x_expm}")
print(f"  ODE solver:   {sol_ode}")
print(f"  Error: {np.linalg.norm(x_expm - sol_ode):.2e}")

Phase Plane Analysis: Classification of 2×2 Systems

Eigenvalues Determine Everything

For 2×2 real matrices, based on the properties of eigenvalues, phase portraits fall into several types:

Case 1: Two distinct real eigenvalues

Letboth be real.

• If: stable node (all trajectories approach origin)
• If: unstable node (all trajectories leave origin)
• If: saddle point (unstable, has stable/unstable manifolds)

Case 2: Complex eigenvalues - If: stable spiral (spirals toward origin) - If: unstable spiral (spirals away from origin) - If: center (closed orbits, periodic motion)

Case 3: Repeated eigenvalues

• If there are two linearly independent eigenvectors: star node
• If there is only one eigenvector: degenerate node (need generalized eigenvectors)

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

def plot_phase_portrait(A, ax, title, xlim=(-3, 3), ylim=(-3, 3)):
    """Draw phase portrait"""
    def system(x, t):
        return A @ x
    
    # Grid points
    x_range = np.linspace(xlim[0], xlim[1], 20)
    y_range = np.linspace(ylim[0], ylim[1], 20)
    X, Y = np.meshgrid(x_range, y_range)
    
    # Compute direction field
    U = A[0, 0] * X + A[0, 1] * Y
    V = A[1, 0] * X + A[1, 1] * Y
    
    # Normalize
    M = np.sqrt(U**2 + V**2)
    M[M == 0] = 1
    U_norm = U / M
    V_norm = V / M
    
    # Plot direction field
    ax.quiver(X, Y, U_norm, V_norm, M, cmap='coolwarm', alpha=0.5)
    
    # Plot trajectories
    t = np.linspace(0, 10, 500)
    t_back = np.linspace(0, -10, 500)
    
    # Select initial values
    angles = np.linspace(0, 2*np.pi, 8, endpoint=False)
    r = 2.5
    for angle in angles:
        x0 = [r * np.cos(angle), r * np.sin(angle)]
        
        # Forward integration
        sol = odeint(system, x0, t)
        ax.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=1, alpha=0.7)
        
        # Backward integration
        sol_back = odeint(system, x0, t_back)
        ax.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=1, alpha=0.7)
    
    # Draw eigenvectors
    eigenvalues, eigenvectors = np.linalg.eig(A)
    for i, (ev, vec) in enumerate(zip(eigenvalues, eigenvectors.T)):
        if np.isreal(ev) and np.isreal(vec).all():
            vec = vec.real
            scale = 2.5
            ax.arrow(0, 0, scale*vec[0], scale*vec[1], 
                    head_width=0.15, head_length=0.1, fc='red', ec='red',
                    linewidth=2)
            ax.arrow(0, 0, -scale*vec[0], -scale*vec[1], 
                    head_width=0.15, head_length=0.1, fc='red', ec='red',
                    linewidth=2)
    
    ax.plot(0, 0, 'ko', markersize=8)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    ax.set_xlabel('x', fontsize=11)
    ax.set_ylabel('y', fontsize=11)
    ax.set_title(title, fontsize=12, fontweight='bold')
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    # Display eigenvalues
    ev_str = ', '.join([f'{ev:.2f}' for ev in eigenvalues])
    ax.text(0.05, 0.95, f'λ = {ev_str}', transform=ax.transAxes, 
            fontsize=9, verticalalignment='top',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# Create 6 typical cases
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 1. Stable node (λ1, λ2 < 0)
A1 = np.array([[-2, 0], [0, -1]])
plot_phase_portrait(A1, axes[0, 0], 'Stable Node')

# 2. Unstable node (λ1, λ2 > 0)
A2 = np.array([[2, 0], [0, 1]])
plot_phase_portrait(A2, axes[0, 1], 'Unstable Node')

# 3. Saddle point (λ1 < 0 < λ2)
A3 = np.array([[-1, 0], [0, 2]])
plot_phase_portrait(A3, axes[0, 2], 'Saddle Point')

# 4. Stable spiral (α < 0)
A4 = np.array([[-0.5, 1], [-1, -0.5]])
plot_phase_portrait(A4, axes[1, 0], 'Stable Spiral')

# 5. Unstable spiral (α > 0)
A5 = np.array([[0.5, 1], [-1, 0.5]])
plot_phase_portrait(A5, axes[1, 1], 'Unstable Spiral')

# 6. Center (α = 0)
A6 = np.array([[0, 1], [-1, 0]])
plot_phase_portrait(A6, axes[1, 2], 'Center')

plt.tight_layout()
plt.savefig('phase_portraits.png', dpi=150, bbox_inches='tight')
plt.show()

Fundamental Matrix

Definition

A fundamental matrixis anmatrix whose columns arelinearly independent solutions:

Properties: 1.$'(t) = A(t)((t)) t(t) = (t)$is a constant vector

Relationship with Matrix Exponential

Theorem: If, then.

More generally,.

Liouville's Formula

Theorem (Liouville/Abel):For the constant coefficient case:

Physical significance: The rate of evolution of volume elements in phase space equals.

If, volume contracts (dissipative system); If, volume is conserved (Hamiltonian system); If, volume expands.

Repeated Eigenvalues and Generalized Eigenvectors

The Problem

When eigenvalues have repeated roots, there may not be enough eigenvectors. For example:Eigenvalue(repeated), but there is only one eigenvector.

Generalized Eigenvectors

Definition: Ifbut, thenis a generalized eigenvector of rank for.

For the example above:Take, then.

Corresponding Solutions

For a generalized eigenvector of rank, the corresponding solution contains polynomial factors:For the example above, the two independent solutions are:

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

# Example with repeated eigenvalue
A = np.array([[3, 1], [0, 3]])

print("Matrix A:")
print(A)
print("\nEigenvalues:", np.linalg.eigvals(A))

# Phase portrait
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

def system(x, t):
    return A @ x

# Trajectories
ax1 = axes[0]
t = np.linspace(0, 1, 200)
t_back = np.linspace(0, -1, 200)

initial_conditions = [
    [1, 0.5], [-1, 0.5], [1, -0.5], [-1, -0.5],
    [0.1, 1], [0.1, -1], [-0.1, 1], [-0.1, -1]
]

for x0 in initial_conditions:
    sol = odeint(system, x0, t)
    sol_back = odeint(system, x0, t_back)
    ax1.plot(sol[:, 0], sol[:, 1], 'b-', linewidth=1.5, alpha=0.7)
    ax1.plot(sol_back[:, 0], sol_back[:, 1], 'b-', linewidth=1.5, alpha=0.7)
    ax1.scatter(x0[0], x0[1], s=50, zorder=3)

# Eigenvector direction
ax1.arrow(0, 0, 1, 0, head_width=0.1, head_length=0.05, fc='red', ec='red', linewidth=2)
ax1.arrow(0, 0, -1, 0, head_width=0.1, head_length=0.05, fc='red', ec='red', linewidth=2)

ax1.plot(0, 0, 'ko', markersize=8)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Degenerate Node (Repeated Eigenvalue λ=3)', fontsize=14, fontweight='bold')
ax1.set_xlim(-2, 2)
ax1.set_ylim(-2, 2)
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)

# Time evolution
ax2 = axes[1]
t_long = np.linspace(-1, 1, 500)
x0 = [0.1, 1]
sol_long = odeint(system, x0, t_long)

ax2.plot(t_long, sol_long[:, 0], 'b-', linewidth=2, label='x(t)')
ax2.plot(t_long, sol_long[:, 1], 'r-', linewidth=2, label='y(t)')

# Analytical solution comparison
t_pos = t_long[t_long >= 0]
# x(t) = (c1 + c2*t) * e^(3t), y(t) = c2 * e^(3t)
# Initial values x(0)=0.1, y(0)=1 => c1=0.1, c2=1
c1, c2 = 0.1, 1
x_analytical = (c1 + c2*t_long) * np.exp(3*t_long)
y_analytical = c2 * np.exp(3*t_long)

ax2.plot(t_long, x_analytical, 'b--', linewidth=1.5, alpha=0.7, label='x(t) analytical')
ax2.plot(t_long, y_analytical, 'r--', linewidth=1.5, alpha=0.7, label='y(t) analytical')

ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Value', fontsize=12)
ax2.set_title('Time Evolution with Polynomial Growth', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-5, 10)

plt.tight_layout()
plt.savefig('repeated_eigenvalue.png', dpi=150, bbox_inches='tight')
plt.show()

Non-homogeneous Systems: Variation of Parameters

The Problem

Solve the non-homogeneous system:

Variation of Parameters

Idea: Assume the solution has the form, whereis the fundamental matrix.

Substituting into the equation:Since:$

$$

Final solution:For the constant coefficient case ():This is the famous Duhamel's formula.

Example: Forced Harmonic Oscillator

Consider a spring system driven by periodic force:Writing as a first-order system (let):

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

# Forced harmonic oscillator
omega0 = 1.0  # Natural frequency
F0 = 0.5  # Driving force amplitude
gamma = 0.1  # Damping coefficient

def forced_oscillator(X, t, omega_drive):
    x, v = X
    dxdt = v
    dvdt = -omega0**2 * x - 2*gamma*v + F0 * np.cos(omega_drive * t)
    return [dxdt, dvdt]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 100, 5000)
X0 = [1, 0]

# Different driving frequencies
driving_freqs = [0.5, 1.0, 1.5, 2.0]

# Steady-state amplitude vs driving frequency (theory)
omega_range = np.linspace(0.1, 2.5, 200)
amplitude_theory = F0 / np.sqrt((omega0**2 - omega_range**2)**2 + (2*gamma*omega_range)**2)

ax1 = axes[0, 0]
ax1.plot(omega_range, amplitude_theory, 'b-', linewidth=2)
ax1.axvline(x=omega0, color='r', linestyle='--', label=f'ω₀ = {omega0}')
ax1.set_xlabel('Driving Frequency ω', fontsize=12)
ax1.set_ylabel('Steady-State Amplitude', fontsize=12)
ax1.set_title('Resonance Curve', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Time evolution: resonance case
ax2 = axes[0, 1]
omega_resonance = omega0
sol_res = odeint(forced_oscillator, X0, t, args=(omega_resonance,))
ax2.plot(t, sol_res[:, 0], 'b-', linewidth=1)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Displacement x', fontsize=12)
ax2.set_title(f'At Resonance: ω = ω₀ = {omega0}', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Time evolution: off resonance
ax3 = axes[1, 0]
omega_off = 2.0
sol_off = odeint(forced_oscillator, X0, t, args=(omega_off,))
ax3.plot(t, sol_off[:, 0], 'g-', linewidth=1)
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Displacement x', fontsize=12)
ax3.set_title(f'Off Resonance: ω = {omega_off}', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)

# Phase space
ax4 = axes[1, 1]
# Phase space trajectories at steady state
t_steady = t[t > 50]
sol_res_steady = sol_res[t > 50]
sol_off_steady = sol_off[t > 50]

ax4.plot(sol_res_steady[:, 0], sol_res_steady[:, 1], 'b-', linewidth=1, 
         label=f'ω = {omega_resonance} (resonance)', alpha=0.7)
ax4.plot(sol_off_steady[:, 0], sol_off_steady[:, 1], 'g-', linewidth=1, 
         label=f'ω = {omega_off} (off resonance)', alpha=0.7)
ax4.set_xlabel('x', fontsize=12)
ax4.set_ylabel('v', fontsize=12)
ax4.set_title('Phase Space (Steady State)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_aspect('equal')

plt.tight_layout()
plt.savefig('forced_oscillator.png', dpi=150, bbox_inches='tight')
plt.show()

Application: Coupled Oscillators

Problem Setup

Two masses connected by springs:

Wall ----/\/\/\---- [m1] ----/\/\/\---- [m2] ----/\/\/\---- Wall
       k1              k2              k3

Letbe the displacements of the two masses from their equilibrium positions.

Equations of motion:$

$ Simplifying (let,):$

$

Matrix Form

Let:

Normal Modes

Eigenvalue analysis reveals the system's normal modes— patterns where the entire system oscillates synchronously at specific frequencies.

For a symmetric coupled system, normal modes have simple physical interpretations: - Symmetric mode: Both masses move in the same direction () - Antisymmetric mode: Masses move in opposite directions ()

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.linalg import eig

# Coupled oscillator system
# Set m1 = m2 = 1, k1 = k2 = k3 = 1

def coupled_oscillators(X, t):
    x1, v1, x2, v2 = X
    dx1 = v1
    dv1 = -2*x1 + x2
    dx2 = v2
    dv2 = x1 - 2*x2
    return [dx1, dv1, dx2, dv2]

# System matrix
A = np.array([
    [0, 1, 0, 0],
    [-2, 0, 1, 0],
    [0, 0, 0, 1],
    [1, 0, -2, 0]
])

# Eigenvalue analysis
eigenvalues, eigenvectors = eig(A)
print("Eigenvalues:", eigenvalues)
print("\nNormal frequencies (absolute value of imaginary parts):")
freqs = np.abs(eigenvalues.imag)
unique_freqs = np.unique(np.round(freqs, 6))
print(unique_freqs[unique_freqs > 0])

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 30, 1000)

# Normal mode 1: symmetric mode (x1(0) = x2(0) = 1)
ax1 = axes[0, 0]
X0_sym = [1, 0, 1, 0]
sol_sym = odeint(coupled_oscillators, X0_sym, t)
ax1.plot(t, sol_sym[:, 0], 'b-', linewidth=2, label='$x_1(t)$')
ax1.plot(t, sol_sym[:, 2], 'r--', linewidth=2, label='$x_2(t)$')
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Displacement', fontsize=12)
ax1.set_title('Symmetric Mode:$x_1(0) = x_2(0) = 1$', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)

# Normal mode 2: antisymmetric mode (x1(0) = 1, x2(0) = -1)
ax2 = axes[0, 1]
X0_antisym = [1, 0, -1, 0]
sol_antisym = odeint(coupled_oscillators, X0_antisym, t)
ax2.plot(t, sol_antisym[:, 0], 'b-', linewidth=2, label='$x_1(t)$')
ax2.plot(t, sol_antisym[:, 2], 'r--', linewidth=2, label='$x_2(t)$')
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Displacement', fontsize=12)
ax2.set_title('Antisymmetric Mode:$x_1(0) = 1, x_2(0) = -1$', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

# General initial condition: beat phenomenon
ax3 = axes[1, 0]
X0_general = [1, 0, 0, 0]  # Only excite the first mass
sol_general = odeint(coupled_oscillators, X0_general, t)
ax3.plot(t, sol_general[:, 0], 'b-', linewidth=1.5, label='$x_1(t)$')
ax3.plot(t, sol_general[:, 2], 'r-', linewidth=1.5, label='$x_2(t)$')
ax3.set_xlabel('Time', fontsize=12)
ax3.set_ylabel('Displacement', fontsize=12)
ax3.set_title('General Initial Condition: Beat Phenomenon', fontsize=14, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(True, alpha=0.3)

# Energy exchange
ax4 = axes[1, 1]
# Kinetic + potential energy
KE1 = 0.5 * sol_general[:, 1]**2
KE2 = 0.5 * sol_general[:, 3]**2
PE1 = 0.5 * sol_general[:, 0]**2  # Only considering boundary springs
PE2 = 0.5 * sol_general[:, 2]**2
PE_coupling = 0.5 * (sol_general[:, 2] - sol_general[:, 0])**2

E1 = KE1 + PE1 + 0.5*PE_coupling
E2 = KE2 + PE2 + 0.5*PE_coupling
E_total = KE1 + KE2 + PE1 + PE2 + PE_coupling

ax4.plot(t, E1, 'b-', linewidth=1.5, label='Energy in mass 1')
ax4.plot(t, E2, 'r-', linewidth=1.5, label='Energy in mass 2')
ax4.plot(t, E_total, 'k--', linewidth=1.5, label='Total energy')
ax4.set_xlabel('Time', fontsize=12)
ax4.set_ylabel('Energy', fontsize=12)
ax4.set_title('Energy Transfer Between Oscillators', fontsize=14, fontweight='bold')
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('coupled_oscillators.png', dpi=150, bbox_inches='tight')
plt.show()

# Print normal frequencies
omega1 = 1  # Symmetric mode frequency sqrt(k1/m) = sqrt(1)
omega2 = np.sqrt(3)  # Antisymmetric mode frequency sqrt((k1+2k2)/m) = sqrt(3)
print(f"\nNormal mode frequencies:")
print(f"  Symmetric mode: ω₁ = {omega1:.4f}")
print(f"  Antisymmetric mode: ω₂ = {omega2:.4f}")
print(f"  Beat frequency: Δω = {omega2 - omega1:.4f}")

Application: Circuit Networks

RLC Circuit

Consider a simple series RLC circuit:where.

Writing as a first-order system (state variables: chargeand current):

Eigenvalues and Circuit Behavior

Let(resonant frequency),(damping coefficient).

Eigenvalues: - Underdamped (): Complex eigenvalues, oscillatory decay - Critically damped (): Repeated real eigenvalue, fastest approach to steady state - Overdamped (): Two real eigenvalues, monotonic decay

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

# RLC circuit response
L = 1.0  # Inductance
C = 1.0  # Capacitance
omega0 = 1 / np.sqrt(L * C)

def rlc_circuit(X, t, R, V_func):
    Q, I = X
    dQdt = I
    dIdt = (V_func(t) - R*I - Q/C) / L
    return [dQdt, dIdt]

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

t = np.linspace(0, 20, 1000)
X0 = [1, 0]  # Initial charge 1, current 0
V_func = lambda t: 0  # No external voltage, studying free response

# Different damping cases
R_values = [0.5, 2.0, 4.0]  # Underdamped, critically damped, overdamped
labels = ['Underdamped', 'Critically damped', 'Overdamped']
colors = ['b', 'g', 'r']

# Charge response
ax1 = axes[0, 0]
for R, label, color in zip(R_values, labels, colors):
    gamma = R / (2*L)
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax1.plot(t, sol[:, 0], color=color, linewidth=2, 
             label=f'{label} (R={R}, γ/ω₀={gamma/omega0:.2f})')

ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('Time', fontsize=12)
ax1.set_ylabel('Charge Q', fontsize=12)
ax1.set_title('RLC Circuit: Charge Response', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Current response
ax2 = axes[0, 1]
for R, label, color in zip(R_values, labels, colors):
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax2.plot(t, sol[:, 1], color=color, linewidth=2, label=f'{label}')

ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('Time', fontsize=12)
ax2.set_ylabel('Current I', fontsize=12)
ax2.set_title('RLC Circuit: Current Response', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# Phase space
ax3 = axes[1, 0]
for R, label, color in zip(R_values, labels, colors):
    sol = odeint(rlc_circuit, X0, t, args=(R, V_func))
    ax3.plot(sol[:, 0], sol[:, 1], color=color, linewidth=2, label=f'{label}')
    ax3.scatter(X0[0], X0[1], s=100, color=color, zorder=3)

ax3.plot(0, 0, 'ko', markersize=8)
ax3.set_xlabel('Charge Q', fontsize=12)
ax3.set_ylabel('Current I', fontsize=12)
ax3.set_title('Phase Space', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_aspect('equal')

# Frequency response
ax4 = axes[1, 1]
omega_range = np.linspace(0.1, 3, 200)
R = 0.5  # Underdamped

# Steady-state amplitude (theory)
# |H(j ω)| = 1 / sqrt((1 - (ω/ω0)²)² + (2γω/ω0²)²)
gamma = R / (2*L)
amplitude = 1 / np.sqrt((1 - (omega_range/omega0)**2)**2 + (2*gamma*omega_range/omega0**2)**2)

ax4.plot(omega_range, amplitude, 'b-', linewidth=2)
ax4.axvline(x=omega0, color='r', linestyle='--', label=f'ω₀ = {omega0}')
ax4.set_xlabel('Driving Frequency ω', fontsize=12)
ax4.set_ylabel('Amplitude Response', fontsize=12)
ax4.set_title('Frequency Response (R=0.5)', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('rlc_circuit.png', dpi=150, bbox_inches='tight')
plt.show()

Higher-Dimensional Systems and Jordan Normal Form

Jordan Normal Form

When a matrix is not diagonalizable, it can be reduced to Jordan normal form:whereis a Jordan matrix composed of Jordan blocks.

AJordan block:

Exponential of a Jordan Block

Key observation: Even if all eigenvalues have negative real parts, polynomial factors may cause short-term growth, but long-term decay still occurs.

Introduction to Stability

Stability of the Zero Solution

Consider the zero solutionof.

Theorem: 1. If all eigenvalues ofhave real parts, then the zero solution is asymptotically stable 2. Ifhas any eigenvalue with real part, then the zero solution is unstable 3. If all eigenvalues have real parts, and eigenvalues with real partare all simple (algebraic multiplicity equals geometric multiplicity), then the zero solution is stable but not asymptotically stable

This sets the stage for detailed stability theory in the next chapter.

Numerical Methods

Numerical Computation of Matrix Exponential

There are several methods for computing:

Method 1: Pad é approximation (scipy.linalg.expm)

The most stable general-purpose method.

Method 2: Diagonalization

Ifis diagonalizable, this is the fastest method.

Method 3: Series truncation

Directly using the definition, suitable for small matrices and small.

import numpy as np
from scipy.linalg import expm, eig
import time

def matrix_exp_series(A, n_terms=50):
    """Compute matrix exponential using series"""
    result = np.eye(len(A))
    A_power = np.eye(len(A))
    factorial = 1.0
    
    for k in range(1, n_terms):
        A_power = A_power @ A
        factorial *= k
        result = result + A_power / factorial
    
    return result

def matrix_exp_diag(A):
    """Compute matrix exponential using diagonalization"""
    eigenvalues, P = eig(A)
    exp_D = np.diag(np.exp(eigenvalues))
    return P @ exp_D @ np.linalg.inv(P)

# Test
A = np.array([[1, 2], [3, 4]], dtype=float)

print("Matrix A:")
print(A)

# Three methods
exp_A_scipy = expm(A)
exp_A_series = matrix_exp_series(A, 100)
exp_A_diag = matrix_exp_diag(A)

print("\nexp(A) - scipy.linalg.expm:")
print(exp_A_scipy)

print("\nexp(A) - series expansion:")
print(exp_A_series)

print("\nexp(A) - diagonalization:")
print(np.real(exp_A_diag))

print("\nError (series vs scipy):", np.linalg.norm(exp_A_series - exp_A_scipy))
print("Error (diagonalization vs scipy):", np.linalg.norm(np.real(exp_A_diag) - exp_A_scipy))

# Performance comparison
n = 100
A_large = np.random.randn(n, n)

start = time.time()
for _ in range(10):
    expm(A_large)
print(f"\n{n}x{n} matrix, scipy.expm: {(time.time()-start)/10*1000:.2f} ms")

start = time.time()
for _ in range(10):
    matrix_exp_series(A_large, 30)
print(f"{n}x{n} matrix, series expansion: {(time.time()-start)/10*1000:.2f} ms")

ODE Solvers

For large systems or long-time integration, directly using matrix exponential may not be the best choice.

scipy.integrate.odeint and solve_ivp provide various adaptive step-size methods:

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

# Comparing different ODE solvers
A = np.array([[-1, 100], [0, -1]])  # Stiff system

def system(t, x):
    return A @ x

x0 = [1, 1]
t_span = (0, 10)
t_eval = np.linspace(0, 10, 1000)

# Different methods
methods = ['RK45', 'RK23', 'DOP853', 'Radau', 'BDF', 'LSODA']

fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()

for ax, method in zip(axes, methods):
    try:
        sol = solve_ivp(system, t_span, x0, method=method, t_eval=t_eval)
        ax.plot(sol.t, sol.y[0], 'b-', linewidth=1.5, label='x ₁(t)')
        ax.plot(sol.t, sol.y[1], 'r-', linewidth=1.5, label='x ₂(t)')
        ax.set_title(f'{method} (nfev={sol.nfev})', fontsize=12, fontweight='bold')
    except Exception as e:
        ax.text(0.5, 0.5, f'Error: {e}', transform=ax.transAxes, ha='center')
        ax.set_title(method, fontsize=12)
    
    ax.set_xlabel('Time', fontsize=10)
    ax.set_ylabel('Value', fontsize=10)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('ode_solvers.png', dpi=150, bbox_inches='tight')
plt.show()

Summary

Key Points of This Chapter

1. Linear systems:, solution is

2. Matrix exponential:

   • Definition: - Computation: diagonalization, Jordan normal form, Pad é approximation

3. Phase space analysis:

   • Eigenvalues determine equilibrium point types
   • Nodes, spirals, saddle points, centers

4. Fundamental matrix: Matrix whose columns are independent solutions

5. Non-homogeneous systems: Variation of parameters, Duhamel's formula

6. Applications: Coupled oscillators, circuits, ecology models

Preview of Next Chapter

In "Stability Theory," we will: - Study Lyapunov stability in depth - Analyze local behavior of nonlinear systems - Learn the Lyapunov function method - Explore the beginnings of bifurcation and chaos

Exercises

Basic Problems

1. Compute the exponential of matrix, i.e.,.

2. Solve the initial value problem:3. Determine the type of the origin (node/spiral/saddle/center) for the following systems: - - -$A =

   (e^A) = e^{(A)}$.

3. For(repeated eigenvalue), find the fundamental matrix.

Advanced Problems

6. Coupled harmonic oscillators:

   • Find the normal mode frequencies of the system
   • Prove conservation of energy
   • Analyze the beat frequency period

7. RLC circuit:

   • Prove the critical damping condition - Calculate the physical meaning of quality factor - Design a filter with bandwidth 10 Hz and center frequency 1000 Hz

8. Prove that if all eigenvalues ofhave negative real parts, then.

9. Chemical reaction: Consider the reaction:

• Write in matrix form
• Prove conservation of total mass
• Find steady-state concentrations

Programming Problems

10. Implement function matrix_exp(A, t) to computeusing three methods and compare accuracy.

11. Write a program to plot a complete phase portrait for 2×2 systems, including:

    • Direction field
    • Eigenvectors
    • Multiple trajectories
    • Automatic detection and labeling of equilibrium point type

12. Simulate a simplified three-body problem: three masses connected by springs oscillating on a line.

    • Find the three normal modes
    • Visualize energy flow among the three masses

13. Implement a circuit simulator that can analyze the response of arbitrary RLC networks.

Thought Questions

14. Why is(when)? Give a specific counterexample.

15. Why can't trajectories in phase space cross (except at equilibrium points)? What is the relationship with uniqueness of solutions?

16. How can we intuitively understand that "zero trace means phase space volume is conserved"?

References

1. Hirsch, M. W., Smale, S., & Devaney, R. L. (2012). Differential Equations, Dynamical Systems, and an Introduction to Chaos. Academic Press.
2. Strang, G. (2019). Linear Algebra and Learning from Data. Wellesley-Cambridge Press.
3. Perko, L. (2001). Differential Equations and Dynamical Systems. Springer.
4. MIT 18.03SC: Differential Equations, Unit III: Linear Systems.
5. Moler, C., & Van Loan, C. (2003). Nineteen dubious ways to compute the exponential of a matrix. SIAM Review.

---

Next Chapter: Stability Theory →

---

This is Chapter 6 of the "World of Ordinary Differential Equations" series.