ODE Solvers in MATLAB中ODE求解器的基本选择

博客文章：ODE Solvers in MATLAB中ODE求解器的基本选择

Ordinary Differential Equations (ODEs)

常微分方程（ODE，Ordinary Differential Equations）一个或者多个变量及其针对单个自变量的导数的方程。ODE在物理、工程、经济学等领域中广泛应用。通常，这一个自变量是时间，当然也可能是空间或其他变量。这个导数的阶数可以是任意的，通常是整数。

一般而言，$y^{\prime }$表示$y$对自变量$t$的导数，$y^{\prime \prime }$表示二阶导数，依此类推。

举例，下面是一个二阶的ODE：

$$y^{\prime \prime }=9y$$

在为微分方程求解时，一类是初值问题（IVP，Initial Value Problem），即给定初始条件$y(t_0)=y_0$和$y^{\prime }(t_0)=y_0^{\prime }$，求解在$t_0$附近的解。通常，用数值方法得到的结果表达在离散的时间点上：

$$\begin{gathered} t=[t_0,t_1,\dots ,t_n] \\ y=[y_0,y_1,\dots ,y_n] \end{gathered}$$

另一类是边值问题（BVP，Boundary Value Problem），即给定边界条件$y(t_0)=y_0$和$y(t_n)=y_n$，求解在$t_0$到$t_n$之间的解。

Matlab提供了多种ODE求解器，适用于不同类型的ODE问题。

ODE分类

Matlab的ODE求解器主要求解如下类型的一阶ODE：

• 显式ODE，形如：$y^{\prime }=f(t,y)$
• 线性隐式ODE，形如：$M(t,y)y^{\prime }=f(t,y)$，这里的$M(t,y)$是一个满秩的矩阵，称为质量矩阵。质量矩阵可能依赖时间，也可能依赖状态变量，也可以是一个常数矩阵。这个矩阵中的导数$y^{\prime }$与质量矩阵形成线性关系，构成隐式ODE。
  线性隐式ODE可以转化为显式ODE的形式。$y^{\prime }=M^{-1}(t,y)f(t,y)$，其中$M^{-1}(t,y)$是质量矩阵的逆矩阵。直接求解线性隐式ODE可以避免直接求逆矩阵的计算，特别是这种计算可能非常复杂或不稳定的情况。
• 如果在某些其工况下，方程组的某个方程中不包含$y^{\prime }$，则成为微分代数方程（DAE，Differential-Algebraic Equation）。DAE的求解器通常需要更多的计算资源和时间，因为它们需要同时处理代数方程和微分方程。
• 全隐式ODE，形如$f(t,y,y^{\prime })=0$，这种形式的ODE通常需要使用更复杂的数值方法进行求解。

Matlab对几种特殊类型ODE的支持

微分方程组

Matlab的ODE求解器可以处理多个变量的微分方程组。对于一个$n$维的微分方程组，可以将其写成向量形式：

$$\left[\begin{array}{c}{y}_1^{\prime }\\ y_2^{\prime }\\ ⋮\\ y_n^{\prime }\end{array}\right]=\left[\begin{array}{c}{f}_1(t,y_1,y_2,\dots ,y_n)\\ f_2(t,y_1,y_2,\dots ,y_n)\\ ⋮\\ f_n(t,y_1,y_2,\dots ,y_n)\end{array}\right]$$

写成向量形式后，可以使用Matlab的ODE求解器直接求解。

function dydt = odefun(t, y)
    dydt = zeros(2, 1); % 假设有两个变量
    dydt(1) = y(2); % y1' = y2
    dydt(2) = -9 * y(1); % y2' = -9 * y1
end

高阶ODE

Matlab的ODE求解器主要处理一阶ODE，但高阶ODE可以通过引入新的变量将其转化为一阶ODE组来求解。例如，考虑一个二阶ODE：

$$y^{\prime \prime }=-9y$$

可以引入一个新的变量$v=y^{\prime }$，将其转化为一阶ODE组：

$$\left[\begin{array}{c}{y}^{\prime }\\ v^{\prime }\end{array}\right]=\left[\begin{array}{c}v\\ -9y\end{array}\right]$$
然后使用Matlab的ODE求解器进行求解。

复数ODE

复数ODE可以按照照实部和虚部分别处理的方法进行求解。将复数变量$y=y_r+i{y}_i$，其中$y_r$是实部，$y_i$是虚部。然后将复数ODE转化为实数ODE组：

$$y^{\prime }=f(t,y)$$

可以转化为：

$$\left[\begin{array}{c}{y}_r^{\prime }\\ y_i^{\prime }\end{array}\right]=\left[\begin{array}{c}\text{Re}\left(f(t,y_r+i{y}_i)\right)\\ \text{Im}\left(f(t,y_r+i{y}_i)\right)\end{array}\right]$$

例如：

function f = complexf(t, y)
    f = y.*t + 2 * i;
end

写成ODE函数：

function fv = imaginaryODE(t, yv)
    y = yv(1) + 1i * yv(2); % 复数变量
    f = complexf(t, y); % 计算复数函数值
    fv = [real(f); imag(f)]; % 返回实部和虚部
end

在求解时，需要将初始条件也转化为实部和虚部：

y0 = 1 + i;
yv0 = [real(y0); imag(y0)];
[t, yv] = ode45(@imaginaryODE, [0 1], yv0);
y = yv(:, 1) + 1i * yv(:, 2); % 将结果转化为复数形式

刚性、非刚性

其实，对于ODE的求解来说，上面的分类是很容易理解的，但并不是特别重要。更重要的是ODE的刚性（Stiffness）问题。

刚性并没有一个严格的数学定义，但通常指的是ODE中存在多个相差较大的（时间）尺度，导致某些数值方法在求解时需要非常小的步长以保持稳定性。刚性ODE通常需要使用隐式方法进行求解。

举个例子，考虑如下的ODE，van der Pol方程。

$$\begin{gathered} y_1^{\prime }=y_2 \\ {y}_2^{\prime }=\mu (1-y_1^2)y_2-y_1 \end{gathered}$$

当$\mu$很大时，这个ODE就是刚性的。后面，还会专门介绍这个方程以及它的求解。

在Matlab中选择求解器时，首先就是要考虑ODE的刚性问题。

求解器基本选择

Matlab提供了多种ODE求解器，适用于不同类型的ODE问题。以下是一些常用的ODE求解器及其适用场景。

非刚性ODE求解器

• ode45：这是最常用的非刚性ODE求解器，基于Dormand-Prince方法（Runge-Kutta方法的一种）。适用于大多数非刚性问题，具有良好的精度和效率。
• ode23：适用于精度要求不高的非刚性问题，计算速度较快。
• ode113：基于多步法的求解器，适用于精度要求较高的非刚性问题，对于ODE函数计算代价较高时常用。
• ode78：高阶Runge-Kutta方法，适用于需要高精度的非刚性问题，但计算成本较高。
• ode89：更高阶的Runge-Kutta方法，适用于需要极高精度的非刚性问题，但计算成本更高。当需要积分非常精确的光华结果或者积分时间长度非常长时，可以考虑使用。

刚性ODE求解器

• ode15s：这是最常用的刚性ODE求解器，基于数值微分代数方程（DDAE）的方法。适用于大多数刚性问题，具有良好的稳定性和效率。
• ode23s：适用于精度要求不高的刚性问题，计算速度较快。
• ode23t：当需要求解隐式ODE或微分代数方程（DAE）时使用，适用于刚性问题。精度较低。
• ode23tb：基于Trapezoidal rule和Backward differentiation formula的方法，适用于刚性问题。

其他特殊ODE求解器

• ode15i：用于求解全隐式ODE，适用于刚性问题，但需要提供初始条件的导数值。

总结

总之，先把问题的微分方程转化为一阶向量的形式，然后首先用ode45试试，如果不行，再考虑刚性问题，换用ode15s。如果还不行，再考虑更高阶的求解器，或者更复杂的ODE形式。