【原】14.4 一阶微分方程

14.4 一阶微分方程

在现实生活中，我们可能希望对变量随时间变化的情况进行建模或定量表示，尤其是涉及到生物、电气或机械系统的情况。如果这些变化是连续的，通常可以使用涉及依赖变量导数的方程来表示系统。这样的方程称为微分方程（Differential Equations）。建模的主要目的是能够编写一组微分方程（DEs），以尽可能准确地描述所研究的系统。很少有微分方程可以通过解析方法求解，因此需要使用数值方法。在本节中，我们将介绍最简单的数值解法：欧拉方法（Euler's method，读作“欧拉”，音近于'boiler’）。我们还会简要探讨如何改进它。

微分方程是描述变量随时间变化的数学工具，它在物理学、化学、工程学等领域得到广泛应用。通过将系统拆分成微小的时间步长，欧拉方法基于初始条件和微分方程的导数信息来逼近连续系统的解。其基本思想是通过将连续的时间域离散化为离散的时间步长，然后根据微分方程求得在每个时间步长上的变量的近似值。

然而，欧拉方法有时会引入较大的近似误差，尤其对于具有复杂动力学行为的系统。为了改进欧拉方法的准确性，科学家和工程师们开发了更高阶的数值方法，如改进的欧拉方法、龙格-库塔方法等。这些方法使用更复杂的算法来逼近微分方程的解，并提供更精确的结果。

---

14.4.1 欧拉法

一般来说，我们想要求解一阶 ode (严格来说是普通ode)的形式:

欧拉数值求解该微分方程的方法是将dy/dx替换为牛顿商，使微分方程变为：

稍微重新排列一下后，我们得到

在科学和工程中，数值求解微分方程是一个非常重要且常见的问题，因此在这里介绍一些通用的符号表示是值得的。假设我们想要在区间 （通常情况下 a = 0）到 上进行微分方程的积分。我们将这个区间分为长度为 的 个步长，因此有

（这与第11章的更新过程中使用的符号表示方式相同，只是这里用更通用的 替代了之前的 ）。

为了与 MATLAB 的下标表示一致，我们定义 为 （在步骤 开始时的欧拉估计），其中 。则在步骤 结束时，有 。我们可以通过迭代方案将方程 替换为：

其中 。从第11章的回顾中可以知道，这种符号表示使得我们能够生成一个向量 ，并且可以对其进行绘图。同时请注意方程 与第11章中表示更新过程的方程之间的惊人相似性。这种相似性并非巧合。更新过程通常由微分方程建模，而欧拉方法为这类微分方程提供了一个近似解。

方程 表明，我们可以通过使用前一步的估计值 和微分方程 计算出下一步的估计值 。这种迭代的过程在每个步骤中不断更新估计值，逼近微分方程的解。欧拉方法是一种简单而直接的数值求解微分方程的方法，它在实践中非常有用。

14.4.2 示例：细菌生长

假设一个由1000个细菌组成的群落以每个个体每小时产生平均0.8个后代的速率 不受限制地持续增长。在经过10小时后，群落中会有多少个细菌？我们可以使用微分方程来模拟这种增长，微分方程表示为：

其中 表示时间 时的群落大小。这个过程称为指数增长。上述方程可以通过解析解得到著名的指数增长公式 .

为了对方程 进行数值求解，我们对其应用欧拉算法得到其中初始值 。

编写欧拉方法非常简单。下面的脚本实现了方程, 其中步长 设置为 0.5。它还计算了精确解以进行比较。

h = 0.5;
r = 0.8;
a = 0;
b = 10;
m = (b - a) / h;
N = zeros(1, m+1);
N(1) = 1000;
t = a:h:b;
for i = 1:m
    N(i+1) = N(i) + r * h * N(i);
end
Nex = N(1) * exp(r * t);
format bank
disp([t' N' Nex'])
plot(t, N), xlabel('Hours'), ylabel('Bacteria')
hold on
plot(t, Nex)
hold off


             0       1000.00       1000.00
          0.50       1400.00       1491.82
          1.00       1960.00       2225.54
          1.50       2744.00       3320.12
          2.00       3841.60       4953.03
          2.50       5378.24       7389.06
          3.00       7529.54      11023.18
          3.50      10541.35      16444.65
          4.00      14757.89      24532.53
          4.50      20661.05      36598.23
          5.00      28925.47      54598.15
          5.50      40495.65      81450.87
          6.00      56693.91     121510.42
          6.50      79371.48     181272.24
          7.00     111120.07     270426.41
          7.50     155568.10     403428.79
          8.00     217795.33     601845.04
          8.50     304913.47     897847.29
          9.00     426878.85    1339430.76
          9.50     597630.40    1998195.90
         10.00     836682.55    2980957.99

这个脚本首先设置了变量和数组的初始值，然后通过循环计算欧拉方法得到的近似解 N。通过指数增长公式，脚本计算了对应的精确解 Nex。最后，使用 disp 函数输出了时间、近似解和精确解的值，并使用 plot 函数绘制了时间和细菌数量的图像。

(蓝色)欧拉法; (红色)精确解

结果显示在表和图中。欧拉方法得到的解并不太好。实际上，误差在每一步都会变得更糟，经过10个小时的细菌时间后，误差约为72%。如果将步长 设得更小，数值解将会改善，但总会存在某个时间点 上误差超出可接受的限度。

在某些情况下，欧拉方法的表现可能会比这里好一些，但还有其他数值方法始终比欧拉方法表现更好。其中两种方法将在下面讨论。更复杂的方法可以在大多数数值分析的教科书中找到。然而，只要你意识到可能会出现误差，欧拉方法始终可以作为第一次近似解。

14.4.3 可选的下标符号

方程 是一种有限差分方案的示例。常规的有限差分符号表示中，初值通常表示为 ，即下标为 。 则表示在第 步结束时的估计值。如果你想要 MATLAB 中的下标与有限差分表示中的下标相同，那么初值 必须表示为 MATLAB 的标量 ，并且在 for 循环开始之前需要单独计算 。由于初始值不再包括在 MATLAB 的向量 t、N 和 Nex 中（现在只有 m 个元素，而不是 m + 1 个元素），你还需要单独显示或绘制初始值。以下是一个完整的脚本，使用有限差分表示生成欧拉解：

h = 0.5;
r = 0.8;
a = 0;
b = 10;
m = (b - a) / h;
N = zeros(1, m); % 现在少了一个元素
N0 = 1000;
N(1) = N0 + r * h * N0; % 不再是'自启动’
for i = 2:m
    N(i) = N(i-1) + r * h * N(i-1); % 有限差分表示
end
t = a+h:h:b; % 排除初始时间为 a
Nex = N0 * exp(r * t);
disp([a N0 N0]) % 单独显示初始值
disp([t' N' Nex'])
plot(a, N0) % 单独绘制初始值
hold on
plot(t, N), xlabel('Hours'), ylabel('Bacteria')
plot(t, Nex), hold off

这个脚本使用有限差分表示生成欧拉方法的解。它首先设置了变量和数组的初始值，然后使用循环根据有限差分表示计算了欧拉方法的近似解 N。初始值 的计算不再自启动，需要单独计算 。在生成时间数组 t 时，排除了初始时间 a。精确解 Nex 的计算与以前相同。脚本通过 disp 函数单独显示了初始值的内容，并通过 plot 函数单独绘制了初始值的图像。整个脚本生成了使用有限差分表示的欧拉方法解，并将其与精确解进行了比较。

14.4.4 预测校正法

预测校正法（predictor-corrector method）是一种数值方法，用于求解微分方程的初值问题。它的基本思想是先使用一个预测步骤估计下一个时间步长的值，然后使用校正步骤来修正这个估计值，得到更准确的近似解。

在预测步骤中，预测器根据已知的点和近似解的斜率来计算下一个时间步长的近似值。通常情况下，欧拉方法或其他较简单的数值方法用于生成这个预测值。

在校正步骤中，校正器采用预测步骤得到的近似值，根据更精确的斜率计算修正值。这个更精确的斜率可以使用更精细的数值方法获得，例如改进的欧拉方法或龙格-库塔方法。

预测校正法的一般步骤如下：

1. 初始化时间步长、初值和时间网格。
2. 在预测步骤中，使用一个简单的数值方法（如欧拉方法）估计下一个时间步长的近似值。
3. 在校正步骤中，根据更精确的数值方法（如改进的欧拉方法或龙格-库塔方法）计算修正值。
4. 使用修正值更新近似解，并更新时间网格。
5. 重复步骤2-4，直到达到所需的时间点或时间步数。

预测校正法的优点是可以提高数值解的精确性，尤其在需要较高精度的情况下。它综合使用不同的数值方法，结合它们的优点来得到更准确的解。然而，它也更复杂，需要更多的计算量。在实际应用中，选择合适的预测器和校正器以及适当的步长是非常重要的，以确保稳定性和准确性。

---

对于一阶微分方程 , 已知，数值求解的一个改进是使用预测校正法。我们可以将欧拉逼近表示为：

但是这个公式更偏向于在右侧计算 时使用旧的 值。因此，我们可以采用以下公式：

，其中，因为这个公式在计算右侧的 时也涉及新的值吗？当然，问题在于目前还未知，所以我们不能在方程右侧使用它。但是我们可以使用欧拉方法估计（预测）从方程得到的，然后使用方程纠正预测值，计算出一个更好的，我们将其称为 。因此，完整的过程如下：

重复以下步骤直到需要的次数为止：使用欧拉方法进行预测：然后使用公式对进行修正：

这被称为预测校正法。下面的脚本可以很容易地适应这个问题。相关的代码行如下：

for i = 1:m % m个长度为dt的步骤
    ne(i+1) = ne(i) + r * h * ne(i);
    np = nc(i) + r * h * nc(i);
    nc(i+1) = nc(i) + r * h * (np + nc(i))/2;
    disp([t(i + 1) ne(i + 1) nc(i + 1) nex(i + 1)])
end;

其中，ne 代表未校正的欧拉解，np是欧拉预测器（由于这是一个中间结果，因此np不需要向量），nc是校正器。最大误差现在只有15％。这比未校正的欧拉解好得多，尽管仍有改进的空间。

---

以下是一个完整的MATLAB代码示例，用于使用预测校正法求解一阶微分方程：

h = 0.5;
r = 0.8;
a = 0;
b = 10;
m = (b - a) / h;
t = a:h:b;
ne = zeros(1, m + 1); % 未校正的欧拉解
nc = zeros(1, m + 1); % 校正器

% 设置初值
ne(1) = 1000;
nc(1) = 1000;

% 预测校正法的主循环
for i = 1:m
    % 预测器
    np = ne(i) + r * h * ne(i);
    
    % 校正器
    nc(i+1) = ne(i) + h * (f(t(i+1), np) + f(t(i), ne(i))) / 2;
    
    % 更新未校正欧拉解
    ne(i+1) = ne(i) + r * h * ne(i);
    
    % 显示结果
    disp([t(i+1) ne(i+1) nc(i+1)])
end

% 函数：f(x, y)
%f = @(x, y) x^3; % 根据具体问题的 f(x, y) 函数定义进行修改

% 方程右端的 f(x, y) 函数定义
%function dy = f(x, y)
%    dy = x^3; % 根据具体问题的 f(x, y) 函数定义进行修改
%end

在此代码中，我们使用h表示步长，r表示方程中的系数，a和b表示所求解的时间范围。m是根据步长计算的步数。ne和nc分别表示未校正的欧拉解和校正器的结果。

我们首先设置初始值ne(1)和nc(1)。然后，使用预测校正法的主循环，从i=1到i=m进行迭代。在每次迭代中，我们使用预测器np来预测未校正欧拉解的下一个值。然后，使用校正器nc(i+1)来纠正预测值，并计算出更准确的校正解nc(i+1)。最后，我们使用未校正欧拉解的公式更新未校正欧拉解的下一个值。