跳转至

MATLAB 常微分方程 (ODE) 求解

常微分方程 (Ordinary Differential Equation, ODE) 是描述一个或多个因变量关于单个自变量（通常是时间）的变化率的方程。它们是模拟物理系统、化学反应、生物种群动态等众多动态系统的基础。

MATLAB 提供了世界一流的 ODE 求解器套件，能够高效、准确地求解各种类型的初值问题 (Initial Value Problem, IVP)。

ODE 求解的基本流程

求解一个形如 \(y' = f(t, y)\) 且具有初始条件 \(y(t_0) = y_0\) 的 ODE 问题，在 MATLAB 中通常遵循以下三步：

  1. 定义函数: 创建一个 MATLAB 函数，用于计算在任意给定时间和状态下，系统的导数值 \(f(t, y)\)。此函数的签名必须是 dydt = odefun(t, y)。
  2. 设置求解选项: 使用 odeset 函数创建一个选项结构体，用于控制求解器的行为，如误差容限、步长等（此步可选，但推荐用于高级控制）。
  3. 调用求解器: 从 MATLAB 提供的求解器套件（如 ode45）中选择一个，传入定义好的函数、求解时间区间和初始条件。

数值求解器 (Numerical Solvers)

MATLAB 的 ODE 求解器都是基于数值方法的。它们并不寻找解析解，而是从初始点出发，通过一系列离散的时间步长，逐步计算出解在未来时间点上的近似值。

数学原理：龙格-库塔方法 (Runge-Kutta Methods)

大多数 MATLAB 求解器，特别是 ode45，都基于龙格-库塔 (RK) 族算法。其核心思想是在每个时间步内，通过在几个“探测点”计算导数值，并将它们进行加权平均，来获得一个比简单的欧拉法更精确的下一步估计。

例如，一个经典的四阶 RK 方法的迭代公式为：

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ k_4 &= f(t_n + h, y_n + hk_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} \]

ode45 使用的是一种更高级的嵌入式龙格-库塔 (4,5) 对（也称 Dormand-Prince 对）。它在每一步同时计算一个四阶和一个五阶的近似解。这两个解之间的差异可以用来估计局部截断误差。求解器会根据这个误差估计，自适应地调整步长 \(h\)：如果误差太大，就减小步长重新计算；如果误差很小，就增大步长以提高效率。

非刚性问题求解器 (Nonstiff Solvers)

适用于大多数常规问题，其中解的变化率与步长处于同一数量级。

ode45 - 通用非刚性问题求解器

  • 功能 首选，通用求解器。它通常是尝试解决任何 ODE 问题的第一个选择。它在大多数情况下都能在速度和精度之间取得最佳平衡。

  • 语法 

    [t, y] = ode45(odefun, tspan, y0)
[t, y] = ode45(odefun, tspan, y0, options)

  • 语法说明

    • odefun: 定义微分方程 \(y' = f(t, y)\) 的函数句柄。
    • tspan: 求解的时间区间，形如 [t0, tf]。
    • y0: 初始条件向量，即在 t0 时的 y 值。
    • options: 由 odeset 创建的选项结构体。
    • t: 返回的时间点列向量。
    • y: 返回的解矩阵，每一行对应一个时间点，每一列对应一个状态变量。

  • 示例: 洛伦兹吸引子 

    % 1. 定义洛伦兹方程组
% dx/dt = sigma * (y - x)
% dy/dt = x * (rho - z) - y
% dz/dt = x * y - beta * z
sigma = 10;
rho = 28;
beta = 8/3;

lorenz_fun = @(t, xyz) [sigma * (xyz(2) - xyz(1));
                       xyz(1) * (rho - xyz(3)) - xyz(2);
                       xyz(1) * xyz(2) - beta * xyz(3)];

% 2. 设置时间和初始条件
tspan = [0, 50];
y0 = [1, 1, 1];

% 3. 调用 ode45 求解
[t, y] = ode45(lorenz_fun, tspan, y0);

% 4. 绘图
figure;
plot3(y(:,1), y(:,2), y(:,3));
title('ode45 求解洛伦兹吸引子');
xlabel('x'); ylabel('y'); zlabel('z');

ode23, ode113 - 其他非刚性求解器

  • ode23: 基于龙格-库塔 (2,3) 对，精度低于 ode45，但对于误差容限要求不高的问题，计算成本更低。
  • ode113: 变阶 Adams-Bashforth-Moulton PECE 求解器。适用于对计算成本敏感且要求高精度的问题。

刚性问题求解器 (Stiff Solvers)

数学原理：刚性问题 (Stiff Problem)

一个 ODE 系统是刚性的，如果其解中包含变化速度差异极大的分量。例如，一个化学反应中，某些中间产物的半衰期可能是纳秒级的，而最终产物的变化是分钟级的。如果使用 ode45 这样的显式求解器，为了捕捉到纳秒级的快速变化，它会被迫采用极小的步长来维持数值稳定性，即使在解已经进入平滑变化的阶段也是如此。这会导致计算时间极长。

刚性求解器（如 ode15s）使用隐式算法。在每一步，它们需要求解一个非线性方程组来确定下一步的解，这使得每一步的计算成本更高。但作为回报，它们可以用比显式求解器大得多的步长来保持稳定，从而在求解刚性问题时总计算时间大大减少。

ode15s - 通用刚性问题求解器

  • 功能 首选，刚性问题通用求解器。它是一个变阶求解器，基于数值微分公式 (NDFs)。当你怀疑一个问题是刚性的，或者 ode45 求解非常缓慢时，应首先尝试 ode15s。

  • 语法 (与 ode45 完全相同) 

    [t, y] = ode15s(odefun, tspan, y0, options)

  • 示例: 范德波尔方程 (刚性模式) 
    % 当 mu 很大时，范德波尔方程是刚性的
mu = 1000;
vdp_fun = @(t, y) [y(2); mu*(1-y(1)^2)*y(2) - y(1)];

tspan = [0, 3000];
y0 = [2; 0];

% 使用 ode15s 求解
tic; % 开始计时
[t, y] = ode15s(vdp_fun, tspan, y0);
toc; % 结束计时

figure;
plot(t, y(:,1));
title('ode15s 求解刚性范德波尔方程');
xlabel('时间 t');
ylabel('解 y_1');
% 尝试用 ode45 求解此问题，你会发现它会耗时非常长

ode23s, ode23t, ode23tb - 其他刚性求解器

  • ode23s: 基于修正的二阶 Rosenbrock 公式，精度较低，但有时比 ode15s 更高效。
  • ode23t: 梯形规则求解器，适用于中度刚性问题。
  • ode23tb: 隐式龙格-库塔公式，效率可能高于 ode15s，但精度较低。

odeset - 设置求解器选项

  • 功能 创建或修改 ODE 求解器的选项结构体，用于精细控制求解过程。

  • 语法 

    options = odeset('Name1', Value1, 'Name2', Value2, ...)

  • 核心参数说明

    • 'RelTol': 相对误差容限（默认 1e-3）。控制解的相对精度。
    • 'AbsTol': 绝对误差容限（默认 1e-6）。控制解的绝对精度。当解接近于零时，此选项非常重要。
    • 'Jacobian': 提供雅可比矩阵 \( \partial f / \partial y \)。对于刚性问题，手动提供雅可比矩阵可以显著提高求解速度。
    • 'Events': 事件函数。用于检测解的特定事件（如穿过零点、达到某个阈值）并终止求解。
    • 'OutputFcn': 输出函数。在每个成功的积分步后调用，可用于实时绘图或数据处理。

  • 示例: 使用选项提高精度 

    lorenz_fun = @(t, xyz) [10 * (xyz(2) - xyz(1)); ...
                       xyz(1) * (28 - xyz(3)) - xyz(2); ...
                       xyz(1) * xyz(2) - (8/3) * xyz(3)];

% 创建高精度选项
opts = odeset('RelTol', 1e-8, 'AbsTol', 1e-11);

[t, y] = ode45(lorenz_fun, [0, 20], [1, 1, 1], opts);

figure;
plot3(y(:,1), y(:,2), y(:,3));
title('使用 odeset 提高求解精度');

符号求解器 (Symbolic Math Toolbox)

dsolve - 符号常微分方程求解器

  • 功能 尝试寻找常微分方程 (ODE) 的精确解析解。

  • 语法 

    sol = dsolve(eqn)
sol = dsolve(eqns, conds)

  • 语法说明

    • eqn(s): 一个或多个符号方程。使用 diff 函数来表示导数，使用 == 定义方程。
    • conds: 初始或边界条件。

  • 示例: 求解二阶线性 ODE 

    % 求解 y'' + y = 0, y(0)=1, y'(0)=0
syms y(t)

% 定义导数
Dy = diff(y, t);
D2y = diff(y, t, 2);

% 定义方程和初始条件
eqn = D2y + y == 0;
cond1 = y(0) == 1;
cond2 = Dy(0) == 0;

% 求解
ySol(t) = dsolve(eqn, [cond1, cond2]);

disp('ODE 的解析解是:');
disp(ySol(t)); % 结果是 cos(t)

% 绘图验证
fplot(ySol(t), [0, 2*pi]);
title('dsolve - 求解 y'''' + y = 0');