常微分方程（二）一阶微分方程

一阶微分方程是微分方程中最基本也是最常用的类型。从银行复利计算到药物代谢、从化学反应到电路分析，一阶 ODE 无处不在。掌握这些方程的求解技巧，是理解更复杂系统的基础。

可分离方程：最自然的形式

什么是可分离方程？

考虑一阶 ODE 的一般形式：

如果 可以写成 的函数和 的函数的乘积：

这就是可分离方程。

核心思想：把所有含 的项移到等式左边，含 的项移到右边，然后两边分别积分。

经典例子：指数增长与衰减

例 1：求解（ 为常数）

解：这是可分离方程，，。

分离变量：

两边积分：

令（包括 的情况），得到通解：

物理意义： -：指数增长（人口增长、细菌繁殖、复利） -：指数衰减（放射性衰变、药物代谢、咖啡冷却）

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 不同 k 值的指数函数
x = np.linspace(0, 5, 200)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 指数增长
for k in [0.5, 1.0, 1.5, 2.0]:
    axes[0].plot(x, np.exp(k * x), label=f'k = {k}', linewidth=2)
axes[0].set_xlabel('x', fontsize=12)
axes[0].set_ylabel('y', fontsize=12)
axes[0].set_title('Exponential Growth (k > 0)', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(0, 30)

# 指数衰减
for k in [-0.5, -1.0, -1.5, -2.0]:
    axes[1].plot(x, np.exp(k * x), label=f'k = {k}', linewidth=2)
axes[1].set_xlabel('x', fontsize=12)
axes[1].set_ylabel('y', fontsize=12)
axes[1].set_title('Exponential Decay (k < 0)', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('第 2 章-一阶微分方程/fig1_exponential.png', dpi=150, bbox_inches='tight')

![指数增长与衰减对比](./常微分方程（二）一阶微分方程/fig1_exponential.png)
plt.close()

银行复利问题

实际场景：你在银行存入 元，年利率为（连续复利），求 年后的本金。

建模：设 为 时刻的本金，变化率等于利息：

初始条件： 求解：

实例计算：存入 10000 元，年利率 5%

def compound_interest():
    P0 = 10000  # 初始本金
    r = 0.05    # 年利率 5%
    t = np.linspace(0, 30, 300)
    
    # 连续复利
    P_continuous = P0 * np.exp(r * t)
    
    # 年复利对比
    P_yearly = P0 * (1 + r) ** t
    
    # 月复利
    P_monthly = P0 * (1 + r/12) ** (12 * t)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(t, P_continuous, 'b-', linewidth=2.5, label='Continuous Compounding')
    ax.plot(t, P_monthly, 'g--', linewidth=2, label='Monthly Compounding')
    ax.plot(t, P_yearly, 'r:', linewidth=2, label='Yearly Compounding')
    
    ax.set_xlabel('Years', fontsize=12)
    ax.set_ylabel('Principal ($)', fontsize=12)
    ax.set_title(f'Compound Interest Comparison (r = {r*100}%)', fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    
    # 标注 10 年和 30 年的值
    for t_mark in [10, 30]:
        idx = np.abs(t - t_mark).argmin()
        ax.annotate(f'${P_continuous[idx]:.0f}', 
                   xy=(t_mark, P_continuous[idx]),
                   xytext=(t_mark+2, P_continuous[idx]+2000),
                   fontsize=10, arrowprops=dict(arrowstyle='->', color='blue'))
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig2_compound_interest.png', dpi=150, bbox_inches='tight')

![连续复利与离散复利对比](./常微分方程（二）一阶微分方程/fig2_compound_interest.png)
    plt.close()

compound_interest()

关键洞察：

时间	年复利	连续复利	差额
10 年	$16,288.95	$16,487.21	$198.26
30 年	$43,219.42	$44,816.89	$1,597.47

连续复利看似只多一点点利息，但长期来看差异显著。这就是 的神奇之处！

药物代谢：半衰期概念

背景：你吃了一片药，身体以固定比例代谢药物。多久药效减半？

模型：设 为血液中药物浓度

解： 半衰期：浓度降为一半的时间

实例：布洛芬的半衰期约为 2 小时，如果初始服用 400mg

def drug_metabolism():
    C0 = 400  # 初始浓度 mg
    t_half = 2  # 半衰期 hours
    k = np.log(2) / t_half
    
    t = np.linspace(0, 12, 200)
    C = C0 * np.exp(-k * t)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(t, C, 'b-', linewidth=2.5, label='Drug Concentration')
    ax.axhline(y=200, color='r', linestyle='--', label='50% (1 half-life)')
    ax.axhline(y=100, color='orange', linestyle='--', label='25% (2 half-lives)')
    ax.axhline(y=50, color='green', linestyle='--', label='12.5% (3 half-lives)')
    
    # 标注半衰期
    for i in range(1, 4):
        t_i = i * t_half
        C_i = C0 / (2**i)
        ax.plot(t_i, C_i, 'ko', markersize=8)
        ax.text(t_i + 0.2, C_i + 15, f't = {t_i}h\nC = {C_i}mg', fontsize=9)
    
    ax.set_xlabel('Time (hours)', fontsize=12)
    ax.set_ylabel('Drug Concentration (mg)', fontsize=12)
    ax.set_title('Drug Metabolism: Ibuprofen (Half-life = 2h)', fontsize=14, fontweight='bold')
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 12)
    ax.set_ylim(0, 450)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig3_drug_metabolism.png', dpi=150, bbox_inches='tight')

![药物代谢的半衰期曲线](./常微分方程（二）一阶微分方程/fig3_drug_metabolism.png)
    plt.close()

drug_metabolism()

临床意义： - 通常 3-5 个半衰期后，药物几乎完全代谢 - 用药间隔通常基于半衰期设计 - 布洛芬：每 4-6 小时服用一次（约 2-3 个半衰期）

Logistic 方程：资源限制下的增长

马尔萨斯模型的问题：指数增长 会导致人口无限增长，这显然不现实。

改进： Pierre Fran ç ois Verhulst（ 1838 年）提出 Logistic 模型 -：内禀增长率 -：环境容纳量（ carrying capacity）

直觉： - 当：（指数增长） - 当：（增长停止） - 当：（人口下降）

求解：这是可分离方程

用部分分数分解：

积分得到：

整理得到Logistic 函数：

def logistic_growth():
    K = 1000  # 环境容纳量
    r = 0.3   # 增长率
    
    t = np.linspace(0, 30, 300)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 不同初始值的 Logistic 曲线
    P0_values = [10, 50, 100, 500, 1200]
    colors = plt.cm.viridis(np.linspace(0, 1, len(P0_values)))
    
    for P0, color in zip(P0_values, colors):
        P = K / (1 + (K/P0 - 1) * np.exp(-r * t))
        axes[0].plot(t, P, color=color, linewidth=2, label=f'P ₀ = {P0}')
    
    axes[0].axhline(y=K, color='red', linestyle='--', linewidth=2, label=f'K = {K}')
    axes[0].set_xlabel('Time', fontsize=12)
    axes[0].set_ylabel('Population', fontsize=12)
    axes[0].set_title('Logistic Growth: Different Initial Conditions', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 增长率 vs 人口
    P_range = np.linspace(0, 1.5 * K, 200)
    dP_dt = r * P_range * (1 - P_range / K)
    
    axes[1].plot(P_range, dP_dt, 'b-', linewidth=2.5)
    axes[1].axvline(x=K/2, color='green', linestyle='--', label='Maximum growth at P = K/2')
    axes[1].axvline(x=K, color='red', linestyle='--', label='Carrying capacity K')
    axes[1].axhline(y=0, color='black', linewidth=0.5)
    axes[1].fill_between(P_range[P_range <= K], dP_dt[P_range <= K], alpha=0.3, color='blue')
    axes[1].fill_between(P_range[P_range > K], dP_dt[P_range > K], alpha=0.3, color='red')
    
    axes[1].set_xlabel('Population P', fontsize=12)
    axes[1].set_ylabel('Growth Rate dP/dt', fontsize=12)
    axes[1].set_title('Growth Rate vs Population', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig4_logistic.png', dpi=150, bbox_inches='tight')

![Logistic增长模型的S形曲线](./常微分方程（二）一阶微分方程/fig4_logistic.png)
    plt.close()

logistic_growth()

S 形曲线的特征：

1. 拐点：在 处，增长最快
2. 渐近线：（永远不会超过太多）
3. 对称性：绕拐点中心对称

现实应用： - 细菌在培养皿中的生长 - 新技术的采用曲线 - 传染病的传播（早期） - 市场饱和模型

一阶线性方程：积分因子法

标准形式

一阶线性 ODE 的标准形式：

特点： 和都是一次的，且系数可以是 的函数。

积分因子法的动机

考虑乘积的导数：

如果能找到一个函数，使得：

那么方程变得容易积分！

条件：比较系数，需要

这本身是可分离方程：

这就是积分因子！

通解公式

将 乘到原方程两边：

积分：

通解：

例题详解

例 2：求解 识别：， 计算积分因子：

乘以积分因子：

左边是：

积分： （用换元）

通解：

def linear_ode_example():
    x = np.linspace(-3, 3, 300)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # 不同 C 值的解曲线
    C_values = [-2, -1, -0.5, 0, 0.5, 1, 2]
    colors = plt.cm.coolwarm(np.linspace(0, 1, len(C_values)))
    
    for C, color in zip(C_values, colors):
        y = 0.5 + C * np.exp(-x**2)
        ax.plot(x, y, color=color, linewidth=2, label=f'C = {C}')
    
    # 平衡解
    ax.axhline(y=0.5, color='black', linestyle='--', linewidth=2, label='Equilibrium y = 1/2')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title("Solution Family: y' + 2xy = x", fontsize=14, fontweight='bold')
    ax.legend(loc='upper right')
    ax.grid(True, alpha=0.3)
    ax.set_ylim(-2, 3)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig5_linear_ode.png', dpi=150, bbox_inches='tight')

![一阶线性方程的解族](./常微分方程（二）一阶微分方程/fig5_linear_ode.png)
    plt.close()

linear_ode_example()

观察：所有解曲线都趋向于平衡解！这是因为 当。

RC 电路：一阶线性 ODE 的经典应用

场景：电阻 和电容 串联，接入电压源。求电容上的电压。

电路方程（基尔霍夫电压定律）：

由于：

整理为标准形式：

定义时间常数：

情况 1：恒定电压源

积分因子： 通解：

物理意义： - 电容从初始电压 指数趋近于 - 时间常数 决定了充放电速度 - 时，完成约 63.2%的充电 - 时，几乎完全充电（ 99.3%）

def rc_circuit():
    R = 1000  # 1k Ω
    C = 0.001  # 1mF
    tau = R * C  # 1s
    V0 = 5  # 电源电压 5V
    
    t = np.linspace(0, 5 * tau, 300)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 充电过程（初始 0V）
    Vc_charge = V0 * (1 - np.exp(-t / tau))
    i_charge = (V0 / R) * np.exp(-t / tau)
    
    axes[0].plot(t, Vc_charge, 'b-', linewidth=2.5, label='Capacitor Voltage$V_C$')
    axes[0].plot(t, V0 - Vc_charge, 'r--', linewidth=2, label='Resistor Voltage$V_R$')
    axes[0].axhline(y=V0, color='green', linestyle=':', label=f'Source Voltage = {V0}V')
    axes[0].axvline(x=tau, color='gray', linestyle='--', alpha=0.5)
    axes[0].text(tau, 0.2, f'τ = {tau}s', fontsize=10)
    
    # 标注 63.2%点
    axes[0].plot(tau, V0 * 0.632, 'ko', markersize=8)
    axes[0].annotate('63.2%', xy=(tau, V0*0.632), xytext=(tau+0.5, V0*0.5),
                    fontsize=10, arrowprops=dict(arrowstyle='->', color='black'))
    
    axes[0].set_xlabel('Time (s)', fontsize=12)
    axes[0].set_ylabel('Voltage (V)', fontsize=12)
    axes[0].set_title('RC Circuit: Charging', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 放电过程（初始 5V，断开电源）
    Vc_discharge = V0 * np.exp(-t / tau)
    
    axes[1].plot(t, Vc_discharge, 'b-', linewidth=2.5, label='Capacitor Voltage$V_C$')
    axes[1].axhline(y=0, color='green', linestyle=':', label='Ground')
    axes[1].axvline(x=tau, color='gray', linestyle='--', alpha=0.5)
    
    # 标注 36.8%点
    axes[1].plot(tau, V0 * 0.368, 'ko', markersize=8)
    axes[1].annotate('36.8%', xy=(tau, V0*0.368), xytext=(tau+0.5, V0*0.5),
                    fontsize=10, arrowprops=dict(arrowstyle='->', color='black'))
    
    axes[1].set_xlabel('Time (s)', fontsize=12)
    axes[1].set_ylabel('Voltage (V)', fontsize=12)
    axes[1].set_title('RC Circuit: Discharging', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig6_rc_circuit.png', dpi=150, bbox_inches='tight')

![RC电路的充放电过程](./常微分方程（二）一阶微分方程/fig6_rc_circuit.png)
    plt.close()

rc_circuit()

恰当方程与积分因子

什么是恰当方程？

考虑形如

或等价地写成微分形式：

如果存在函数 使得：

即，，则方程是恰当的。

恰当性检验：

这是因为（混合偏导相等）。

求解恰当方程

如果方程恰当，通解就是。

步骤：

1. 验证
2. 从 对 积分，得
3. 用 确定
4. 通解为

例 3：求解 验证恰当性： -， -， 相等！方程恰当。

求：

利用：

所以，（常数可忽略）。

通解：

使方程变恰当：积分因子

当方程不恰当时（），可以尝试找一个积分因子 使得：

变成恰当方程。

特殊情况 1：（只依赖于）

需要 只依赖于。

此时：

特殊情况 2：（只依赖于）

需要 只依赖于。

此时：

伯努利方程

形式与化简

伯努利方程：

其中（否则就是线性方程）。

技巧：做变量替换 则 原方程两边除以：

即：

这是关于 的线性方程！

例 4：求解 这是 的伯努利方程。

令，则 原方程除以：

这是线性方程！积分因子

用分部积分：

回代：

混合问题：槽中的盐水

问题设定

一个容量为 1000 升的水槽，初始装有纯水。现在以 5 升/分钟的速率注入含盐 2 克/升的盐水，同时以 5 升/分钟的速率排出混合液。求任意时刻槽中盐的含量。

建模

设 为 时刻槽中盐的克数。

盐的流入速率： 克/分钟

盐的流出速率： 克/分钟

（因为流出的浓度 = 当前槽中浓度 = 克/升）

微分方程：

初始条件：（纯水）

求解

这是一阶线性方程。整理为标准形式：

积分因子：

用：，所以 解：

def salt_tank():
    t = np.linspace(0, 1500, 300)
    Q = 2000 * (1 - np.exp(-t / 200))
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # 盐含量随时间变化
    axes[0].plot(t, Q, 'b-', linewidth=2.5, label='Salt Amount Q(t)')
    axes[0].axhline(y=2000, color='red', linestyle='--', linewidth=2, label='Equilibrium = 2000g')
    
    # 标注关键时刻
    for t_mark in [200, 400, 600, 1000]:
        Q_mark = 2000 * (1 - np.exp(-t_mark / 200))
        axes[0].plot(t_mark, Q_mark, 'ko', markersize=6)
        axes[0].annotate(f't={t_mark}min\nQ={Q_mark:.0f}g', 
                        xy=(t_mark, Q_mark), xytext=(t_mark+50, Q_mark-200),
                        fontsize=9)
    
    axes[0].set_xlabel('Time (minutes)', fontsize=12)
    axes[0].set_ylabel('Salt Amount (grams)', fontsize=12)
    axes[0].set_title('Salt Tank Problem', fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 浓度变化
    C = Q / 1000  # 克/升
    axes[1].plot(t, C, 'g-', linewidth=2.5, label='Concentration C(t)')
    axes[1].axhline(y=2, color='red', linestyle='--', linewidth=2, label='Inlet Concentration = 2 g/L')
    
    axes[1].set_xlabel('Time (minutes)', fontsize=12)
    axes[1].set_ylabel('Concentration (g/L)', fontsize=12)
    axes[1].set_title('Salt Concentration Over Time', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig7_salt_tank.png', dpi=150, bbox_inches='tight')

![盐水混合问题的浓度变化](./常微分方程（二）一阶微分方程/fig7_salt_tank.png)
    plt.close()

salt_tank()

关键洞察： - 槽中盐量趋向于 2000 克（浓度趋向 2 克/升，等于进入浓度） - 时间常数 分钟 - 约 1000 分钟（约 17 小时）后，基本达到平衡

变体：水槽容量变化

如果流入速率 流出速率，水槽水量会变化！

例：流入 6 升/分钟（ 2 克/升），流出 4 升/分钟，初始 500 升纯水。

水量： 升

微分方程：

这需要更复杂的积分因子，但原理相同。

存在唯一性定理

Picard-Lindel ö f 定理

定理：考虑初值问题

如果 在包含 的矩形区域 上满足： 1. 关于 连续 2. 存在且关于 连续（ Lipschitz 条件）

则在附近存在唯一解。

为什么重要？

正面：告诉我们什么时候可以放心求解

反面：当条件不满足时，可能出现： - 无解 - 无穷多解 - 解只存在于有限区间

例 5：， 检验： 在 处无界！

这个初值问题有无穷多解： -（平凡解） - 对任意（非平凡解） - - 分段函数...

def non_unique_solution():
    x = np.linspace(-1, 3, 300)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    # y = 0（平凡解）
    ax.axhline(y=0, color='blue', linewidth=2.5, label='y = 0 (trivial)')
    
    # y = x^3
    ax.plot(x, x**3, 'r-', linewidth=2, label='y = x ³')
    
    # y = (x-1)^3 for x >= 1, 0 otherwise
    y1 = np.where(x >= 1, (x-1)**3, 0)
    ax.plot(x, y1, 'g--', linewidth=2, label='y = (x-1)³ for x ≥ 1')
    
    # y = (x-2)^3 for x >= 2, 0 otherwise
    y2 = np.where(x >= 2, (x-2)**3, 0)
    ax.plot(x, y2, 'm:', linewidth=2, label='y = (x-2)³ for x ≥ 2')
    
    ax.plot(0, 0, 'ko', markersize=10, label='Initial point (0, 0)')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title("Non-uniqueness: dy/dx = 3y^(2/3), y(0) = 0", fontsize=14, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_xlim(-1, 3)
    ax.set_ylim(-2, 10)
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig8_non_unique.png', dpi=150, bbox_inches='tight')
    plt.close()

non_unique_solution()

数值方法入门

欧拉方法（ Euler's Method）

思想：用切线近似曲线

其中 是步长。

几何意义：从 出发，沿斜率 走一小步。

def euler_method_demo():
    # ODE: dy/dx = y, y(0) = 1
    # 精确解: y = e^x
    
    def f(x, y):
        return y
    
    def euler(f, x0, y0, h, n_steps):
        x = [x0]
        y = [y0]
        for _ in range(n_steps):
            y_new = y[-1] + h * f(x[-1], y[-1])
            x_new = x[-1] + h
            x.append(x_new)
            y.append(y_new)
        return np.array(x), np.array(y)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    x_exact = np.linspace(0, 2, 200)
    y_exact = np.exp(x_exact)
    
    # 不同步长对比
    step_sizes = [0.5, 0.2, 0.1, 0.05]
    colors = ['red', 'orange', 'green', 'purple']
    
    for h, color in zip(step_sizes, colors):
        n_steps = int(2 / h)
        x_euler, y_euler = euler(f, 0, 1, h, n_steps)
        axes[0].plot(x_euler, y_euler, 'o-', color=color, markersize=4, 
                    linewidth=1.5, label=f'h = {h}')
    
    axes[0].plot(x_exact, y_exact, 'b-', linewidth=2.5, label='Exact:$y = e^x$')
    axes[0].set_xlabel('x', fontsize=12)
    axes[0].set_ylabel('y', fontsize=12)
    axes[0].set_title("Euler's Method: Different Step Sizes", fontsize=14, fontweight='bold')
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    # 误差分析
    step_sizes_log = np.logspace(-2, 0, 20)
    errors = []
    
    for h in step_sizes_log:
        n_steps = int(2 / h)
        x_euler, y_euler = euler(f, 0, 1, h, n_steps)
        error = np.abs(y_euler[-1] - np.exp(2))
        errors.append(error)
    
    axes[1].loglog(step_sizes_log, errors, 'bo-', linewidth=2, markersize=6, label='Error')
    axes[1].loglog(step_sizes_log, step_sizes_log, 'r--', linewidth=2, label='O(h) reference')
    axes[1].set_xlabel('Step Size h', fontsize=12)
    axes[1].set_ylabel('Error at x = 2', fontsize=12)
    axes[1].set_title('Error vs Step Size (Log-Log)', fontsize=14, fontweight='bold')
    axes[1].legend()
    axes[1].grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig9_euler.png', dpi=150, bbox_inches='tight')
    plt.close()

euler_method_demo()

误差分析： - 局部截断误差：（每步） - 全局误差：（累积）

欧拉方法是一阶方法。

改进的欧拉方法（ Heun's Method）

思想：用斜率的平均值

这是二阶方法，全局误差。

四阶 Runge-Kutta 方法（ RK4）

最常用的数值方法：

全局误差：

def rk4_comparison():
    def f(x, y):
        return y
    
    def euler(f, x0, y0, h, n_steps):
        x, y = [x0], [y0]
        for _ in range(n_steps):
            y.append(y[-1] + h * f(x[-1], y[-1]))
            x.append(x[-1] + h)
        return np.array(x), np.array(y)
    
    def rk4(f, x0, y0, h, n_steps):
        x, y = [x0], [y0]
        for _ in range(n_steps):
            k1 = f(x[-1], y[-1])
            k2 = f(x[-1] + h/2, y[-1] + h*k1/2)
            k3 = f(x[-1] + h/2, y[-1] + h*k2/2)
            k4 = f(x[-1] + h, y[-1] + h*k3)
            y.append(y[-1] + h/6 * (k1 + 2*k2 + 2*k3 + k4))
            x.append(x[-1] + h)
        return np.array(x), np.array(y)
    
    h = 0.2
    n_steps = 10
    
    x_euler, y_euler = euler(f, 0, 1, h, n_steps)
    x_rk4, y_rk4 = rk4(f, 0, 1, h, n_steps)
    x_exact = np.linspace(0, 2, 200)
    y_exact = np.exp(x_exact)
    
    fig, ax = plt.subplots(figsize=(10, 6))
    
    ax.plot(x_exact, y_exact, 'b-', linewidth=2.5, label='Exact')
    ax.plot(x_euler, y_euler, 'r--o', linewidth=2, markersize=6, label='Euler')
    ax.plot(x_rk4, y_rk4, 'g-s', linewidth=2, markersize=6, label='RK4')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(f'Euler vs RK4 (h = {h})', fontsize=14, fontweight='bold')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 显示误差
    euler_error = np.abs(y_euler[-1] - np.exp(2))
    rk4_error = np.abs(y_rk4[-1] - np.exp(2))
    ax.text(0.5, 5, f'Euler error: {euler_error:.4f}\nRK4 error: {rk4_error:.6f}', 
           fontsize=11, bbox=dict(boxstyle='round', facecolor='wheat'))
    
    plt.tight_layout()
    plt.savefig('第 2 章-一阶微分方程/fig10_rk4.png', dpi=150, bbox_inches='tight')
    plt.close()

rk4_comparison()

Python 实战：综合应用

使用 SciPy 求解 ODE

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

# 定义 ODE：冷却定律
def cooling(t, T, k, T_env):
    return -k * (T - T_env)

# 参数
k = 0.1
T_env = 20
T0 = 90

# 求解
sol = solve_ivp(
    fun=lambda t, T: cooling(t, T, k, T_env),
    t_span=[0, 60],
    y0=[T0],
    method='RK45',
    dense_output=True
)

# 绘图
t = np.linspace(0, 60, 300)
T = sol.sol(t)[0]

plt.figure(figsize=(10, 6))
plt.plot(t, T, 'b-', linewidth=2.5, label='Numerical Solution')
plt.plot(t, T_env + (T0 - T_env) * np.exp(-k * t), 'r--', 
         linewidth=2, label='Analytical Solution')
plt.xlabel('Time (min)', fontsize=12)
plt.ylabel('Temperature (° C)', fontsize=12)
plt.title('Cooling Law: SciPy solve_ivp', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('第 2 章-一阶微分方程/fig11_scipy.png', dpi=150, bbox_inches='tight')
plt.close()

方向场可视化

def plot_direction_field(f, x_range, y_range, ax=None, density=20):
    """绘制一阶 ODE 的方向场"""
    if ax is None:
        fig, ax = plt.subplots(figsize=(10, 8))
    
    x = np.linspace(*x_range, density)
    y = np.linspace(*y_range, density)
    X, Y = np.meshgrid(x, y)
    
    # 计算斜率
    dY = f(X, Y)
    dX = np.ones_like(dY)
    
    # 归一化
    N = np.sqrt(dX**2 + dY**2)
    dX, dY = dX/N, dY/N
    
    ax.quiver(X, Y, dX, dY, N, cmap='coolwarm', alpha=0.6)
    return ax

# 例： dy/dx = sin(x) - y
def f(x, y):
    return np.sin(x) - y

fig, ax = plt.subplots(figsize=(10, 8))
plot_direction_field(f, (-2*np.pi, 2*np.pi), (-3, 3), ax)

# 绘制几条解曲线
from scipy.integrate import odeint

x_span = np.linspace(-2*np.pi, 2*np.pi, 300)
for y0 in [-2, -1, 0, 1, 2]:
    sol = odeint(lambda y, x: np.sin(x) - y, y0, x_span)
    ax.plot(x_span, sol, 'k-', linewidth=1.5, alpha=0.7)

ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title("Direction Field: dy/dx = sin(x) - y", fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('第 2 章-一阶微分方程/fig12_direction_field.png', dpi=150, bbox_inches='tight')
plt.close()

总结

本章我们系统学习了一阶微分方程的主要类型和求解方法：

类型	标准形式	求解方法
可分离方程		分离变量后积分
一阶线性方程		积分因子
恰当方程	，	找势函数
伯努利方程		变量替换

应用领域： - 复利与投资（指数增长） - 药物代谢（指数衰减） - 种群动力学（ Logistic 方程） - 电路分析（ RC 电路） - 化学反应（混合问题）

数值方法： - 欧拉方法：简单但精度低 - RK4 方法：工程中最常用 - SciPy：现代 Python 工具

练习题

基础题

1. 求解可分离方程，。

2. 求解线性方程。

3. 验证 是否恰当，若不恰当求积分因子。

4. 一个水池初始有 1000 升水含 50kg 盐。纯水以 10 升/分钟流入，混合液以 10 升/分钟流出。求盐量。

5. RC 电路中，，初始电容电压为 0，接入 10V 电源。求。

进阶题

6. 求解伯努利方程。

7. 证明：对于 型方程，令 可化为可分离方程。

8. 讨论， 的解的存在性和唯一性。

9. 细菌数量满足 Logistic 方程，。求： - 的表达式

   • 何时？
   • 增长最快的时刻

10. 求解方程（提示：用极坐标）

编程题

11. 用 Python 实现改进的欧拉方法（ Heun's method），并与欧拉方法比较精度。

12. 对 Logistic 方程绘制方向场，并叠加不同初值的解曲线。

13. 模拟连续服药问题：每 6 小时服药 500mg，药物半衰期 4 小时，绘制血药浓度曲线。

14. 实现自适应步长的 RK45 方法，并与固定步长方法比较效率。

参考资料

1. Boyce, W. E., & DiPrima, R. C. (2012). Elementary Differential Equations and Boundary Value Problems. Wiley.
2. Zill, D. G. (2017). A First Course in Differential Equations with Modeling Applications. Cengage.
3. Edwards, C. H., & Penney, D. E. (2014). Differential Equations and Boundary Value Problems. Pearson.
4. SciPy Documentation: scipy.integrate

---

下一章预告：《高阶线性微分方程》 — 弹簧振动、 RLC 电路、谐振现象