Ordinary Differential Equations (18): Advanced Topics and Summary
Chen Kai BOSS
  2019-06-30 11:15:00 2019-06-30 11:15      3.5k Words  22 Mins  

Our differential equations journey is nearing its end, but mathematics itself is endless. This chapter introduces several active research frontiers — Neural ODEs, delay differential equations, stochastic differential equations, and fractional differential equations — which are profoundly changing our understanding of dynamical systems. Finally, we review the core content of the entire series and provide readers with a complete learning roadmap.

Neural ODEs

From ResNet to Continuous Depth

Residual networks (ResNet) in deep learning can be written as:

As the number of layers approaches infinity and changes per layer become infinitesimal, this becomes an ODE:This is the core idea of Neural ODEs.

Advantages of Neural ODEs

  1. Memory efficiency: No need to store all intermediate states
  2. Adaptive computation: Complex inputs can use more "layers"
  3. Continuous-time modeling: Naturally handles irregular time series
  4. Reversibility: Both forward and backward computation are exact

Backpropagation and the Adjoint Method

Traditional backpropagation requires storing all intermediate activations. Neural ODEs use the adjoint method:

Define the adjoint state, which satisfies:Parameter gradients: 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

def neural_ode_demo():
    """Neural ODE demonstration: learning spiral dynamics"""
    
    # Generate spiral data
    def true_dynamics(t, y):
        """True spiral dynamics"""
        return np.array([
            -0.1 * y[0] + 2.0 * y[1],
            -2.0 * y[0] - 0.1 * y[1]
        ])
    
    # Generate training data
    np.random.seed(42)
    t_span = np.linspace(0, 10, 100)
    y0 = np.array([2.0, 0.0])
    
    # True trajectory
    sol = solve_ivp(true_dynamics, [0, 10], y0, t_eval=t_span, method='RK45')
    true_y = sol.y.T
    
    # Add noise as observations
    noise = np.random.normal(0, 0.1, true_y.shape)
    observed_y = true_y + noise
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # True trajectory
    ax1 = axes[0, 0]
    ax1.plot(true_y[:, 0], true_y[:, 1], 'b-', linewidth=2, label='True trajectory')
    ax1.scatter(observed_y[::5, 0], observed_y[::5, 1], c='r', s=30, alpha=0.5, label='Observations')
    ax1.plot(y0[0], y0[1], 'go', markersize=10, label='Start')
    ax1.set_xlabel('$y_1$', fontsize=12)
    ax1.set_ylabel('$y_2$', fontsize=12)
    ax1.set_title('Spiral Dynamics', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_aspect('equal')
    
    # Time series
    ax2 = axes[0, 1]
    ax2.plot(t_span, true_y[:, 0], 'b-', linewidth=2, label='$y_1$(true)')
    ax2.plot(t_span, true_y[:, 1], 'r-', linewidth=2, label='$y_2$(true)')
    ax2.scatter(t_span[::5], observed_y[::5, 0], c='b', s=20, alpha=0.5)
    ax2.scatter(t_span[::5], observed_y[::5, 1], c='r', s=20, alpha=0.5)
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Value', fontsize=12)
    ax2.set_title('Time Series', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # ResNet vs Neural ODE conceptual diagram
    ax3 = axes[1, 0]
    
    # ResNet (discrete)
    n_layers = 10
    y_discrete = [y0.copy()]
    y_current = y0.copy()
    dt = 1.0
    for i in range(n_layers):
        dy = true_dynamics(i*dt, y_current) * dt
        y_current = y_current + dy
        y_discrete.append(y_current.copy())
    y_discrete = np.array(y_discrete)
    
    ax3.plot(y_discrete[:, 0], y_discrete[:, 1], 'ro-', linewidth=2, markersize=8, 
             label=f'ResNet ({n_layers} layers)')
    ax3.plot(true_y[:, 0], true_y[:, 1], 'b-', linewidth=2, alpha=0.5, label='Neural ODE (continuous)')
    ax3.set_xlabel('$y_1$', fontsize=12)
    ax3.set_ylabel('$y_2$', fontsize=12)
    ax3.set_title('ResNet vs Neural ODE', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_aspect('equal')
    
    # Vector field
    ax4 = axes[1, 1]
    
    y1_range = np.linspace(-2.5, 2.5, 20)
    y2_range = np.linspace(-2.5, 2.5, 20)
    Y1, Y2 = np.meshgrid(y1_range, y2_range)
    
    DY1 = -0.1 * Y1 + 2.0 * Y2
    DY2 = -2.0 * Y1 - 0.1 * Y2
    
    M = np.sqrt(DY1**2 + DY2**2)
    M[M == 0] = 1
    DY1_norm = DY1 / M
    DY2_norm = DY2 / M
    
    ax4.quiver(Y1, Y2, DY1_norm, DY2_norm, M, cmap='coolwarm', alpha=0.6)
    ax4.plot(true_y[:, 0], true_y[:, 1], 'b-', linewidth=2)
    ax4.set_xlabel('$y_1$', fontsize=12)
    ax4.set_ylabel('$y_2$', fontsize=12)
    ax4.set_title('Learned Vector Field', fontsize=14, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    ax4.set_aspect('equal')
    
    plt.tight_layout()
    plt.savefig('neural_ode_demo.png', dpi=150, bbox_inches='tight')
    plt.show()

neural_ode_demo()

Continuous Normalizing Flows

An important application of Neural ODEs is Continuous Normalizing Flows.

Given an initial distribution, through ODE transformation:The probability density change satisfies:This allows transforming simple distributions (like Gaussian) into complex distributions for generative models.

Delay Differential Equations (DDEs)

Basic Concepts

Many systems' current rate of change depends not only on the current state but also on past states:whereis the delay time.

Examples: - Population dynamics: Adult reproduction depends on individuals bornyears ago - Control systems: Sensor and actuator response delays - Neural networks: Signal transmission delays

Delay Logistic EquationThis is Hutchinson's equation, describing delayed density-dependent growth of a single species.

Surprising result: When, the equilibrium becomes unstable and periodic oscillations emerge!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def delay_differential_equations():
    """Delay differential equation examples"""
    
    from collections import deque
    
    def solve_dde_simple(f, history, t_span, tau, dt=0.01):
        """Simple DDE solver (Euler method)"""
        t_start, t_end = t_span
        n_steps = int((t_end - t_start) / dt)
        n_delay = int(tau / dt)
        
        # Initialize history
        t_history = np.arange(t_start - tau, t_start + dt, dt)
        if callable(history):
            x_history = [history(t) for t in t_history]
        else:
            x_history = [history] * len(t_history)
        
        x_buffer = deque(x_history, maxlen=n_delay + 1)
        
        t_values = [t_start]
        x_values = [x_buffer[-1]]
        
        t = t_start
        for _ in range(n_steps):
            x_current = x_buffer[-1]
            x_delayed = x_buffer[0]
            
            dx = f(t, x_current, x_delayed)
            x_new = x_current + dt * dx
            
            x_buffer.append(x_new)
            t += dt
            
            t_values.append(t)
            x_values.append(x_new)
        
        return np.array(t_values), np.array(x_values)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Delay Logistic equation
    r = 1.0
    K = 1.0
    
    def delay_logistic(t, N, N_delayed, r, K):
        return r * N * (1 - N_delayed / K)
    
    ax1 = axes[0, 0]
    
    tau_values = [0.5, 1.5, 2.0, 2.5]  # Critical value ~π/(2r) ≈ 1.57
    
    for tau in tau_values:
        f = lambda t, N, N_d: delay_logistic(t, N, N_d, r, K)
        t, N = solve_dde_simple(f, 0.5, [0, 50], tau, dt=0.01)
        ax1.plot(t, N, linewidth=2, label=f'τ = {tau}')
    
    ax1.axhline(y=K, color='k', linestyle='--', alpha=0.5, label='K')
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('Population N', fontsize=12)
    ax1.set_title('Delay Logistic Equation: Hutchinson Model', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Stability diagram
    ax2 = axes[0, 1]
    
    tau_range = np.linspace(0, 3, 100)
    r_tau = r * tau_range
    
    ax2.fill_between(tau_range, 0, np.pi/2, alpha=0.3, color='green', label='Stable')
    ax2.fill_between(tau_range, np.pi/2, 3, alpha=0.3, color='red', label='Oscillatory')
    ax2.axhline(y=np.pi/2, color='k', linestyle='--', linewidth=2)
    ax2.text(2.5, np.pi/2 + 0.1, '$r\\tau = \\pi/2$', fontsize=12)
    
    ax2.set_xlabel('Delay τ', fontsize=12)
    ax2.set_ylabel('r τ', fontsize=12)
    ax2.set_title('Stability Diagram (r = 1)', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 3)
    ax2.set_ylim(0, 3)
    
    # Mackey-Glass equation (chaotic)
    ax3 = axes[1, 0]
    
    def mackey_glass(t, x, x_delayed, a=0.2, b=0.1, c=10, n=10):
        return a * x_delayed / (1 + x_delayed**n) - b * x
    
    tau_mg = 17  # This delay value produces chaos
    t_mg, x_mg = solve_dde_simple(mackey_glass, 1.2, [0, 1000], tau_mg, dt=0.1)
    
    # Only plot steady-state portion
    start_idx = len(t_mg) // 2
    ax3.plot(t_mg[start_idx:], x_mg[start_idx:], 'b-', linewidth=0.5)
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('x', fontsize=12)
    ax3.set_title(f'Mackey-Glass Equation (τ = {tau_mg}, chaotic)', fontsize=14, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    # Phase space reconstruction
    ax4 = axes[1, 1]
    
    # Use delay coordinate embedding
    embedding_delay = int(tau_mg / 0.1)
    x_t = x_mg[start_idx:-embedding_delay]
    x_t_tau = x_mg[start_idx + embedding_delay:]
    
    ax4.plot(x_t[::10], x_t_tau[::10], 'b-', linewidth=0.5, alpha=0.7)
    ax4.set_xlabel('$x(t)$', fontsize=12)
    ax4.set_ylabel('$x(t-\\tau)$', fontsize=12)
    ax4.set_title('Phase Space Reconstruction (Attractor)', fontsize=14, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('delay_differential_equations.png', dpi=150, bbox_inches='tight')
    plt.show()

delay_differential_equations()

Stochastic Differential Equations (SDEs)

Brownian Motion and Wiener Process

Brownian motionis a stochastic process satisfying: 1.$W(0) = 0W(t) - W(s)N(0, t-s)$3. Increments over different time intervals are independent

It ô Equation

The standard form of stochastic differential equations: -: Drift term (deterministic part) -: Diffusion term (stochastic part) -: Wiener increment,

Geometric Brownian Motion

The fundamental model in financial mathematics:Solution:This is the foundation of Black-Scholes option pricing.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
def stochastic_differential_equations():
    """Stochastic differential equation examples"""
    
    np.random.seed(42)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Brownian motion
    ax1 = axes[0, 0]
    
    n_paths = 10
    n_steps = 1000
    T = 1.0
    dt = T / n_steps
    
    for _ in range(n_paths):
        dW = np.random.normal(0, np.sqrt(dt), n_steps)
        W = np.concatenate([[0], np.cumsum(dW)])
        t = np.linspace(0, T, n_steps + 1)
        ax1.plot(t, W, linewidth=1, alpha=0.7)
    
    ax1.axhline(y=0, color='k', linewidth=0.5)
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('$W_t$', fontsize=12)
    ax1.set_title('Brownian Motion (Wiener Process)', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # Geometric Brownian motion
    ax2 = axes[0, 1]
    
    mu = 0.1   # Drift rate
    sigma = 0.3  # Volatility
    S0 = 100   # Initial price
    
    for _ in range(n_paths):
        dW = np.random.normal(0, np.sqrt(dt), n_steps)
        W = np.concatenate([[0], np.cumsum(dW)])
        t = np.linspace(0, T, n_steps + 1)
        
        # Analytical solution
        S = S0 * np.exp((mu - 0.5*sigma**2)*t + sigma*W)
        ax2.plot(t, S, linewidth=1, alpha=0.7)
    
    ax2.axhline(y=S0, color='k', linestyle='--', alpha=0.5)
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('$S_t$', fontsize=12)
    ax2.set_title(f'Geometric Brownian Motion (μ={mu}, σ={sigma})', fontsize=14, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # Ornstein-Uhlenbeck process (mean reversion)
    ax3 = axes[1, 0]
    
    theta = 1.0   # Reversion speed
    mu_ou = 0.0   # Long-term mean
    sigma_ou = 0.5
    
    def euler_maruyama_ou(X0, T, n_steps, theta, mu, sigma):
        dt = T / n_steps
        X = np.zeros(n_steps + 1)
        X[0] = X0
        
        for i in range(n_steps):
            dW = np.random.normal(0, np.sqrt(dt))
            X[i+1] = X[i] + theta*(mu - X[i])*dt + sigma*dW
        
        return X
    
    t = np.linspace(0, 5, 1000)
    
    for X0 in [-2, -1, 0, 1, 2]:
        X = euler_maruyama_ou(X0, 5, 999, theta, mu_ou, sigma_ou)
        ax3.plot(t, X, linewidth=1, alpha=0.7)
    
    ax3.axhline(y=mu_ou, color='k', linestyle='--', label='Long-term mean')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('$X_t$', fontsize=12)
    ax3.set_title('Ornstein-Uhlenbeck Process (Mean Reversion)', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Probability density evolution
    ax4 = axes[1, 1]
    
    # GBM distribution at different times
    times = [0.1, 0.5, 1.0]
    S_range = np.linspace(50, 200, 200)
    
    for t in times:
        # Log-normal distribution
        mean_log = np.log(S0) + (mu - 0.5*sigma**2)*t
        std_log = sigma * np.sqrt(t)
        
        pdf = 1/(S_range * std_log * np.sqrt(2*np.pi)) * \
              np.exp(-(np.log(S_range) - mean_log)**2 / (2*std_log**2))
        
        ax4.plot(S_range, pdf, linewidth=2, label=f't = {t}')
    
    ax4.axvline(x=S0, color='k', linestyle='--', alpha=0.5)
    ax4.set_xlabel('$S_t$', fontsize=12)
    ax4.set_ylabel('PDF', fontsize=12)
    ax4.set_title('GBM: Probability Density Evolution', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('stochastic_differential_equations.png', dpi=150, bbox_inches='tight')
    plt.show()

stochastic_differential_equations()

Fokker-Planck Equation

The probability densityof an SDE satisfies the Fokker-Planck equation (or Kolmogorov forward equation):This is a partial differential equation connecting SDEs and PDEs.

Fractional Differential Equations

Fractional Derivatives

Riemann-Liouville fractional integral:

Caputo fractional derivative:where.

Physical Meaning

Fractional derivatives have memory effects— the current state depends on the entire history, not just the instantaneous rate of change.

Applications: - Viscoelastic materials: Between purely elastic and purely viscous - Anomalous diffusion: Diffusion rate proportional to() - Porous media flow: Non-Darcy flow

Fractional Relaxation

Standard first-order relaxation:, solution is(exponential decay)

Fractional relaxation:The solution involves the Mittag-Leffler function:where 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def fractional_differential_equations():
    """Fractional differential equation examples"""
    
    from scipy.special import gamma
    
    def mittag_leffler(alpha, z, n_terms=100):
        """Compute Mittag-Leffler function E_α(z)"""
        result = 0
        for k in range(n_terms):
            result += z**k / gamma(alpha*k + 1)
        return result
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Fractional relaxation
    ax1 = axes[0, 0]
    
    t = np.linspace(0, 10, 200)
    lambda_val = 1.0
    
    # Different orders
    alphas = [0.5, 0.75, 1.0, 1.5]
    
    for alpha in alphas:
        if alpha == 1.0:
            y = np.exp(-lambda_val * t)
            label = f'α = 1 (exponential)'
        else:
            y = np.array([mittag_leffler(alpha, -lambda_val * ti**alpha) for ti in t])
            label = f'α = {alpha}'
        ax1.plot(t, y, linewidth=2, label=label)
    
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('y(t)', fontsize=12)
    ax1.set_title('Fractional Relaxation', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(0, 1.1)
    
    # Anomalous diffusion
    ax2 = axes[0, 1]
    
    # Mean square displacement <x ²> ∝ t^α
    t = np.linspace(0.1, 10, 100)
    
    for alpha in [0.5, 1.0, 1.5, 2.0]:
        msd = t**alpha
        if alpha < 1:
            label = f'α = {alpha} (subdiffusion)'
        elif alpha == 1:
            label = f'α = {alpha} (normal)'
        elif alpha < 2:
            label = f'α = {alpha} (superdiffusion)'
        else:
            label = f'α = {alpha} (ballistic)'
        ax2.plot(t, msd, linewidth=2, label=label)
    
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Mean Square Displacement', fontsize=12)
    ax2.set_title('Anomalous Diffusion', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_yscale('log')
    ax2.set_xscale('log')
    
    # Viscoelastic material response
    ax3 = axes[1, 0]
    
    # Strain response under step stress
    t = np.linspace(0, 5, 200)
    
    # Elastic
    ax3.axhline(y=1, color='b', linestyle='--', linewidth=2, label='Elastic (α=0)')
    
    # Viscous
    ax3.plot(t, t, 'r--', linewidth=2, label='Viscous (α=1)')
    
    # Fractional (viscoelastic)
    for alpha in [0.3, 0.5, 0.7]:
        strain = t**alpha
        ax3.plot(t, strain, linewidth=2, label=f'Viscoelastic (α={alpha})')
    
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Strain (normalized)', fontsize=12)
    ax3.set_title('Viscoelastic Response (Step Stress)', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    ax3.set_xlim(0, 5)
    ax3.set_ylim(0, 5)
    
    # Mittag-Leffler function
    ax4 = axes[1, 1]
    
    z = np.linspace(-5, 0, 200)
    
    for alpha in [0.5, 0.75, 1.0, 1.25, 1.5]:
        E = np.array([mittag_leffler(alpha, zi) for zi in z])
        if alpha == 1.0:
            label = f'α = 1 (exp(z))'
        else:
            label = f'α = {alpha}'
        ax4.plot(z, E, linewidth=2, label=label)
    
    ax4.set_xlabel('z', fontsize=12)
    ax4.set_ylabel('$E_α(z)$', fontsize=12)
    ax4.set_title('Mittag-Leffler Function', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_ylim(0, 1.5)
    
    plt.tight_layout()
    plt.savefig('fractional_differential_equations.png', dpi=150, bbox_inches='tight')
    plt.show()

fractional_differential_equations()

Modern Developments in Numerical Methods

Structure-Preserving Algorithms

Traditional numerical methods may destroy important properties of physical systems. Symplectic integrators preserve the phase space volume of Hamiltonian systems:

St ö rmer-Verlet algorithm:$

$

Exponential Integrators

For stiff systems, exponential integrators use matrix exponentials to exactly handle the linear part:where.

Physics-Informed Neural Networks (PINNs)

Combining neural networks with physical laws:

Loss function:

Complete Series Knowledge Review

Core Knowledge Summary

First-Order Equations

Type Standard Form Solution Method
Separable Separate variables and integrate
Linear Integrating factor
Exact , Potential function method
Bernoulli Substitution
Riccati Transform after finding particular solution

Higher-Order Linear Equations

Constant coefficient homogeneous equations: Characteristic roots determine solution form

Characteristic Root Corresponding Solution
Real root
Repeated root(-fold) ,
Complex roots ,

Nonhomogeneous equations: - Method of undetermined coefficients (special right-hand sides) - Variation of parameters (general right-hand sides) - Laplace transform method

Linear Systems

  • General solution:
  • Matrix exponential: Diagonalization or Jordan normal form
  • Phase plane analysis: Nodes, foci, saddles, centers

Stability Theory

Lyapunov stability: - Stable: Small perturbation keeps trajectory in neighborhood - Asymptotically stable: Trajectory tends toward equilibrium - Unstable: Trajectory leaves equilibrium

Determination methods: - Linearization (Hartman-Grobman theorem) - Lyapunov functions - Routh-Hurwitz criterion

Nonlinear Dynamics

Qualitative analysis: - Equilibrium classification - Limit cycles - Bifurcation theory (saddle-node, Hopf, period-doubling) - Chaos

Important theorems: - Poincar é-Bendixson theorem (planar systems) - Lyapunov exponents (chaos determination)

Problem-Solving Strategy Summary

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
Encountering an ODE problem
    │
    ├─ Is it first-order?
    │   ├─ Yes → Determine type (separable/linear/exact/Bernoulli)
    │   └─ No → Continue
    │
    ├─ Is it constant coefficient linear?
    │   ├─ Yes → Characteristic equation or Laplace transform
    │   └─ No → Continue
    │
    ├─ Are there singular points?
    │   ├─ Yes → Frobenius series solution
    │   └─ No → Continue
    │
    ├─ Is it in system form?
    │   ├─ Yes → Matrix method or phase plane analysis
    │   └─ No → Continue
    │
    └─ Numerical methods
        ├─ Non-stiff → RK45
        ├─ Stiff → Implicit methods (BDF/Radau)
        └─ Long-time integration → Symplectic methods

Common Pitfalls and Cautions

  1. Missing singular solutions: Envelopes may be solutions
  2. Implicit solutions: Not all ODEs have explicit solutions
  3. Numerical instability: Stiff systems require implicit methods
  4. Chaotic sensitivity: Tiny initial changes lead to completely different trajectories
  5. Stability determination: Linearization only valid near equilibrium
  6. Physical meaning: Mathematical solutions may not be physically realizable

Study Recommendations and Advanced Pathways

Beginner Path

  1. Solid foundation (1-2 months)
    • Master all basic types of first-order equations
    • Skillfully use integrating factors
    • Understand basic phase plane concepts
  2. Core methods (2-3 months)
    • Higher-order constant coefficient equations
    • Laplace transforms
    • Linear systems and matrix exponentials
  3. Qualitative theory (2 months)
    • Stability analysis
    • Phase portrait plotting
    • Introduction to bifurcation

Advanced Directions

Direction A: Applied Mathematics - Partial differential equations - Numerical analysis - Optimization and control

Direction B: Dynamical Systems - Ergodic theory - Advanced bifurcation theory - Chaos and complexity

Direction C: Stochastic Analysis - Stochastic differential equations - Stochastic processes - Financial mathematics

Direction D: Machine Learning - Neural ODEs - Physics-informed neural networks - Scientific computing

Textbooks: 1. Boyce & DiPrima - Elementary Differential Equations (introductory) 2. Strogatz - Nonlinear Dynamics and Chaos (nonlinear) 3. Hirsch, Smale, Devaney - Differential Equations, Dynamical Systems, and Chaos (comprehensive) 4. Ø ksendal - Stochastic Differential Equations (stochastic)

Online Courses: 1. MIT 18.03 - Differential Equations 2. MIT 18.06 - Linear Algebra (matrix method foundation) 3. Coursera - Introduction to Dynamical Systems and Chaos

Software Tools: 1. Python: scipy.integrate, diffeqpy 2. Julia: DifferentialEquations.jl (currently most powerful) 3. MATLAB: ode45, pdepe 4. Mathematica: DSolve, NDSolve

Exercises

Comprehensive Problems

  1. Neural ODE: Implement a simple Neural ODE model to learn spiral data dynamics. Compare parameter efficiency with standard neural networks.

  2. Delay effects: Consider a feedback control system with delay:Analyze stability for different delaysand design a Smith predictor for compensation.

  3. Stochastic resonance: Add noise to a double-well potential system:Explore the effect of noise intensityon transition behavior.

  4. Fractional oscillator: Analyze the fractional harmonic oscillator:Compare oscillation behavior for differentwith the standard oscillator.

  5. Physics-informed network: Use PINN to solve the boundary value problem:Analyze convergence and accuracy.

Theoretical Problems

  1. Prove that symplectic integrators preserve phase space volume (discrete version of Liouville's theorem).

  2. Explain why It ô and Stratonovich integrals give different results and when each applies.

  3. Derive the Fokker-Planck equation and explain its relationship to SDEs.

  4. Prove that the Mittag-Leffler function reduces to the exponential function when.

  5. Discuss the universal approximation properties of Neural ODEs — what kinds of functions can they represent?

Application Problems

  1. Financial options: Simulate stock prices using geometric Brownian motion and calculate Monte Carlo prices for European call options.

  2. Epidemic prediction: Introduce stochasticity into the SIR model and analyze uncertainty in epidemic prediction.

  3. Neuroscience: Implement the Hodgkin-Huxley model and analyze neural action potentials.

  4. Climate modeling: Implement the simplified Lorenz-63 model and analyze the butterfly effect and predictability timescales.

Programming Projects

  1. Implement a complete ODE solver library including:
    • Explicit methods: Euler, RK4, RK45
    • Implicit methods: Backward Euler, Crank-Nicolson
    • Symplectic methods: St ö rmer-Verlet
    • Adaptive step size control
    • Stiffness detection
  2. Build an interactive phase plane analysis tool that can:
    • Plot vector fields
    • Trace trajectories
    • Automatically identify and classify equilibria
    • Compute Lyapunov exponents
  3. Implement a Neural ODE framework supporting:
    • Arbitrary neural network architectures
    • Adjoint method backpropagation
    • Adaptive ODE solvers
    • Continuous normalizing flows

Conclusion

Differential equations are the language for describing the laws of change in nature. From Newton's celestial mechanics to modern machine learning, they have always stood at the frontier of science.

Our journey began with the simplest first-order equations, passed through the elegant framework of linear theory, traversed the complex landscape of nonlinear worlds, and finally arrived at the new continents of stochastic and fractional calculus. This is not just an accumulation of mathematical techniques but training in ways of thinking — how to see patterns in change, how to capture dynamics with equations.

But this is just the beginning. True understanding comes from application — when you model a real problem with ODEs, when you debug numerical code to make solutions converge, when you stare at a chaotic attractor thinking about the nature of predictability — in those moments, differential equations become not symbols in a textbook but tools for understanding the world.

I hope this series can be a valuable stop on your mathematical journey. Happy learning, and never stop exploring!

References

  1. Chen, R. T. Q., et al. (2018). Neural Ordinary Differential Equations. NeurIPS.
  2. Kloeden, P. E., & Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Springer.
  3. Diethelm, K. (2010). The Analysis of Fractional Differential Equations. Springer.
  4. Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric Numerical Integration. Springer.
  5. Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). Physics-informed neural networks. Journal of Computational Physics.
  6. Smith, H. (2011). An Introduction to Delay Differential Equations. Springer.

Series Complete

Thank you for reading the "Erta Ordinary Differential Equations" series. This is Chapter 18, also the final chapter of this series.

  • Post title：Ordinary Differential Equations (18): Advanced Topics and Summary
  • Post author：Chen Kai
  • Create time：2019-06-30 11:15:00
  • Post link：https://www.chenk.top/ode-chapter-18-advanced-topics-summary/
  • Copyright Notice：All articles in this blog are licensed under BY-NC-SA unless stating additionally.
 Comments