MATLAB 编程/高级主题/数值运算/常微分方程

所有微分方程都具有相同的语法，您必须使用该语法，以及相同的输入和输出参数。所有求解器都需要三个输入参数：您要解决的微分方程的函数句柄、要积分的时间间隔和一组初始条件。假设我们要解决简单的微分方程 y' = y，y(0) = 1，其真解为 y(t) = e^t。假设我们希望在 t = 0 到 t = 1 之间得到解。

要使用函数句柄，您必须首先创建一个包含该函数的 M 文件，如下所示

function Yprime = func(t, y)
Yprime = 0;
Yprime = y;

请注意，即使时间参数未在微分方程中使用，您也必须包含它。将 Yprime 初始化为一个零数组将在您尝试解决多个函数时为您节省麻烦；Yprime 必须作为垂直数组返回，但是，如果您没有将其初始化为垂直数组（或在最后进行转置），它默认情况下会返回一个水平数组。

还要注意，除了 ode15i 之外，必须显式地针对 y' 求解该函数。

创建此文件后，按顺序（函数、时间间隔、初始条件）使用参数调用 ODE 函数

>> ode45(@func, [0,1], 1)

另一种方法是使用匿名函数，这仅在您只有一个函数时有用（否则您必须使用函数句柄）。您还必须根据因变量和时间来声明匿名函数

>> func = @(t,y) y;
>> ode45(func, [0,1], 1)

在没有输入参数的情况下调用 ODE 函数会生成解的图形。使用一个输出参数调用它会返回一个结构数组

>> Struct = ode45(func, [0,1], 1)
 Struct = solver: 'ode45'
          extdata: [1x1 struct]
          x: [1x11 double]
          y: [1x11 double]
          stats: [1x1 struct]
          idata: [1x1 struct]

此数据主要用于将来调用“deval”函数以在您想要的时间获取答案

>> Solution = deval(Struct, 0.5)
 Solution = 1.6487

Deval 还可以通过包含第二个输出参数来返回感兴趣点的导数。

>> [Solution, Derivative] = deval(Struct, 0.5)
Solution = 1.6487
Derivative = 1.6487

由于 e^x 的导数是其本身，因此导数和解在这里相同是有道理的。

使用两个输出参数调用 ode45 会返回两个数据列表；首先是 t，然后是适当大小矩阵中的自变量。

ODE 选项

MATLAB 为您提供了相当多的选项来修改它求解微分方程的方式。帮助文件对所有选项都做了很好的描述，因此这里将不再详细介绍。要获取列表，请使用帮助函数

>>help odeset

要获取不同选项名称的列表以及要传递给它的内容，只需在命令提示符中键入“odeset”。它返回一个数据类型或一个有限的选项列表。它还以括号形式列出了所有选项的默认值。

要设置特定选项或选项列表，请先键入选项名称，然后键入您想要的选项值。例如，假设您要将默认相对容差从 10^-3 调整为 10^-4。您将按如下方式调用“odeset”

>> option = odeset('RelTol', 10^-4);

请注意，选项名称必须作为字符串传递，否则您很可能会收到“未定义变量”错误。可以通过将所有选项都放在一个列表中来一次传递多个选项

>> option = odeset('RelTol', 10^-4, 'AbsTol', 10^-7, 'MassSingular', 'no');

然后可以将选项结构作为第四个（可选）输入参数传递给 ode 函数

>> [T,X] = ode45(@func, [0,1], 1, option);

这将返回比默认值更准确的值，因为误差容差更严格。它还应更快地计算，因为 MATLAB 不会检查这是否是微分代数方程（这就是 MassSinglular 选项的作用；它通常设置为“maybe”，因此 MATLAB 会自行检查）。

ODE 求解器

MATLAB 并不只有一种 ODE 求解器，截至 MATLAB 7.1 版本，它有八种。有些求解器更适合解决某些问题，而不是其他问题，这就是为什么所有这些求解器都被包含在内的原因。MATLAB 帮助文件中列出了每个函数可以执行的操作，但这里提供一个快速摘要，大致按照您应该尝试它们的顺序（除非您已经知道问题的分类）

大多数问题都可以使用 ode45 求解，由于这是速度和精度的最佳折衷方案，因此它应该是您首先使用的求解器。

如果您需要非常严格的误差容差或大量数据点，请使用 ode113。

如果您需要相对宽松的误差容差或问题在 ode45 中很慢，请尝试使用 ode23。

如果问题确实是刚性的，而 ode45 失败，请对于宽松的容差使用 ode23tb，对于稍微严格的容差（具有恒定质量矩阵）使用 ode23s，对于更严格的容差或非恒定质量矩阵使用 ode15s。

如果您需要在刚性问题上进行数值阻尼的解，请使用 ode23t。

s 表示该算法适用于刚性问题。

还有另一种 ODE 求解器，它是特殊的

只能使用 ode15i 求解隐式问题。

由于 ode15i 在语法方面有一些区别，因此将在下一节中进行讨论。

隐式 ODE：ODE15i

由于 ode15i 是唯一一个求解隐式方程的 ODE 求解器，因此它在如何输入函数方面必须有一些特殊的语法规则。

如果您使用的是 ode15i，请按如下方式声明该函数

function Z = func(t,y,Yprime)
Z = 0;
Z = y - Yprime;

请注意，对于 ode15i，您必须将函数置于标准形式（针对 0 求解），而对于所有其他 ODE 函数，您必须显式地针对 y' 求解。还要注意，您必须声明三个输入参数，而不是通常的两个参数。

当您调用 ode15i 时，您不仅必须包含 y 的初始条件，还必须包含 Yprime 的初始条件。Yprime 的初始条件进入第四个参数

>> [t,x] = ode15i(@func, [0,1], 1, 1);

这将返回与 ode45 相似的结果（用于显式方程），但数据点更少。

选项结构作为可选的第五个参数传递给 ode15i。ode15i 的所有输出选项与其他 ODE 求解器相同。