常微分方程（十七）物理与工程应用

微分方程不是纯数学的游戏——它是理解物理世界的语言。从天体运动到电路响应，从弹簧振动到化学反应，几乎所有动态系统的行为都可以用微分方程来描述。本章我们将深入探索常微分方程在物理和工程中的核心应用，建立从数学到实践的桥梁。

经典力学中的微分方程

牛顿运动定律的数学表述

牛顿第二定律是力学的基石，但它本质上是一个微分方程：

m {dt^2} = (, , t) $$

这里是位置向量，力可能依赖于位置、速度和时间。

重要观察：知道了力的表达式，物体的运动完全由初始位置和初始速度决定。这体现了牛顿力学的决定论本质。

自由落体与空气阻力

考虑一个从高处下落的物体，同时受到重力和空气阻力：

m = mg - kv^2 $$

其中是阻力系数，是速度（向下为正）。

终端速度：当时，加速度为零，此时达到终端速度：

v_{terminal} = $$

求解过程：这是一个可分离变量的一阶 ODE：

利用部分分式分解：

积分后得到：

v(t) = v_{terminal} ( ) $$ （假设初始速度为零）

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

def falling_with_drag():
    """带空气阻力的自由落体"""
    m = 70  # 质量 kg（人体）
    g = 9.8  # 重力加速度
    k = 0.25  # 阻力系数（与体型和姿势有关）
    
    v_terminal = np.sqrt(m * g / k)
    print(f"终端速度: {v_terminal:.1f} m/s ({v_terminal * 3.6:.1f} km/h)")
    
    def drag_ode(v, t):
        return g - (k/m) * v**2
    
    t = np.linspace(0, 30, 500)
    v0 = 0
    
    # 数值解
    v_numerical = odeint(drag_ode, v0, t).flatten()
    
    # 解析解
    v_analytical = v_terminal * np.tanh(g * t / v_terminal)
    
    # 位置（积分速度）
    x_numerical = np.cumsum(v_numerical) * (t[1] - t[0])
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # 速度-时间
    ax1 = axes[0]
    ax1.plot(t, v_numerical, 'b-', linewidth=2, label='Numerical')
    ax1.plot(t, v_analytical, 'r--', linewidth=2, label='Analytical')
    ax1.axhline(y=v_terminal, color='k', linestyle=':', label=f'Terminal: {v_terminal:.1f} m/s')
    ax1.set_xlabel('Time (s)', fontsize=12)
    ax1.set_ylabel('Velocity (m/s)', fontsize=12)
    ax1.set_title('Velocity vs Time', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 位置-时间
    ax2 = axes[1]
    ax2.plot(t, x_numerical, 'g-', linewidth=2)
    ax2.set_xlabel('Time (s)', fontsize=12)
    ax2.set_ylabel('Distance (m)', fontsize=12)
    ax2.set_title('Distance vs Time', fontsize=14, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    
    # 不同阻力系数的比较
    ax3 = axes[2]
    k_values = [0.1, 0.25, 0.5, 1.0]
    for k_val in k_values:
        v_term = np.sqrt(m * g / k_val)
        v = v_term * np.tanh(g * t / v_term)
        ax3.plot(t, v, linewidth=2, label=f'k = {k_val}')
    
    ax3.set_xlabel('Time (s)', fontsize=12)
    ax3.set_ylabel('Velocity (m/s)', fontsize=12)
    ax3.set_title('Effect of Drag Coefficient', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('falling_with_drag.png', dpi=150, bbox_inches='tight')
    plt.show()

falling_with_drag()

行星运动与开普勒问题

两体问题是经典力学中最重要的问题之一。考虑质量为的行星绕质量为的恒星运动（）：

m {dt^2} = - $$

引入极坐标，并利用角动量守恒，可得到Binet 方程：

其中。

这是一个简谐振子方程！通解为：

u = = (1 + e(- _0)) $$

即圆锥曲线的极坐标方程，是离心率： - ：椭圆（束缚轨道） - ：抛物线（临界逃逸） - ：双曲线（逃逸轨道）

def kepler_orbits():
    """开普勒轨道的可视化"""
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 不同离心率的轨道
    ax1 = axes[0]
    theta = np.linspace(0, 2*np.pi, 1000)
    
    eccentricities = [0, 0.3, 0.6, 0.9]
    p = 1  # 半正焦距
    
    for e in eccentricities:
        r = p / (1 + e * np.cos(theta))
        # 只画束缚部分
        valid = r > 0
        x = r[valid] * np.cos(theta[valid])
        y = r[valid] * np.sin(theta[valid])
        ax1.plot(x, y, linewidth=2, label=f'e = {e}')
    
    # 抛物线和双曲线
    for e, style in [(1.0, '--'), (1.5, ':')]:
        theta_range = np.linspace(-np.pi*0.8, np.pi*0.8, 500)
        r = p / (1 + e * np.cos(theta_range))
        valid = r > 0
        x = r[valid] * np.cos(theta_range[valid])
        y = r[valid] * np.sin(theta_range[valid])
        ax1.plot(x, y, linestyle=style, linewidth=2, label=f'e = {e}')
    
    ax1.plot(0, 0, 'yo', markersize=15, label='Sun')
    ax1.set_xlim(-4, 2)
    ax1.set_ylim(-3, 3)
    ax1.set_aspect('equal')
    ax1.set_xlabel('x', fontsize=12)
    ax1.set_ylabel('y', fontsize=12)
    ax1.set_title('Kepler Orbits: Different Eccentricities', fontsize=14, fontweight='bold')
    ax1.legend(loc='upper right')
    ax1.grid(True, alpha=0.3)
    
    # 数值模拟轨道
    ax2 = axes[1]
    
    def kepler_ode(t, state, GM):
        x, y, vx, vy = state
        r = np.sqrt(x**2 + y**2)
        ax = -GM * x / r**3
        ay = -GM * y / r**3
        return [vx, vy, ax, ay]
    
    GM = 1.0
    
    # 不同初始条件
    initial_conditions = [
        ([1, 0, 0, 1.0], 'Circle'),
        ([1, 0, 0, 0.8], 'Ellipse'),
        ([1, 0, 0, 1.414], 'Parabola-like'),
    ]
    
    for (state0, name) in initial_conditions:
        t_span = (0, 20)
        t_eval = np.linspace(0, 20, 2000)
        
        sol = solve_ivp(kepler_ode, t_span, state0, args=(GM,), 
                       t_eval=t_eval, method='DOP853', rtol=1e-10)
        
        ax2.plot(sol.y[0], sol.y[1], linewidth=2, label=name)
    
    ax2.plot(0, 0, 'yo', markersize=15)
    ax2.set_xlim(-3, 2)
    ax2.set_ylim(-2, 2)
    ax2.set_aspect('equal')
    ax2.set_xlabel('x', fontsize=12)
    ax2.set_ylabel('y', fontsize=12)
    ax2.set_title('Numerical Orbit Simulation', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('kepler_orbits.png', dpi=150, bbox_inches='tight')
    plt.show()

kepler_orbits()

振动系统的深入分析

多自由度振动系统

实际的振动系统往往有多个自由度。考虑个质量通过弹簧连接：

M + C + K = (t) $$

其中、、分别是质量矩阵、阻尼矩阵和刚度矩阵。

模态分析：无阻尼自由振动假设解的形式为，代入得特征值问题：

特征值对应固有频率，特征向量对应振型。

三自由度弹簧系统

def multi_dof_vibration():
    """三自由度振动系统的模态分析"""
    from scipy.linalg import eigh
    
    # 质量矩阵（假设单位质量）
    M = np.array([[1, 0, 0],
                  [0, 1, 0],
                  [0, 0, 1]])
    
    # 刚度矩阵（三个质量由弹簧连接，两端固定）
    k = 1.0
    K = np.array([[2*k, -k, 0],
                  [-k, 2*k, -k],
                  [0, -k, 2*k]])
    
    # 求解广义特征值问题
    eigenvalues, eigenvectors = eigh(K, M)
    frequencies = np.sqrt(eigenvalues)
    
    print("固有频率:", frequencies)
    print("\n 振型（列向量）:")
    print(eigenvectors)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 绘制振型
    ax1 = axes[0, 0]
    x_pos = [0, 1, 2, 3]  # 平衡位置
    
    for i, mode in enumerate(eigenvectors.T):
        # 归一化振型
        mode_normalized = mode / np.max(np.abs(mode))
        displacement = np.concatenate([[0], mode_normalized, [0]])  # 边界固定
        ax1.plot(x_pos + [4], displacement, 'o-', linewidth=2, markersize=10,
                label=f'Mode {i+1}: ω = {frequencies[i]:.3f}')
    
    ax1.axhline(y=0, color='k', linewidth=0.5)
    ax1.set_xlabel('Mass position', fontsize=12)
    ax1.set_ylabel('Displacement (normalized)', fontsize=12)
    ax1.set_title('Mode Shapes', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 时域响应（激励第一个质量）
    ax2 = axes[0, 1]
    
    def three_dof_system(t, state):
        x = state[:3]
        v = state[3:]
        a = -np.linalg.solve(M, K @ x)
        return np.concatenate([v, a])
    
    # 初始条件：只位移第一个质量
    state0 = [1, 0, 0, 0, 0, 0]
    t_span = (0, 50)
    t_eval = np.linspace(0, 50, 2000)
    
    sol = solve_ivp(three_dof_system, t_span, state0, t_eval=t_eval, method='RK45')
    
    for i in range(3):
        ax2.plot(sol.t, sol.y[i], linewidth=1.5, label=f'Mass {i+1}')
    
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Displacement', fontsize=12)
    ax2.set_title('Free Vibration Response', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 能量转移
    ax3 = axes[1, 0]
    
    # 计算每个质量的动能
    KE = np.zeros((3, len(sol.t)))
    for i in range(3):
        KE[i] = 0.5 * M[i, i] * sol.y[i+3]**2
    
    for i in range(3):
        ax3.plot(sol.t, KE[i], linewidth=1.5, label=f'Mass {i+1}')
    
    ax3.plot(sol.t, np.sum(KE, axis=0), 'k--', linewidth=2, label='Total KE')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Kinetic Energy', fontsize=12)
    ax3.set_title('Energy Distribution', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 相空间轨迹（第一个质量）
    ax4 = axes[1, 1]
    ax4.plot(sol.y[0], sol.y[3], 'b-', linewidth=1, alpha=0.7)
    ax4.plot(sol.y[0][0], sol.y[3][0], 'go', markersize=10, label='Start')
    ax4.set_xlabel('$ x_1$', fontsize=12)
    ax4.set_ylabel('$ v_1$', fontsize=12)
    ax4.set_title('Phase Portrait (Mass 1)', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('multi_dof_vibration.png', dpi=150, bbox_inches='tight')
    plt.show()

multi_dof_vibration()

非线性振动： Duffing 振子

当位移较大时，弹簧可能表现出非线性：

这就是著名的Duffing 方程。

：硬化弹簧（越拉越硬）
：软化弹簧（越拉越软）

Duffing 系统展现出丰富的动力学行为： - 多个稳态解（跳跃现象） - 次谐波和超谐波共振 - 混沌运动

def duffing_oscillator():
    """Duffing 振子的动力学"""
    
    def duffing(t, state, delta, alpha, beta, gamma, omega):
        x, v = state
        return [v, -delta*v - alpha*x - beta*x**3 + gamma*np.cos(omega*t)]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 参数
    delta = 0.3  # 阻尼
    alpha = -1.0  # 线性项（双井势）
    beta = 1.0   # 非线性项
    omega = 1.2  # 驱动频率
    
    # 不同驱动幅度
    gamma_values = [0.2, 0.3, 0.37, 0.5]
    
    for ax, gamma in zip(axes.flatten(), gamma_values):
        t_span = (0, 500)
        t_eval = np.linspace(400, 500, 2000)  # 只看稳态
        
        sol = solve_ivp(duffing, t_span, [0.1, 0], args=(delta, alpha, beta, gamma, omega),
                       t_eval=t_eval, method='RK45')
        
        ax.plot(sol.y[0], sol.y[1], 'b-', linewidth=0.5, alpha=0.7)
        ax.set_xlabel('x', fontsize=12)
        ax.set_ylabel('v', fontsize=12)
        ax.set_title(f'γ = {gamma}', fontsize=14, fontweight='bold')
        ax.grid(True, alpha=0.3)
    
    plt.suptitle('Duffing Oscillator: Transition to Chaos', fontsize=16, fontweight='bold', y=1.02)
    plt.tight_layout()
    plt.savefig('duffing_oscillator.png', dpi=150, bbox_inches='tight')
    plt.show()

duffing_oscillator()

电路与电磁学

一般 RLC 网络

复杂的电路网络可以用微分方程组描述。考虑状态空间方法：

选择电容电压和电感电流作为状态变量，则：

C = i_C $

利用基尔霍夫定律建立和与状态变量的关系，得到：

其中，是输入（电压源或电流源）。

耦合电路：变压器

变压器通过磁耦合连接两个电路：

v_1 = L_1 + M $

其中是互感系数，满足。

def coupled_circuits():
    """耦合电路（变压器）的分析"""
    
    # 参数
    L1 = 1e-3   # 初级电感 1mH
    L2 = 4e-3   # 次级电感 4mH
    M = 1.5e-3  # 互感
    R1 = 10     # 初级电阻
    R2 = 100    # 次级电阻（负载）
    
    # 耦合系数
    k = M / np.sqrt(L1 * L2)
    print(f"耦合系数 k = {k:.2f}")
    
    # 匝比
    n = np.sqrt(L2 / L1)
    print(f"匝比 n = {n:.1f}")
    
    def transformer(t, state, V0, omega):
        i1, i2 = state
        # v1 = V0*sin(ω t) = L1*di1/dt + M*di2/dt + R1*i1
        # 0 = M*di1/dt + L2*di2/dt + R2*i2
        
        det = L1*L2 - M**2
        v1 = V0 * np.sin(omega * t)
        
        di1_dt = (L2*(v1 - R1*i1) + M*R2*i2) / det
        di2_dt = (-M*(v1 - R1*i1) - L1*R2*i2) / det
        
        return [di1_dt, di2_dt]
    
    V0 = 10  # 输入电压幅值
    omega = 2*np.pi*1000  # 1kHz
    
    t_span = (0, 0.02)  # 20ms
    t_eval = np.linspace(0, 0.02, 2000)
    
    sol = solve_ivp(transformer, t_span, [0, 0], args=(V0, omega),
                   t_eval=t_eval, method='RK45')
    
    # 计算电压
    v1 = V0 * np.sin(omega * sol.t)
    v2 = R2 * sol.y[1]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 电流
    ax1 = axes[0, 0]
    ax1.plot(sol.t*1000, sol.y[0]*1000, 'b-', linewidth=2, label='$ i_1$(primary)')
    ax1.plot(sol.t*1000, sol.y[1]*1000, 'r-', linewidth=2, label='$ i_2$(secondary)')
    ax1.set_xlabel('Time (ms)', fontsize=12)
    ax1.set_ylabel('Current (mA)', fontsize=12)
    ax1.set_title('Transformer Currents', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 电压
    ax2 = axes[0, 1]
    ax2.plot(sol.t*1000, v1, 'b-', linewidth=2, label='$ v_1$(input)')
    ax2.plot(sol.t*1000, v2, 'r-', linewidth=2, label='$ v_2$(output)')
    ax2.set_xlabel('Time (ms)', fontsize=12)
    ax2.set_ylabel('Voltage (V)', fontsize=12)
    ax2.set_title('Transformer Voltages', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 功率
    ax3 = axes[1, 0]
    p1 = v1 * sol.y[0]  # 输入功率
    p2 = v2 * sol.y[1]  # 输出功率
    ax3.plot(sol.t*1000, p1*1000, 'b-', linewidth=2, label='Input power')
    ax3.plot(sol.t*1000, p2*1000, 'r-', linewidth=2, label='Output power')
    ax3.set_xlabel('Time (ms)', fontsize=12)
    ax3.set_ylabel('Power (mW)', fontsize=12)
    ax3.set_title('Power Transfer', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 频率响应
    ax4 = axes[1, 1]
    freqs = np.logspace(2, 5, 200)
    gain = []
    
    for f in freqs:
        omega_i = 2*np.pi*f
        # 稳态分析（使用相量法）
        # 简化：假设稳态，计算增益
        s = 1j*omega_i
        Z11 = R1 + s*L1
        Z12 = s*M
        Z21 = s*M
        Z22 = R2 + s*L2
        
        # 输出电压/输入电压
        H = -Z21*R2 / (Z11*Z22 - Z12*Z21)
        gain.append(np.abs(H))
    
    ax4.semilogx(freqs, 20*np.log10(gain), 'g-', linewidth=2)
    ax4.set_xlabel('Frequency (Hz)', fontsize=12)
    ax4.set_ylabel('Gain (dB)', fontsize=12)
    ax4.set_title('Frequency Response', fontsize=14, fontweight='bold')
    ax4.grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('coupled_circuits.png', dpi=150, bbox_inches='tight')
    plt.show()

coupled_circuits()

非线性电路：二极管和晶体管

实际电路元件往往是非线性的。例如，二极管的 I-V 特性：

i = I_s(e^{v/V_T} - 1) $$

其中是反向饱和电流，是热电压。

包含二极管的电路方程：

这类非线性 ODE 通常只能数值求解。

热传导与扩散

牛顿冷却定律

物体与环境的热交换遵循：

解为：

T(t) = T_{env} + (T_0 - T_{env})e^{-kt} $$

应用：法医学中估计死亡时间。

多区域热传导

考虑两个通过界面接触的物体：

C_1 = -h(T_1 - T_2) $

这是一个耦合的线性系统。

def thermal_conduction():
    """热传导问题"""
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # 简单牛顿冷却
    ax1 = axes[0]
    T_env = 20  # 环境温度
    T0 = 90     # 初始温度
    k_values = [0.05, 0.1, 0.2, 0.5]
    t = np.linspace(0, 60, 500)
    
    for k in k_values:
        T = T_env + (T0 - T_env) * np.exp(-k * t)
        ax1.plot(t, T, linewidth=2, label=f'k = {k}')
    
    ax1.axhline(y=T_env, color='k', linestyle='--', label='Environment')
    ax1.set_xlabel('Time (min)', fontsize=12)
    ax1.set_ylabel('Temperature (° C)', fontsize=12)
    ax1.set_title('Newton\'s Law of Cooling', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 两区域热传导
    ax2 = axes[1]
    
    def two_region(t, state):
        T1, T2 = state
        C1, C2 = 100, 50  # 热容
        h = 5    # 界面传热系数
        k = 1    # 与环境传热系数
        T_env = 20
        
        dT1 = -h/C1 * (T1 - T2)
        dT2 = h/C2 * (T1 - T2) - k/C2 * (T2 - T_env)
        return [dT1, dT2]
    
    sol = solve_ivp(two_region, [0, 120], [100, 50], t_eval=np.linspace(0, 120, 500))
    
    ax2.plot(sol.t, sol.y[0], 'r-', linewidth=2, label='Region 1')
    ax2.plot(sol.t, sol.y[1], 'b-', linewidth=2, label='Region 2')
    ax2.axhline(y=20, color='k', linestyle='--', label='Environment')
    ax2.set_xlabel('Time (min)', fontsize=12)
    ax2.set_ylabel('Temperature (° C)', fontsize=12)
    ax2.set_title('Two-Region Heat Transfer', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 法医学应用：死亡时间估计
    ax3 = axes[2]
    
    # 逆向计算
    T_body = 37  # 正常体温
    T_env = 20
    T_measured = 30  # 测量温度
    k = 0.1  # 假设的冷却常数
    
    # t = -ln((T_measured - T_env)/(T_body - T_env))/k
    t_death = -np.log((T_measured - T_env)/(T_body - T_env))/k
    print(f"估计死亡时间：{t_death:.1f}小时前")
    
    t = np.linspace(-t_death-2, 10, 500)
    T = T_env + (T_body - T_env) * np.exp(-k * (t + t_death))
    
    ax3.plot(t, T, 'g-', linewidth=2)
    ax3.axvline(x=0, color='r', linestyle='--', label='Measurement time')
    ax3.axvline(x=-t_death, color='b', linestyle='--', label=f'Death time (~{t_death:.1f}h ago)')
    ax3.axhline(y=T_measured, color='gray', linestyle=':', alpha=0.5)
    ax3.plot(0, T_measured, 'ro', markersize=10)
    ax3.set_xlabel('Time (hours)', fontsize=12)
    ax3.set_ylabel('Body Temperature (° C)', fontsize=12)
    ax3.set_title('Forensic Application', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('thermal_conduction.png', dpi=150, bbox_inches='tight')
    plt.show()

thermal_conduction()

化学反应动力学

一级反应

放射性衰变和许多化学反应遵循一级动力学：

解： 半衰期：

连续反应

考虑：

解： - - -

可逆反应与化学平衡

：

平衡常数。

def chemical_kinetics():
    """化学反应动力学"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 一级反应（半衰期）
    ax1 = axes[0, 0]
    t = np.linspace(0, 10, 500)
    
    half_lives = [1, 2, 4]
    for t_half in half_lives:
        k = np.log(2) / t_half
        N = np.exp(-k * t)
        ax1.plot(t, N, linewidth=2, label=f'$ t_{{ 1/2 }}$= {t_half}')
    
    ax1.axhline(y=0.5, color='k', linestyle='--', alpha=0.5)
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('$N/N_0$', fontsize=12)
    ax1.set_title('First-Order Kinetics', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 连续反应
    ax2 = axes[0, 1]
    
    k1, k2 = 1.0, 0.5
    A0 = 1.0
    
    A = A0 * np.exp(-k1 * t)
    B = k1 * A0 / (k2 - k1) * (np.exp(-k1 * t) - np.exp(-k2 * t))
    C = A0 - A - B
    
    ax2.plot(t, A, 'b-', linewidth=2, label='[A]')
    ax2.plot(t, B, 'g-', linewidth=2, label='[B]')
    ax2.plot(t, C, 'r-', linewidth=2, label='[C]')
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Concentration', fontsize=12)
    ax2.set_title('Consecutive Reactions: A → B → C', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 可逆反应
    ax3 = axes[1, 0]
    
    def reversible(t, state, kf, kr):
        A, B = state
        return [-kf*A + kr*B, kf*A - kr*B]
    
    kf, kr = 1.0, 0.3
    K_eq = kf / kr
    A_eq = 1 / (1 + K_eq)
    B_eq = K_eq / (1 + K_eq)
    
    sol = solve_ivp(reversible, [0, 10], [1, 0], args=(kf, kr),
                   t_eval=np.linspace(0, 10, 500))
    
    ax3.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='[A]')
    ax3.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='[B]')
    ax3.axhline(y=A_eq, color='b', linestyle='--', alpha=0.5)
    ax3.axhline(y=B_eq, color='r', linestyle='--', alpha=0.5)
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Concentration', fontsize=12)
    ax3.set_title(f'Reversible Reaction: $K_{{ eq }}$= {K_eq:.2f}', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # Michaelis-Menten 酶动力学
    ax4 = axes[1, 1]
    
    def michaelis_menten(t, state, Vmax, Km):
        S = state[0]
        return [-Vmax * S / (Km + S)]
    
    Vmax = 1.0
    Km_values = [0.1, 0.5, 2.0]
    
    for Km in Km_values:
        sol = solve_ivp(michaelis_menten, [0, 10], [1.0], args=(Vmax, Km),
                       t_eval=np.linspace(0, 10, 500))
        ax4.plot(sol.t, sol.y[0], linewidth=2, label=f'$K_m$= {Km}')
    
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('[S]', fontsize=12)
    ax4.set_title('Michaelis-Menten Kinetics', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('chemical_kinetics.png', dpi=150, bbox_inches='tight')
    plt.show()

chemical_kinetics()

生物系统与种群动力学

Lotka-Volterra 捕食者-猎物模型

- ：猎物数量 - ：捕食者数量 - ：猎物增长率 - ：捕食率 - ：捕食者死亡率 - ：捕食者转化效率

平衡点分析：

：灭绝（不稳定鞍点）
：共存（中心，周期解）

竞争排斥原理

两个物种竞争同一资源：

竞争系数决定了共存还是排斥。

def population_dynamics():
    """种群动力学模型"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Lotka-Volterra
    ax1 = axes[0, 0]
    
    def lotka_volterra(t, state, alpha, beta, delta, gamma):
        x, y = state
        return [alpha*x - beta*x*y, delta*x*y - gamma*y]
    
    alpha, beta, delta, gamma = 1.0, 0.1, 0.05, 0.5
    
    t_span = (0, 100)
    t_eval = np.linspace(0, 100, 2000)
    
    for x0, y0 in [(10, 2), (20, 5), (5, 10)]:
        sol = solve_ivp(lotka_volterra, t_span, [x0, y0], 
                       args=(alpha, beta, delta, gamma), t_eval=t_eval)
        ax1.plot(sol.y[0], sol.y[1], linewidth=1.5)
    
    # 平衡点
    x_eq, y_eq = gamma/delta, alpha/beta
    ax1.plot(x_eq, y_eq, 'r*', markersize=15, label='Equilibrium')
    
    ax1.set_xlabel('Prey (x)', fontsize=12)
    ax1.set_ylabel('Predator (y)', fontsize=12)
    ax1.set_title('Lotka-Volterra Phase Portrait', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 时间演化
    ax2 = axes[0, 1]
    sol = solve_ivp(lotka_volterra, t_span, [10, 2], 
                   args=(alpha, beta, delta, gamma), t_eval=t_eval)
    
    ax2.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Prey')
    ax2.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Predator')
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Population', fontsize=12)
    ax2.set_title('Population Dynamics', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 竞争模型
    ax3 = axes[1, 0]
    
    def competition(t, state, r1, r2, K1, K2, a12, a21):
        x, y = state
        return [r1*x*(1 - (x + a12*y)/K1),
                r2*y*(1 - (y + a21*x)/K2)]
    
    # 共存情况
    r1, r2 = 1.0, 1.0
    K1, K2 = 100, 100
    a12, a21 = 0.5, 0.5
    
    sol = solve_ivp(competition, [0, 50], [10, 10], 
                   args=(r1, r2, K1, K2, a12, a21),
                   t_eval=np.linspace(0, 50, 500))
    
    ax3.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Species 1')
    ax3.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Species 2')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Population', fontsize=12)
    ax3.set_title('Competition: Coexistence', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 排斥情况
    ax4 = axes[1, 1]
    a12, a21 = 1.5, 1.5
    
    sol = solve_ivp(competition, [0, 50], [10, 9], 
                   args=(r1, r2, K1, K2, a12, a21),
                   t_eval=np.linspace(0, 50, 500))
    
    ax4.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Species 1')
    ax4.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Species 2')
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('Population', fontsize=12)
    ax4.set_title('Competition: Exclusion', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('population_dynamics.png', dpi=150, bbox_inches='tight')
    plt.show()

population_dynamics()

传染病模型： SIR

基本再生数 决定疫情是否爆发。

def sir_model():
    """SIR 传染病模型"""
    
    def sir(t, state, beta, gamma):
        S, I, R = state
        N = S + I + R
        return [-beta*S*I/N, beta*S*I/N - gamma*I, gamma*I]
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # 基本 SIR 动力学
    ax1 = axes[0]
    
    beta = 0.3   # 传播率
    gamma = 0.1  # 恢复率
    N = 1000
    I0 = 1
    R0 = 0
    S0 = N - I0 - R0
    
    R_0 = beta * S0 / (gamma * N)
    print(f"基本再生数 R ₀ = {R_0:.2f}")
    
    sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                   t_eval=np.linspace(0, 200, 1000))
    
    ax1.plot(sol.t, sol.y[0], 'b-', linewidth=2, label='Susceptible')
    ax1.plot(sol.t, sol.y[1], 'r-', linewidth=2, label='Infected')
    ax1.plot(sol.t, sol.y[2], 'g-', linewidth=2, label='Recovered')
    ax1.set_xlabel('Time (days)', fontsize=12)
    ax1.set_ylabel('Population', fontsize=12)
    ax1.set_title(f'SIR Model ($R_0$= {R_0:.1f})', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 不同 R0 的影响
    ax2 = axes[1]
    
    R0_values = [1.5, 2.5, 4.0]
    
    for R0_val in R0_values:
        beta = R0_val * gamma * N / S0
        sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                       t_eval=np.linspace(0, 200, 1000))
        ax2.plot(sol.t, sol.y[1], linewidth=2, label=f'$R_0$= {R0_val}')
    
    ax2.set_xlabel('Time (days)', fontsize=12)
    ax2.set_ylabel('Infected', fontsize=12)
    ax2.set_title('Effect of $R_0$ on Peak', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 相平面
    ax3 = axes[2]
    
    beta = 0.3
    sol = solve_ivp(sir, [0, 200], [S0, I0, R0], args=(beta, gamma),
                   t_eval=np.linspace(0, 200, 1000))
    
    ax3.plot(sol.y[0], sol.y[1], 'b-', linewidth=2)
    ax3.plot(S0, I0, 'go', markersize=10, label='Start')
    ax3.plot(sol.y[0][-1], sol.y[1][-1], 'ro', markersize=10, label='End')
    
    # 阈值线
    S_threshold = gamma * N / beta
    ax3.axvline(x=S_threshold, color='k', linestyle='--', label=f'Threshold S = {S_threshold:.0f}')
    
    ax3.set_xlabel('Susceptible', fontsize=12)
    ax3.set_ylabel('Infected', fontsize=12)
    ax3.set_title('SIR Phase Plane', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('sir_model.png', dpi=150, bbox_inches='tight')
    plt.show()

sir_model()

控制系统

反馈控制的基本原理

考虑一个被控系统：

其中是控制输入。

比例控制：（负反馈）

闭环系统：稳定性条件：

PID 控制器

u(t) = K_p e(t) + K_i _0^t e()d+ K_d $$ - ：比例增益（即时响应） - ：积分增益（消除稳态误差） - ：微分增益（预测趋势，抑制振荡）

def control_systems():
    """控制系统分析"""
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 开环 vs 闭环
    ax1 = axes[0, 0]
    
    a = 1  # 不稳定的开环系统
    b = 1
    
    def open_loop(t, x):
        return a * x
    
    def closed_loop(t, x, K):
        return (a - b*K) * x
    
    t_span = (0, 5)
    t_eval = np.linspace(0, 5, 500)
    x0 = [1]
    
    # 开环（不稳定）
    sol_open = solve_ivp(open_loop, t_span, x0, t_eval=t_eval)
    ax1.plot(sol_open.t, sol_open.y[0], 'r-', linewidth=2, label='Open loop (unstable)')
    
    # 闭环（不同 K）
    for K in [1, 2, 5]:
        sol_closed = solve_ivp(closed_loop, t_span, x0, args=(K,), t_eval=t_eval)
        ax1.plot(sol_closed.t, sol_closed.y[0], linewidth=2, label=f'Closed loop K={K}')
    
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('State x', fontsize=12)
    ax1.set_title('Feedback Stabilization', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim(-1, 5)
    
    # PID 控制器
    ax2 = axes[0, 1]
    
    def plant_with_pid(t, state, Kp, Ki, Kd, setpoint):
        x, integral_e = state
        e = setpoint - x
        # 简单的一阶系统： dx/dt = -x + u
        # PID 控制： u = Kp*e + Ki*integral + Kd*(-dx/dt) ≈ Kp*e + Ki*integral - Kd*v
        # 其中 v 是 x 的导数估计
        u = Kp * e + Ki * integral_e
        dxdt = -x + u
        return [dxdt, e]
    
    setpoint = 1.0
    t_span = (0, 20)
    t_eval = np.linspace(0, 20, 1000)
    
    # P 控制
    sol_p = solve_ivp(plant_with_pid, t_span, [0, 0], args=(2, 0, 0, setpoint), t_eval=t_eval)
    ax2.plot(sol_p.t, sol_p.y[0], 'b-', linewidth=2, label='P only')
    
    # PI 控制
    sol_pi = solve_ivp(plant_with_pid, t_span, [0, 0], args=(2, 1, 0, setpoint), t_eval=t_eval)
    ax2.plot(sol_pi.t, sol_pi.y[0], 'g-', linewidth=2, label='PI')
    
    ax2.axhline(y=setpoint, color='k', linestyle='--', label='Setpoint')
    ax2.set_xlabel('Time', fontsize=12)
    ax2.set_ylabel('Output', fontsize=12)
    ax2.set_title('PID Control Response', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 二阶系统的极点配置
    ax3 = axes[1, 0]
    
    # 二阶系统： x'' + 2 ζω₀ x' + ω₀² x = ω₀² r
    def second_order(t, state, zeta, omega0, r):
        x, v = state
        return [v, omega0**2 * (r - x) - 2*zeta*omega0*v]
    
    omega0 = 2
    r = 1  # 参考输入
    
    for zeta in [0.3, 0.7, 1.0, 2.0]:
        sol = solve_ivp(second_order, [0, 10], [0, 0], args=(zeta, omega0, r),
                       t_eval=np.linspace(0, 10, 500))
        ax3.plot(sol.t, sol.y[0], linewidth=2, label=f'ζ = {zeta}')
    
    ax3.axhline(y=r, color='k', linestyle='--')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Output', fontsize=12)
    ax3.set_title('Second-Order System Response', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 根轨迹概念
    ax4 = axes[1, 1]
    
    # 画不同增益下的极点
    K_values = np.linspace(0, 10, 50)
    
    # 简单系统： 1/(s+1)(s+2)的闭环极点
    # 特征方程： s ² + 3s + 2 + K = 0
    # s = (-3 ± sqrt(9 - 4(2+K)))/2
    
    for K in K_values:
        discriminant = 9 - 4*(2+K)
        if discriminant >= 0:
            s1 = (-3 + np.sqrt(discriminant))/2
            s2 = (-3 - np.sqrt(discriminant))/2
            ax4.plot(s1, 0, 'b.', markersize=3)
            ax4.plot(s2, 0, 'b.', markersize=3)
        else:
            real_part = -1.5
            imag_part = np.sqrt(-discriminant)/2
            ax4.plot(real_part, imag_part, 'b.', markersize=3)
            ax4.plot(real_part, -imag_part, 'b.', markersize=3)
    
    # 开环极点
    ax4.plot(-1, 0, 'rx', markersize=15, markeredgewidth=3, label='Open-loop poles')
    ax4.plot(-2, 0, 'rx', markersize=15, markeredgewidth=3)
    
    ax4.axvline(x=0, color='k', linewidth=0.5)
    ax4.axhline(y=0, color='k', linewidth=0.5)
    ax4.set_xlabel('Real', fontsize=12)
    ax4.set_ylabel('Imaginary', fontsize=12)
    ax4.set_title('Root Locus (Conceptual)', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim(-4, 1)
    ax4.set_ylim(-3, 3)
    
    plt.tight_layout()
    plt.savefig('control_systems.png', dpi=150, bbox_inches='tight')
    plt.show()

control_systems()

流体力学中的 ODE

Navier-Stokes 方程的简化

完整的 Navier-Stokes 方程是偏微分方程，但在某些对称情况下可以简化为 ODE 。

Poiseuille 流（管道流）：稳态层流满足

解为抛物线速度分布：

u(r) = (-)(R^2 - r^2) $$

边界层方程

Blasius 方程描述平板边界层：

f''' + ff'' = 0 $$

边界条件：$ f(0) = f'(0) = 0 $，$ f'() = 1$ 这是一个三阶非线性 BVP（边值问题）。

def fluid_mechanics():
    """流体力学中的 ODE"""
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    # Poiseuille 流
    ax1 = axes[0]
    
    R = 1  # 管道半径
    r = np.linspace(0, R, 100)
    
    # 归一化的速度分布
    u = 1 - (r/R)**2
    
    ax1.plot(u, r, 'b-', linewidth=2)
    ax1.plot(u, -r, 'b-', linewidth=2)
    ax1.fill_betweenx(r, 0, u, alpha=0.3)
    ax1.fill_betweenx(-r, 0, u, alpha=0.3)
    
    ax1.axvline(x=0, color='k', linewidth=2)
    ax1.axhline(y=R, color='k', linewidth=2)
    ax1.axhline(y=-R, color='k', linewidth=2)
    
    ax1.set_xlabel('Velocity u/u_max', fontsize=12)
    ax1.set_ylabel('Radial position r', fontsize=12)
    ax1.set_title('Poiseuille Flow Profile', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    
    # Blasius 边界层
    ax2 = axes[1]
    
    from scipy.integrate import solve_bvp
    
    def blasius_ode(eta, y):
        f, f_prime, f_double_prime = y
        return [f_prime, f_double_prime, -0.5 * f * f_double_prime]
    
    def blasius_bc(ya, yb):
        return [ya[0], ya[1], yb[1] - 1]
    
    eta = np.linspace(0, 10, 100)
    y_init = np.zeros((3, eta.size))
    y_init[1] = eta / 10  # 初始猜测
    
    sol = solve_bvp(blasius_ode, blasius_bc, eta, y_init)
    
    ax2.plot(sol.y[1], sol.x, 'b-', linewidth=2, label="$ f'(\\eta)$= u/U")
    ax2.axhline(y=5, color='k', linestyle='--', alpha=0.5, label='δ₉₉')
    ax2.set_xlabel("$ f'$(normalized velocity)", fontsize=12)
    ax2.set_ylabel("$\\eta$(similarity variable)", fontsize=12)
    ax2.set_title('Blasius Boundary Layer', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 1.2)
    
    # 边界层厚度随 x 变化
    ax3 = axes[2]
    
    x = np.linspace(0.1, 10, 100)
    Re_x = 1000 * x  # 假设某个 Re 数
    delta = 5 * x / np.sqrt(Re_x)  # Blasius 结果
    
    ax3.plot(x, delta, 'g-', linewidth=2)
    ax3.set_xlabel('x (distance from leading edge)', fontsize=12)
    ax3.set_ylabel('δ (boundary layer thickness)', fontsize=12)
    ax3.set_title('Boundary Layer Growth', fontsize=14, fontweight='bold')
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fluid_mechanics.png', dpi=150, bbox_inches='tight')
    plt.show()

fluid_mechanics()

陀螺仪与刚体动力学

欧拉方程

刚体绕固定点转动的欧拉方程：

I_1_1 = (I_2 - I_3)_2_3 + M_1 $

$$ I_2_2 = (I_3 - I_1)_3_1 + M_2 $

这是一组非线性耦合 ODE 。

无力矩对称陀螺

当且无外力矩时：

其中。

解为圆周运动：这就是进动（ precession）。

def gyroscope_dynamics():
    """陀螺仪动力学"""
    
    def euler_equations(t, state, I1, I2, I3):
        omega1, omega2, omega3 = state
        domega1 = (I2 - I3) * omega2 * omega3 / I1
        domega2 = (I3 - I1) * omega3 * omega1 / I2
        domega3 = (I1 - I2) * omega1 * omega2 / I3
        return [domega1, domega2, domega3]
    
    fig = plt.figure(figsize=(14, 10))
    
    # 对称陀螺 (I1 = I2)
    ax1 = fig.add_subplot(2, 2, 1)
    
    I1, I2, I3 = 1.0, 1.0, 0.5  # 扁平陀螺
    omega0 = [0.1, 0, 10]  # 主要绕 z 轴旋转
    
    sol = solve_ivp(euler_equations, [0, 10], omega0, args=(I1, I2, I3),
                   t_eval=np.linspace(0, 10, 1000))
    
    ax1.plot(sol.t, sol.y[0], 'r-', linewidth=2, label='$\\omega_1$')
    ax1.plot(sol.t, sol.y[1], 'g-', linewidth=2, label='$\\omega_2$')
    ax1.plot(sol.t, sol.y[2], 'b-', linewidth=2, label='$\\omega_3$')
    ax1.set_xlabel('Time', fontsize=12)
    ax1.set_ylabel('Angular velocity', fontsize=12)
    ax1.set_title('Symmetric Top (I ₁ = I ₂)', fontsize=14, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 三维轨迹
    ax2 = fig.add_subplot(2, 2, 2, projection='3d')
    ax2.plot(sol.y[0], sol.y[1], sol.y[2], 'b-', linewidth=1)
    ax2.set_xlabel('$\\omega_1$')
    ax2.set_ylabel('$\\omega_2$')
    ax2.set_zlabel('$\\omega_3$')
    ax2.set_title('Angular Velocity Trajectory', fontsize=14, fontweight='bold')
    
    # 非对称情况（ Euler top）
    ax3 = fig.add_subplot(2, 2, 3)
    
    I1, I2, I3 = 1.0, 2.0, 3.0
    omega0 = [1, 0.1, 0.1]
    
    sol = solve_ivp(euler_equations, [0, 50], omega0, args=(I1, I2, I3),
                   t_eval=np.linspace(0, 50, 2000))
    
    ax3.plot(sol.t, sol.y[0], 'r-', linewidth=1, label='$\\omega_1$')
    ax3.plot(sol.t, sol.y[1], 'g-', linewidth=1, label='$\\omega_2$')
    ax3.plot(sol.t, sol.y[2], 'b-', linewidth=1, label='$\\omega_3$')
    ax3.set_xlabel('Time', fontsize=12)
    ax3.set_ylabel('Angular velocity', fontsize=12)
    ax3.set_title('Asymmetric Top (I ₁ ≠ I ₂ ≠ I ₃)', fontsize=14, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 能量和角动量守恒
    ax4 = fig.add_subplot(2, 2, 4)
    
    # 计算能量和角动量
    E = 0.5 * (I1*sol.y[0]**2 + I2*sol.y[1]**2 + I3*sol.y[2]**2)
    L_sq = (I1*sol.y[0])**2 + (I2*sol.y[1])**2 + (I3*sol.y[2])**2
    
    ax4.plot(sol.t, E/E[0], 'b-', linewidth=2, label='Energy (normalized)')
    ax4.plot(sol.t, np.sqrt(L_sq)/np.sqrt(L_sq[0]), 'r--', linewidth=2, label='|L| (normalized)')
    ax4.set_xlabel('Time', fontsize=12)
    ax4.set_ylabel('Conserved quantity', fontsize=12)
    ax4.set_title('Conservation Laws', fontsize=14, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    ax4.set_ylim(0.99, 1.01)
    
    plt.tight_layout()
    plt.savefig('gyroscope_dynamics.png', dpi=150, bbox_inches='tight')
    plt.show()

gyroscope_dynamics()

总结

本章我们深入探索了常微分方程在物理和工程中的广泛应用：

力学系统

质点运动：从简单的自由落体到复杂的行星运动
振动系统：单自由度和多自由度振动，线性和非线性
刚体动力学：欧拉方程和陀螺运动

电磁学

电路分析： RLC 电路、耦合电路、非线性元件
频率响应：传递函数和 Bode 图

热学与化学

热传导：牛顿冷却定律及其应用
反应动力学：一级、连续、可逆反应

生物与生态

种群动力学：捕食者-猎物模型、竞争排斥
传染病模型： SIR 模型和基本再生数

控制系统

反馈控制：稳定性和极点配置
PID 控制：工业应用中最常见的控制策略

流体力学

简化模型： Poiseuille 流、边界层理论

练习题

基础题

一个 kg 的物体从高空落下，空气阻力为。求终端速度和到达终端速度 99%所需的时间。
求解带阻尼的简谐振子，初始条件$ x(0) = 1 $，$ (0) = 0$。
RLC 电路中H，μ F，Ω。判断电路是欠阻尼还是过阻尼，并求自然频率。
一级反应的半衰期为 2 小时。求反应物浓度降至初始值 10%所需时间。
在 SIR 模型中，如果/天，/天，计算并解释其物理意义。

进阶题

双摆系统：推导双摆的运动方程，并数值求解。观察对初始条件的敏感性。
Van der Pol 振子：

分析不同值时的极限环行为。

化学振荡（ Belousov-Zhabotinsky 反应的简化模型）：

探索参数空间中的振荡行为。

航天器姿态控制：设计一个简单的反馈控制律，使卫星的姿态角稳定在零位置。
流行病控制：在 SIR 模型中引入疫苗接种率$ v $，分析如何选择$ v$ 以防止疫情爆发。

编程题

实现一个通用的二阶系统模拟器，可以分析不同阻尼比下的阶跃响应、频率响应和相平面轨迹。
编写程序求解三体问题（限制性），并可视化轨道。
实现 PID 控制器，用于温度控制系统的模拟。包括自整定功能。
模拟多自由度建筑结构在地震激励下的响应，分析不同层间刚度的影响。

思考题

为什么陀螺仪能保持稳定的指向？从角动量守恒的角度解释。
在控制系统中，为什么纯积分控制会导致不稳定或振荡？
生态系统中，为什么 Lotka-Volterra 模型预测的周期振荡在现实中很少观察到？
从物理直觉解释为什么临界阻尼是"最优"的——既不振荡，又最快到达平衡。

参考资料

Strogatz, S. H. (2018). Nonlinear Dynamics and Chaos. CRC Press.
Ogata, K. (2010). Modern Control Engineering. Prentice Hall.
Murray, J. D. (2002). Mathematical Biology. Springer.
Thornton, S. T., & Marion, J. B. (2004). Classical Dynamics of Particles and Systems. Brooks/Cole.
Dorf, R. C., & Bishop, R. H. (2017). Modern Control Systems. Pearson.
Batchelor, G. K. (2000). An Introduction to Fluid Dynamics. Cambridge University Press.

本文是《常微分方程的世界》系列的第 17 章。

经典力学中的微分方程

牛顿运动定律的数学表述

自由落体与空气阻力

行星运动与开普勒问题

振动系统的深入分析

多自由度振动系统

三自由度弹簧系统

非线性振动： Duffing 振子

电路与电磁学

一般 RLC 网络

耦合电路：变压器

非线性电路：二极管和晶体管

热传导与扩散

牛顿冷却定律

多区域热传导

化学反应动力学

一级反应

连续反应

可逆反应与化学平衡

生物系统与种群动力学

Lotka-Volterra 捕食者-猎物模型

竞争排斥原理

传染病模型： SIR

控制系统

反馈控制的基本原理

PID 控制器

流体力学中的 ODE

Navier-Stokes 方程的简化

边界层方程

陀螺仪与刚体动力学

欧拉方程

无力矩对称陀螺

总结

力学系统

电磁学

热学与化学

生物与生态

控制系统

流体力学

练习题

基础题

进阶题

编程题

思考题

参考资料

Batchelor, G. K. (2000). An Introduction to Fluid Dynamics. Cambridge University Press.