14.6 龙格-库塔法
有多种算法可以用于求解一阶常微分方程组，这些算法统称为龙格-库塔（Runge-Kutta）方法。这些方法的公式比较复杂，可以在大多数数值分析的书籍中找到。然而，正如您可能猜到的那样，MATLAB 已经提供了许多用于求解常微分方程的求解器，这些求解器在MATLAB的帮助文档 "MATLAB Help: Mathematics: Differential Equations" 中有详细介绍。其中包括 ode23（二阶/三阶）和 ode45（四阶/五阶），它们实现了龙格-库塔方法。（数值方法的阶数是在主导误差项中 h（即 dt）的幂次。由于 h 通常很小，幂次越高，误差越小）。我们将在这里演示使用 ode23 和 ode45 的用法，首先是使用单个一阶常微分方程，然后是使用方程组。
下面是在MATLAB中使用 ode23 和 ode45 求解一阶常微分方程的示例：
% 定义常微分方程函数 dy/dt = f(t, y)
f = @(t, y) t^2 - y;
% 定义初始时间 t0 和初始条件 y0
t0 = 0;
y0 = 1;
% 定义求解时间段 [tstart, tend]
tstart = 0;
tend = 1;
% 使用 ode23 求解方程
[t, y] = ode23(f, [tstart, tend], y0);
% 显示结果
disp([t, y]);
% 使用 ode45 求解方程
[t, y] = ode45(f, [tstart, tend], y0);
% 显示结果
disp([t, y]);
在这个示例中，我们定义了一个一阶常微分方程 dy/dt = t^2 - y，并使用 ode23 和 ode45 分别对方程进行了求解。ode23 是一个二阶/三阶的龙格-库塔方法，而 ode45 是一个四阶/五阶的龙格-库塔方法。我们指定了初始时间 t0、初始条件 y0，以及求解的时间段 [tstart, tend]。然后，我们通过调用 ode23 和 ode45 函数来求解方程，并将结果保存在结果向量 t 和 y 中。最后，我们使用 disp 函数显示求解的结果。
以上示例演示了如何使用 ode23 和 ode45 求解一阶常微分方程。对于求解更复杂的方程组，可以将其转化为一组一阶方程，并使用相同的方法进行求解（传入一个向量函数 f 即可）。
14.6.1单个微分方程
使用ode23求解细菌生长问题的步骤如下，方程式为：
- 首先，编写一个函数文件来解决微分方程右侧的函数。该函数需要以 和 为输入变量（即微分方程的自变量和因变量），按照这个顺序创建函数文件
f.m，示例如下：
function y = f(t, Nr)
y = 0.8 * Nr;
a = 0;
b = 10;
n0 = 1000;
[t, Nr] = ode23(@f, [a:0.5:b], n0);
[a:0.5:b]: 一个向量(tspan)，指定积分范围。如果tspan有两个元素([a b])，求解器返回在每个积分步骤处评估的解（求解器选择积分步骤并可能对其进行变化）。这种形式适合绘图。然而，如果您想要在固定的时间间隔显示解，就像我们在这里想要的那样，请使用包含三个元素的tspan形式。然后，在tspan中的每个时间点返回求解结果。解的精确度不受所使用的tspan形式的影响。
- 输出参数是两个向量：时间点 t 处的解 Nr。使用ode23求解10小时后的细菌数量为2961338个。从表14.1的精确解中可以看出，这里的误差仅为0.7%。
如果您从ode23获得的解不够准确，您可以使用额外的可选参数请求更高的准确度。请参考帮助文档。
如果您需要更准确的数值解，您可以使用ode45。它给出的细菌最终数量为2981290，误差约为0.01%。
14.6.2 微分方程组:混沌
天气预测具有不确定性和变化性的主要原因，不再被认为是系统本身的复杂性，而是对其建模的微分方程的性质。这类微分方程被称为混沌方程。当初始条件微小变化时，这些方程将产生完全不同的结果。换句话说，准确的天气预测在很大程度上依赖于对初始条件的精确测量。
在1961年，研究气象学的爱德华·洛伦兹发现了这种现象。虽然他最初的方程过于复杂，不便在这里讨论，但下面这个简化的系统具有相同的混沌特性：
可以使用MATLAB的ODE求解器非常容易地求解这个微分方程组。基本思想是使用一定的初始条件求解微分方程组，绘制解的图形，然后轻微改变初始条件，并将新的解与旧解叠加在一起，以观察它们之间的差异程度。
通过这种方法，研究人员可以更好地理解混沌系统的行为。这种对初始条件敏感的特性使得天气预报变得困难，因为即使有很小的误差，初始条件的微小变化也可能导致完全不同的天气演变。因此，在天气预报中准确测量和预测初始条件变得至关重要。
我们首先求解初始条件为 , , 的方程组。
- 写一个函数文件
lorenz.m来表示系统的右手边，如下所示：
function f = lorenz(t, x)
f = zeros(3,1);
f(1) = 10 * (x(2) - x(1));
f(2) = -x(1) * x(3) + 28 * x(1) - x(2);
f(3) = x(1) * x(2) - 8 * x(3) / 3;
MATLAB向量x的三个元素分别表示三个相关的标量变量x、y和z，即x(1)、x(2)和x(3)。向量f的元素表示三个微分方程的右侧。当一个向量被这样的微分方程函数返回时，它必须是一个列向量，因此有以下语句：
f = zeros(3,1);
这行代码将创建一个大小为3x1的列向量f，并初始化为零向量。这是为了确保f有正确的大小和格式，以便与微分方程函数一起使用。
- 使用以下命令求解从 t=0 到 t=10 的系统：
x0 = [-2 -3.5 21]; % 初始值向量
[t, x] = ode45(@lorenz, [0 10], x0);
plot(t,x)
注意，我们现在使用ode45，因为它更准确。你将会看到三个图形，分别对应着x、y和z（以不同的颜色表示）。
- 如果图像中只有一个图形，更容易观察更改初始值的影响。事实上，最好只绘制y(t)的解。
MATLAB的解x实际上是一个具有三列的矩阵（可以从whos命令看到）。我们想要的解y(t)将是第二列，所以使用以下命令单独绘制它：
plot(t,x(:,2),'g')
hold on
然后使用hold命令保留图像在坐标轴上。
现在我们可以观察更改初始值的影响。让我们只改变x(0)的初始值，从-2改为-2.04──这只是三个初始值中的一个，变化只有2%。使用以下命令进行这个操作，求解微分方程，并绘制y(t)的新图表（以不同的颜色表示）：
x0 = [-2.04 -3.5 21];
[t, x] = ode45(@lorenz, [0 10], x0);
plot(t,x(:,2),'r')
你会看到，在t约为1.5之前，两个图形几乎无法区分。随着t的增加，差异会逐渐增加，直到t约为6，解会突然以相反的方向发生巨大变化。随着t的进一步增加，新的解与旧的解毫无相似之处。
现在使用ode23对系统使用原始初始值进行求解：
x0 = [-2 -3.5 21];
[t,x] = ode23(@lorenz, [0 10], x0);
只绘制y(t)的图形，即x(:,2)，然后以不同的颜色绘制使用相同初始值的ode45解。
x0 = [-2.04 -3.5 21];
% 求解微分方程
[t, x] = ode45(@lorenz, [0 10], x0);
plot(t,x(:,2),'g')
hold on
x0 = [-2 -3.5 21];
[t,x] = ode23(@lorenz, [0 10], x0);
plot(t,x(:,2),'b')一个奇怪的现象发生了——在t>1.5时，解开始急剧偏离！初始条件是相同的，唯一的区别是龙格-库塔方法的顺序。
最后，使用 ode23s 求解系统，并将其解与 ode23s 的解叠加在一起（s代表“稳定”。对于一种稳定的微分方程，解的变化在积分区间内与时间相比非常短）。在t>5时，ode45和ode23s的解才开始分歧。
解释是ode23、ode23s和ode45都存在数值误差（如果能将它们与精确解进行比较的话——实际上无法找到精确解）。然而，这三种情况下的数值误差是不同的。这种差异与使用非常微小的不同初始值开始数值解的效果相同。
我们如何知道何时得到了“正确”的数值解？实际上我们不知道——我们所能做的就是提高数值方法的准确性，直到在感兴趣的区间内不再发生剧烈的变化。因此，在我们的例子中，只有对于 t<5 的解（使用ode23s或ode45）我们可以相当确定。如果这还不够好，你就需要找到更准确的微分方程求解器。
所以要小心：解“混沌”的微分方程非常棘手！
顺便说一下，如果你想看到著名的混沌“蝴蝶”图片，只需将 x 与 z 绘制为时间的函数（这样得到的图形被称为相位平面图）。以下命令可以实现这个目标：
plot(x(:,1), x(:,3))
你将会看到解的一个静态二维投影，也就是解随时间的演化。MATLAB Launch Pad中的演示示例可以让你动态地观察轨迹在三维中的演化（示例：图形：洛伦兹吸引子动画）。
14.6.3 向ODE求解器传递附加参数
在上述的MATLAB常微分方程（ODE）求解器的示例中，微分方程的右手边中的系数（例如，方程 中的值28）都是常数。在实际建模中，您很可能经常需要更改这些系数。为了避免每次更改系数时都需要编辑函数文件，您可以将这些系数作为额外的参数传递给ODE求解器，然后由求解器将它们传递给微分方程函数。以下以Lotka-Volterra捕食者-猎物模型为例，展示如何实现这一点：
首先，编写一个名为volterra.m的函数M文件，内容如下：
function f = volterra(t, x, p, q, r, s)
f = zeros(2,1);
f(1) = p*x(1) - q*x(1)*x(2);
f(2) = r*x(1)*x(2) - s*x(2);
然后，在命令行窗口中输入以下语句，生成在图中显示的特征性振荡图形：
p = 0.4; q = 0.04; r = 0.02; s = 2;
[t,x] = ode23(@volterra,[0 10],[105; 8],[],p,q,r,s);
plot(t, x)
Lotka-Volterra模型:(红色)捕食者;(蓝色)猎物注意：
- 额外的参数（p，q，r和s）必须跟随ODE求解器的第四个输入参数。如果没有设置选项参数（如我们的情况），可以使用[ ]作为选项参数的占位符。
- 现在可以在命令行窗口中更改系数，并获得新的解，而无需编辑函数文件。
