常微分方程（十八）前沿专题与总结
Chen Kai BOSS
  2019-06-30 11:15:00 2019-06-30 11:15      6.1k 字  29 分钟  

我们的微分方程之旅即将到达终点，但数学本身永无止境。本章将介绍几个活跃的前沿研究方向——神经网络 ODE 、延迟微分方程、随机微分方程和分数阶微分方程——它们正在深刻改变我们对动态系统的理解。最后，我们将回顾整个系列的核心内容，为读者提供一份完整的学习路线图。

神经网络 ODE (Neural ODEs)

从 ResNet 到连续深度

深度学习中的残差网络（ ResNet）可以写成：

当层数趋于无穷、每层变化趋于无穷小时，这变成了一个 ODE：

这就是Neural ODE的核心思想。

Neural ODE 的优势

  1. 内存效率：不需要存储所有中间状态
  2. 自适应计算：复杂输入可以用更多"层"
  3. 连续时间建模：自然处理不规则时间序列
  4. 可逆性：前向和反向计算都精确

反向传播与伴随方法

传统反向传播需要存储所有中间激活值。 Neural ODE 使用伴随方法：

定义伴随状态，它满足：

参数梯度：

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
120
121
122
123
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
import torch
import torch.nn as nn

class NeuralODEFunc(nn.Module):
    """Neural ODE 的动力学函数"""
    def __init__(self, hidden_dim=64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, 2)
        )
    
    def forward(self, t, y):
        return self.net(y)

def neural_ode_demo():
    """Neural ODE 演示：学习螺旋线动力学"""
    
    # 生成螺旋线数据
    def true_dynamics(t, y):
        """真实的螺旋线动力学"""
        return np.array([
            -0.1 * y[0] + 2.0 * y[1],
            -2.0 * y[0] - 0.1 * y[1]
        ])
    
    # 生成训练数据
    np.random.seed(42)
    t_span = np.linspace(0, 10, 100)
    y0 = np.array([2.0, 0.0])
    
    # 真实轨迹
    sol = solve_ivp(true_dynamics, [0, 10], y0, t_eval=t_span, method='RK45')
    true_y = sol.y.T
    
    # 添加噪声作为观测
    noise = np.random.normal(0, 0.1, true_y.shape)
    observed_y = true_y + noise
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 真实轨迹
    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')
    
    # 时间序列
    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 概念图
    ax3 = axes[1, 0]
    
    # ResNet（离散）
    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')
    
    # 向量场
    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()

连续归一化流

Neural ODE 的一个重要应用是连续归一化流（ Continuous Normalizing Flows）。

给定初始分布，通过 ODE 变换：

概率密度的变化满足：

这允许我们将简单分布（如高斯）变换为复杂分布，用于生成模型。

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
def continuous_normalizing_flow():
    """连续归一化流的概念演示"""
    
    # 简单示例：将高斯分布变换为双峰分布
    
    def flow_dynamics(t, z, target='bimodal'):
        """流的动力学（简化版）"""
        x, y = z
        if target == 'bimodal':
            # 将点推向两个中心
            center1 = np.array([-1.5, 0])
            center2 = np.array([1.5, 0])
            
            d1 = np.sqrt((x - center1[0])**2 + (y - center1[1])**2) + 0.1
            d2 = np.sqrt((x - center2[0])**2 + (y - center2[1])**2) + 0.1
            
            # 向较近的中心移动
            if d1 < d2:
                dx = -(x - center1[0]) / d1 * 0.5
                dy = -(y - center1[1]) / d1 * 0.5
            else:
                dx = -(x - center2[0]) / d2 * 0.5
                dy = -(y - center2[1]) / d2 * 0.5
            
            return [dx, dy]
        else:
            return [0, 0]
    
    fig, axes = plt.subplots(1, 4, figsize=(16, 4))
    
    # 初始分布：标准高斯
    np.random.seed(42)
    n_samples = 500
    z0 = np.random.randn(n_samples, 2) * 0.5
    
    times = [0, 1, 3, 6]
    
    for ax, t_end in zip(axes, times):
        if t_end == 0:
            z_t = z0
        else:
            z_t = []
            for z_init in z0:
                sol = solve_ivp(flow_dynamics, [0, t_end], z_init, method='RK45')
                z_t.append(sol.y[:, -1])
            z_t = np.array(z_t)
        
        ax.scatter(z_t[:, 0], z_t[:, 1], alpha=0.5, s=10)
        ax.set_xlim(-4, 4)
        ax.set_ylim(-3, 3)
        ax.set_xlabel('$ z_1$', fontsize=12)
        ax.set_ylabel('$ z_2$', fontsize=12)
        ax.set_title(f't = {t_end}', fontsize=14, fontweight='bold')
        ax.grid(True, alpha=0.3)
        ax.set_aspect('equal')
    
    plt.suptitle('Continuous Normalizing Flow: Gaussian → Bimodal', fontsize=16, fontweight='bold', y=1.05)
    plt.tight_layout()
    plt.savefig('continuous_normalizing_flow.png', dpi=150, bbox_inches='tight')
    plt.show()

continuous_normalizing_flow()

延迟微分方程 (Delay Differential Equations)

基本概念

许多系统的当前变化率不仅依赖于当前状态，还依赖于过去的状态：

其中 是延迟时间。

例子： - 人口动力学：成年个体的繁殖依赖于 年前出生的个体 - 控制系统：传感器和执行器的响应延迟 - 神经网络：信号传输延迟

延迟 Logistic 方程

这是 Hutchinson 方程，用于描述单一物种的延迟密度依赖增长。

惊人结果：当 时，平衡点变得不稳定，出现周期振荡！

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
120
def delay_differential_equations():
    """延迟微分方程示例"""
    
    from scipy.integrate import solve_ivp
    from collections import deque
    
    def solve_dde_simple(f, history, t_span, tau, dt=0.01):
        """简单的 DDE 求解器（欧拉法）"""
        t_start, t_end = t_span
        n_steps = int((t_end - t_start) / dt)
        n_delay = int(tau / dt)
        
        # 初始化历史
        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))
    
    # 延迟 Logistic 方程
    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]  # 临界值约为 π/(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)
    
    # 稳定性图
    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 方程（混沌）
    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  # 这个延迟值产生混沌
    t_mg, x_mg = solve_dde_simple(mackey_glass, 1.2, [0, 1000], tau_mg, dt=0.1)
    
    # 只画稳态部分
    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)
    
    # 相空间重构
    ax4 = axes[1, 1]
    
    # 使用延迟坐标嵌入
    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()

控制系统中的延迟

反馈控制器的延迟会显著影响稳定性：

其中 是反馈控制。

Smith 预测器：一种补偿延迟的控制策略，通过预测未来状态来提前采取行动。

随机微分方程 (Stochastic Differential Equations)

布朗运动与 Wiener 过程

布朗运动 是一个随机过程，满足： 1. $W(0) = 0增量W(t) - W(s)服从N(0, t-s)$3. 不同时间段的增量相互独立

It ô方程

随机微分方程的标准形式：

$$

dX_t = (X_t, t)dt + (X_t, t)dW_t $$ - ：漂移项（确定性部分） - ：扩散项（随机部分） - $ dW_t：增量， dW_t N(0, dt)$

几何布朗运动

金融数学中的基本模型：

$$

dS_t = S_t dt + S_t dW_t $$

解为：

$$

S_t = S_0 ((- )t + W_t) $$

这是 Black-Scholes 期权定价的基础。

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():
    """随机微分方程示例"""
    
    np.random.seed(42)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 布朗运动
    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)
    
    # 几何布朗运动
    ax2 = axes[0, 1]
    
    mu = 0.1   # 漂移率
    sigma = 0.3  # 波动率
    S0 = 100   # 初始价格
    
    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)
        
        # 解析解
        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 过程（均值回归）
    ax3 = axes[1, 0]
    
    theta = 1.0   # 回归速度
    mu_ou = 0.0   # 长期均值
    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)
    
    # 概率密度演化
    ax4 = axes[1, 1]
    
    # GBM 在不同时间的分布
    times = [0.1, 0.5, 1.0]
    S_range = np.linspace(50, 200, 200)
    
    for t in times:
        # 对数正态分布
        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 方程

SDE 的概率密度 满足Fokker-Planck 方程（或 Kolmogorov 前向方程）：

这是一个偏微分方程，连接了 SDE 和 PDE 。

分数阶微分方程 (Fractional Differential Equations)

分数阶导数

Riemann-Liouville 分数阶积分：

$$

I^f(t) = _0^t (t-s)^{}f(s)ds $$

Caputo 分数阶导数：

$$

D^f(t) = _0^t ds $$

其中。

物理意义

分数阶导数具有记忆效应——当前状态依赖于整个历史，而不仅仅是瞬时变化率。

应用： - 黏弹性材料：介于完全弹性和完全黏性之间 - 异常扩散：扩散速度与Missing superscript or subscript argument t^ 成正比（） - 多孔介质流动：非 Darcy 流

分数阶松弛

标准一阶松弛：，解为（指数衰减）

分数阶松弛： 解涉及Mittag-Leffler 函数：

$$

y(t) = y_0 E_(-t^) $$

其中 

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 fractional_differential_equations():
    """分数阶微分方程示例"""
    
    from scipy.special import gamma
    
    def mittag_leffler(alpha, z, n_terms=100):
        """计算 Mittag-Leffler 函数 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))
    
    # 分数阶松弛
    ax1 = axes[0, 0]
    
    t = np.linspace(0, 10, 200)
    lambda_val = 1.0
    
    # 不同阶数
    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)
    
    # 异常扩散
    ax2 = axes[0, 1]
    
    # 均方位移 <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')
    
    # 黏弹性材料响应
    ax3 = axes[1, 0]
    
    # 阶跃应力下的应变响应
    # 标准弹性：ε = σ/E（瞬时）
    # 标准黏性：ε = σ t/η（线性增长）
    # 分数阶：ε ∝ t^α
    
    t = np.linspace(0, 5, 200)
    
    # 弹性
    ax3.axhline(y=1, color='b', linestyle='--', linewidth=2, label='Elastic (α=0)')
    
    # 黏性
    ax3.plot(t, t, 'r--', linewidth=2, label='Viscous (α=1)')
    
    # 分数阶（黏弹性）
    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 函数
    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()

数值方法的现代发展

结构保持算法

传统数值方法可能破坏物理系统的重要性质。辛积分法保持哈密顿系统的相空间体积：

St ö rmer-Verlet 算法：

$$

p_{n+1/2} = p_n - V(q_n) $

$

$$ q_{n+1} = q_n + h M^{-1}p_{n+1/2} $

$

指数积分法

对于刚性系统，指数积分法利用矩阵指数精确处理线性部分：

$$

y_{n+1} = e^{hA}y_n + h_1(hA)f(y_n) $$

其中。

物理信息神经网络 (PINNs)

结合神经网络和物理定律：

损失函数：

$$

L = L_{data} + L_{physics}

L_{physics} = | {dt} - f(u_{NN}, t)|^2 $$ 

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
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
def modern_numerical_methods():
    """现代数值方法演示"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 辛积分 vs 标准 RK4（哈密顿系统）
    ax1 = axes[0, 0]
    
    # 简谐振子 H = p ²/2 + q ²/2
    def hamiltonian(q, p):
        return 0.5*p**2 + 0.5*q**2
    
    def rk4_oscillator(q0, p0, h, n_steps):
        """标准 RK4"""
        q, p = [q0], [p0]
        for _ in range(n_steps):
            q_n, p_n = q[-1], p[-1]
            
            k1q, k1p = p_n, -q_n
            k2q, k2p = p_n + 0.5*h*k1p, -(q_n + 0.5*h*k1q)
            k3q, k3p = p_n + 0.5*h*k2p, -(q_n + 0.5*h*k2q)
            k4q, k4p = p_n + h*k3p, -(q_n + h*k3q)
            
            q.append(q_n + h/6*(k1q + 2*k2q + 2*k3q + k4q))
            p.append(p_n + h/6*(k1p + 2*k2p + 2*k3p + k4p))
        
        return np.array(q), np.array(p)
    
    def symplectic_verlet(q0, p0, h, n_steps):
        """辛 Verlet 积分"""
        q, p = [q0], [p0]
        for _ in range(n_steps):
            q_n, p_n = q[-1], p[-1]
            
            p_half = p_n - 0.5*h*q_n  # grad V = q
            q_new = q_n + h*p_half
            p_new = p_half - 0.5*h*q_new
            
            q.append(q_new)
            p.append(p_new)
        
        return np.array(q), np.array(p)
    
    q0, p0 = 1.0, 0.0
    h = 0.1
    n_steps = 1000
    
    q_rk4, p_rk4 = rk4_oscillator(q0, p0, h, n_steps)
    q_sym, p_sym = symplectic_verlet(q0, p0, h, n_steps)
    
    H_rk4 = hamiltonian(q_rk4, p_rk4)
    H_sym = hamiltonian(q_sym, p_sym)
    
    t = np.arange(n_steps + 1) * h
    
    ax1.plot(t, (H_rk4 - H_rk4[0])/H_rk4[0] * 100, 'r-', linewidth=2, label='RK4')
    ax1.plot(t, (H_sym - H_sym[0])/H_sym[0] * 100, 'b-', linewidth=2, label='Symplectic Verlet')
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('Energy Error (%)', fontsize=12)
    ax1.set_title('Energy Conservation: RK4 vs Symplectic', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 相空间轨迹
    ax2 = axes[0, 1]
    
    ax2.plot(q_rk4, p_rk4, 'r-', linewidth=1, alpha=0.7, label='RK4')
    ax2.plot(q_sym, p_sym, 'b-', linewidth=1, alpha=0.7, label='Symplectic')
    
    # 真实解（圆）
    theta = np.linspace(0, 2*np.pi, 100)
    ax2.plot(np.cos(theta), np.sin(theta), 'k--', linewidth=2, label='Exact')
    
    ax2.set_xlabel('q', fontsize=12)
    ax2.set_ylabel('p', fontsize=12)
    ax2.set_title('Phase Space (Long-time Integration)', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_aspect('equal')
    
    # 刚性系统：隐式 vs 显式
    ax3 = axes[1, 0]
    
    # 刚性系统示例
    def stiff_system(t, y, eps=0.01):
        return np.array([-y[0]/eps + y[1], y[0] - y[1]])
    
    eps = 0.01
    y0 = [1, 0]
    t_span = (0, 5)
    
    # 使用不同方法
    from scipy.integrate import solve_ivp
    
    sol_explicit = solve_ivp(stiff_system, t_span, y0, method='RK45', 
                            t_eval=np.linspace(0, 5, 500), args=(eps,))
    sol_implicit = solve_ivp(stiff_system, t_span, y0, method='Radau', 
                            t_eval=np.linspace(0, 5, 500), args=(eps,))
    
    ax3.plot(sol_explicit.t, sol_explicit.y[0], 'r-', linewidth=2, label='RK45 (explicit)')
    ax3.plot(sol_implicit.t, sol_implicit.y[0], 'b--', linewidth=2, label='Radau (implicit)')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('$ y_1$', fontsize=12)
    ax3.set_title(f'Stiff System (ε = {eps})', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 函数调用次数比较
    ax4 = axes[1, 1]
    
    print(f"RK45 function evaluations: {sol_explicit.nfev}")
    print(f"Radau function evaluations: {sol_implicit.nfev}")
    
    methods = ['RK23', 'RK45', 'DOP853', 'Radau', 'BDF']
    nfev_list = []
    
    for method in methods:
        sol = solve_ivp(stiff_system, t_span, y0, method=method, args=(eps,))
        nfev_list.append(sol.nfev)
    
    colors = ['red' if 'RK' in m or 'DOP' in m else 'blue' for m in methods]
    bars = ax4.bar(methods, nfev_list, color=colors, alpha=0.7)
    ax4.set_ylabel('Function Evaluations', fontsize=12)
    ax4.set_title('Method Comparison (Stiff System)', fontsize=14, fontweight='bold')
    ax4.grid(True, alpha=0.3, axis='y')
    
    # 添加图例
    from matplotlib.patches import Patch
    legend_elements = [Patch(facecolor='red', alpha=0.7, label='Explicit'),
                      Patch(facecolor='blue', alpha=0.7, label='Implicit')]
    ax4.legend(handles=legend_elements)
    
    plt.tight_layout()
    plt.savefig('modern_numerical_methods.png', dpi=150, bbox_inches='tight')
    plt.show()

modern_numerical_methods()

全系列知识回顾

系列导航图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
第 1 章：起源与直觉
    ↓
第 2 章：一阶方程 ─────────────────────────┐
    ↓                                     │
第 3 章：高阶线性方程 ← 特征方程              │
    ↓                                     │
第 4 章：拉普拉斯变换 ← 代数化                 │ 解法分支
    ↓                                     │
第 5 章：级数解法 ← Frobenius 方法            │
    ↓                                     │
第 6 章：线性微分方程组 ← 矩阵方法 ────────────┘
    ↓
第 7-9 章：稳定性与定性理论
    ↓
第 10-12 章：分岔与混沌
    ↓
第 13 章：偏微分方程引论 → [另一个系列]
    ↓
第 14-16 章：高级专题
    ↓
第 17 章：物理与工程应用
    ↓
第 18 章：前沿专题与总结（本章）

核心知识点总结

一阶方程

类型 标准形式 解法
可分离变量 分离变量后积分
线性 $ y' + P(x)y = Q(x)积分因子 e^{P dx}$
恰当 ， 势函数法
Bernoulli $ y' + P(x)y = Q(x)y^n代换 v = y^{1-n}$
Riccati 已知特解后变换

高阶线性方程

常系数齐次方程：特征根决定解的形式

特征根 对应解
实根 $ e^{rx}重根 r（ k$重）

非齐次方程： - 待定系数法（特殊右端项） - 参数变动法（一般右端项） - 拉普拉斯变换法

线性系统

- 通解： - 矩阵指数：对角化或 Jordan 标准形 - 相平面分析：节点、焦点、鞍点、中心

稳定性理论

Lyapunov 稳定性： - 稳定：小扰动后轨迹保持在邻域内 - 渐近稳定：轨迹趋向平衡点 - 不稳定：轨迹远离平衡点

判定方法： - 线性化（ Hartman-Grobman 定理） - Lyapunov 函数 - Routh-Hurwitz 准则

非线性动力学

定性分析： - 平衡点分类 - 极限环 - 分岔理论（鞍结点、 Hopf 、倍周期） - 混沌

重要定理： - Poincar é-Bendixson 定理（平面系统） - Lyapunov 指数（混沌判定）

解题策略总结

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
遇到 ODE 问题
    │
    ├─ 是否是一阶？
    │   ├─ 是 → 判断类型（可分离/线性/恰当/Bernoulli）
    │   └─ 否 → 继续
    │
    ├─ 是否是常系数线性？
    │   ├─ 是 → 特征方程法或拉普拉斯变换
    │   └─ 否 → 继续
    │
    ├─ 是否有奇点？
    │   ├─ 是 → Frobenius 级数解
    │   └─ 否 → 继续
    │
    ├─ 是否是系统形式？
    │   ├─ 是 → 矩阵方法或相平面分析
    │   └─ 否 → 继续
    │
    └─ 数值方法
        ├─ 非刚性 → RK45
        ├─ 刚性 → 隐式方法（ BDF/Radau）
        └─ 长时间积分 → 辛方法

常见陷阱与注意事项

  1. 奇解的遗漏：包络线可能是解
  2. 隐式解：不是所有 ODE 都有显式解
  3. 数值不稳定：刚性系统需要隐式方法
  4. 混沌敏感性：初值微小变化导致轨迹完全不同
  5. 稳定性判定：线性化只在平衡点邻域有效
  6. 物理意义：数学解不一定物理可实现

学习建议与进阶路线

初学者路线

  1. 扎实基础（ 1-2 个月）
    • 掌握一阶方程的所有基本类型
    • 熟练使用积分因子
    • 理解相平面的基本概念
  2. 核心方法（ 2-3 个月）
    • 高阶常系数方程
    • 拉普拉斯变换
    • 线性系统与矩阵指数
  3. 定性理论（ 2 个月）
    • 稳定性分析
    • 相图绘制
    • 分岔入门

进阶方向

方向 A：应用数学 - 偏微分方程 - 数值分析 - 优化与控制

方向 B：动力系统 - 遍历理论 - 分岔理论深化 - 混沌与复杂性

方向 C：随机分析 - 随机微分方程 - 随机过程 - 金融数学

方向 D：机器学习 - 神经网络 ODE - 物理信息神经网络 - 科学计算

推荐资源

教材： 1. Boyce & DiPrima - Elementary Differential Equations（入门） 2. Strogatz - Nonlinear Dynamics and Chaos（非线性） 3. Hirsch, Smale, Devaney - Differential Equations, Dynamical Systems, and Chaos（综合） 4. Ø ksendal - Stochastic Differential Equations（随机）

在线课程： 1. MIT 18.03 - Differential Equations 2. MIT 18.06 - Linear Algebra（矩阵方法基础） 3. Coursera - Introduction to Dynamical Systems and Chaos

软件工具： 1. Python： scipy.integrate, diffeqpy 2. Julia： DifferentialEquations.jl（目前最强大） 3. MATLAB： ode45, pdepe 4. Mathematica： DSolve, NDSolve

练习题

综合题

  1. 神经网络 ODE：实现一个简单的 Neural ODE 模型，学习螺旋线数据的动力学。比较与标准神经网络的参数效率。

  2. 延迟效应：考虑带延迟的反馈控制系统：

分析不同延迟 下的稳定性，设计 Smith 预测器进行补偿。

  1. 随机共振：在双井势系统中加入噪声： $$

dx = (x - x^3)dt + D , dW$$

探索噪声强度 对跃迁行为的影响。

  1. 分数阶振子：分析分数阶谐振子： $$

D^x + ^2 x = 0, < < 2$$

比较不同 下的振荡行为与标准振子的差异。

  1. 物理信息网络：使用 PINN 求解边值问题： $$

y'' + y = 0, y(0) = 0, y(/2) = 1$$

分析收敛性和精度。

理论题

  1. 证明辛积分器保持相空间体积（ Liouville 定理的离散版本）。

  2. 解释为什么 It ô积分与 Stratonovich 积分给出不同结果，各自适用于什么情况。

  3. 推导 Fokker-Planck 方程，并解释其与 SDE 的关系。

  4. 证明 Mittag-Leffler 函数在 时退化为指数函数。

  5. 讨论 Neural ODE 的通用近似性质——它能表示什么样的函数？

应用题

  1. 金融期权：使用几何布朗运动模拟股票价格，计算欧式看涨期权的蒙特卡洛价格。

  2. 流行病预测：在 SIR 模型中引入随机性，分析疫情预测的不确定性。

  3. 神经科学：实现 Hodgkin-Huxley 模型，分析神经元的动作电位。

  4. 气候模型：实现简化的 Lorenz-63 模型，分析蝴蝶效应和可预测性时间尺度。

编程项目

  1. 实现一个完整的 ODE 求解器库，包括：
    • 显式方法： Euler, RK4, RK45
    • 隐式方法：后向 Euler, Crank-Nicolson
    • 辛方法： St ö rmer-Verlet
    • 自适应步长控制
    • 刚性检测
  2. 构建一个交互式的相平面分析工具，可以：
    • 绘制向量场
    • 追踪轨迹
    • 自动识别平衡点和分类
    • 计算 Lyapunov 指数
  3. 实现 Neural ODE 框架，支持：
    • 任意神经网络结构
    • 伴随方法反向传播
    • 自适应 ODE 求解器
    • 连续归一化流

结语

微分方程是描述自然界变化规律的语言。从牛顿时代的天体力学，到现代的机器学习，它始终站在科学的前沿。

我们的旅程从最简单的一阶方程开始，经过线性理论的优美框架，穿越非线性世界的复杂地貌，最终抵达随机和分数阶的新大陆。这不仅是数学技术的积累，更是思维方式的训练——如何从变化中看到规律，如何用方程捕捉动态。

但这只是开始。真正的理解来自于应用——当你用 ODE 建模一个实际问题，当你调试数值代码让解收敛，当你盯着混沌吸引子思考可预测性的本质——那些时刻，微分方程不再是书本上的符号，而成为你理解世界的工具。

希望这个系列能成为你数学旅途中有价值的一站。祝学习愉快，探索不止！

参考资料

  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.

本系列完结

感谢您阅读《常微分方程的世界》系列。这是第 18 章，也是本系列的最后一章。

  • 本文标题：常微分方程（十八）前沿专题与总结
  • 本文作者：Chen Kai
  • 创建时间：2019-06-30 11:15:00
  • 本文链接：https://www.chenk.top/%E5%B8%B8%E5%BE%AE%E5%88%86%E6%96%B9%E7%A8%8B%EF%BC%88%E5%8D%81%E5%85%AB%EF%BC%89%E5%89%8D%E6%B2%BF%E4%B8%93%E9%A2%98%E4%B8%8E%E6%80%BB%E7%BB%93/
  • 版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
 评论