34、MATLAB 数值求解微分方程全解析

MATLAB 数值求解微分方程全解析

1. MATLAB 微分方程求解概述

MATLAB 提供了众多用于数值求解形如 $\frac{dy}{dt} = f(t, y)$ 这类常微分方程的函数。对于高阶微分方程或微分方程组，需将其重新表述为一阶表达式系统。不同的微分方程可能需要不同的求解技术，MATLAB 包含了多种微分方程求解器，它们具有相同的格式，方便通过更改函数名来尝试不同的求解方法。

常微分方程求解函数	适用问题类型	数值求解方法	注释
ode45	非刚性微分方程	龙格 - 库塔法	若对函数了解不多，是首选的初步猜测技术。使用显式龙格 - 库塔 (4, 5) 公式，即多曼 - 普林斯对。
ode23	非刚性微分方程	龙格 - 库塔法	使用博加基和尚皮内的显式龙格 - 库塔 (2, 3) 对。若函数 “轻度刚性”，可能比 ode45 更好。
ode113	非刚性微分方程	亚当斯法	与 ode45 和 ode23 这两种单步求解器不同，它是多步求解器。
ode15s	刚性微分方程和微分代数方程	NDFs (BDFs)	使用数值微分公式 (NDFs) 或向后微分公式 (BDFs)。难以预测哪种技术对刚性微分方程最有效。
ode23s	刚性微分方程	罗森布罗克法	改进的二阶罗森布罗克公式。
ode23t	中度刚性微分方程和微分代数方程	梯形法则	若需要无数值阻尼的解，很有用。
ode23tb	刚性微分方程	TR - BDF2	使用带梯形法则 (TR) 和二阶向后微分公式 (BDF2) 的隐式龙格 - 库塔公式。
ode15i	完全隐式微分方程	BDF	使用向后差分公式 (BDF) 求解形如 $f (y, y′, t ) = 0$ 的隐式微分方程。

每个求解器至少需要以下三个输入：
- 一个函数句柄，用于描述关于 $t$ 和 $y$ 的一阶微分方程或微分方程组。
- 感兴趣的时间跨度。
- 系统中每个方程的初始条件。

求解器都会返回 $t$ 和 $y$ 值的数组，格式如下：

[t,y] = odesolver(function_handle,[initial_time, final_time], [initial_cond_array])

如果不指定结果数组 [t,y] ，函数将绘制结果图。

2. 函数句柄输入

函数句柄是函数的 “昵称”，它可以指向存储为 M 文件的标准 MATLAB 函数，也可以指向匿名 MATLAB 函数。对于我们讨论的微分方程 $\frac{dy}{dt} = f(t, y)$，函数句柄等同于 $\frac{dy}{dt}$。

例如，对于单个简单的微分方程，匿名函数示例如下：

dydt = @(t,y) 2*t;

对应于 $\frac{dy}{dt} = 2t$。虽然此函数结果中未使用 $y$ 值，但 $y$ 仍需作为输入的一部分，才能被 ode 求解器接受。

如果要指定方程组，定义一个函数 M 文件可能更方便。函数的输出必须是一阶导数值的列向量，例如：

function dy=twofuns(t,y)
    dy(1) = y(2);
    dy(2) = -y(1);
    dy=[dy(1); dy(2)];
end

此函数表示方程组 $\frac{dy}{dt} = x$，$\frac{dx}{dt} = -y$，也可更紧凑地表示为 $y_1’ = y_2$，$y_2’ = -y_1$。

3. 求解问题

将感兴趣的时间跨度和每个方程的初始条件作为向量，与函数句柄一起输入到求解器方程中。以求解方程 $\frac{dy}{dt} = 2t$ 为例：

graph TD;
    A[定义匿名函数 dydt] --> B[指定时间跨度和初始条件];
    B --> C[选择求解器 ode45 求解];
    C --> D[得到 t 和 y 值数组或绘制结果图];

我们在前面创建了该常微分方程的匿名函数 dydt ，要在 $-1$ 到 $1$ 的区间内计算 $y$，并指定初始条件 $y(-1) = 1$。若对自己的方程或方程组特性不了解，首次尝试应使用 ode45 ：

[t,y] = ode45(dydt,[-1,1],1);

此命令返回 $t$ 值数组和对应的 $y$ 值数组。可以自行绘制这些值，或者若不指定输出数组，让求解器函数绘制结果：

ode45(dydt,[-1,1],1);

结果与解析解 $y = t^2$ 一致。

当输入函数或函数系统存储在 M 文件中时，语法略有不同。现有 M 文件的句柄定义为 @m_file_name 。例如，要求解 twofun 中描述的方程组，使用以下命令：

ode45(@twofun,[-1,1],[1,1]);

也可以为函数文件分配一个函数句柄，然后将其用作微分方程求解器的输入：

some_fun = @twofun;
ode45(some_fun,[-1,1],[1,1]);

感兴趣的时间跨度是从 $-1$ 到 $1$，初始条件均为 $1$。注意，系统中每个方程都有一个初始条件。

4. 求解高阶微分方程

ode 系列函数（如 ode45 或 ode23 ）用于求解单个一阶微分方程或一阶微分方程组。若需要求解高阶问题，可通过简单的代换将高阶微分方程表示为一系列方程。

例如，考虑方程 $\frac{d^2y}{dt^2} + \frac{dy}{dt} = y + t$，引入新变量 $z = \frac{dy}{dt}$，则 $\frac{dz}{dt} = \frac{d^2y}{dt^2}$。代入原方程得到 $\frac{dz}{dt} + \frac{dy}{dt} = y + t$，这是一个一阶微分方程。实际上，我们将原方程替换为以下两个方程：
$\frac{dy}{dt} = x$
$\frac{dz}{dt} = y + t - \frac{dy}{dt}$

接下来创建一个函数文件用于 ode 求解器，函数应包含两个输入（通常称为 $t$ 和 $y$），其中 $t$ 是自变量，$y$ 是因变量数组。在这个例子中，$y(1)$ 对应手动推导中的 $y$，$y(2)$ 对应 $z$。包含方程组的函数如下：

function dydt = twoeq(t,y)
    dydt(1) = y(2);
    dydt(2) = y(1) + t - dydt(1);
    dydt = dydt';
end

注意，函数输出已被格式化为列向量，这是 ode 求解器的要求。函数名可以任意，这里 twoeq 具有描述性。

一旦在函数中定义了方程组，就可以将其作为输入用于 ode 求解器。例如，若时间范围定义为 $-1$ 到 $+1$，初始条件定义为 $y = 0$ 和 $z = 0$（即 $y = 0$ 且 $\frac{dy}{dt} = 0$），则命令为：

ode45(@twoeq,[-1,1],[0,0]);

这种已知起始值的问题称为初值问题。

5. 边值问题

初值问题是边值问题的特殊情况。在初值问题中，因变量的起始值已知；而在更一般的边值问题中，因变量在某个不一定是起始点的值已知。

考虑前面描述的两个常微分方程的系统，若不知道 $y$ 和 $\frac{dy}{dt}$ 的初始值，但知道 $t = -1$ 和 $t = 1$ 时 $y$ 的值，这个问题不能用 ode 求解器家族的函数解决，但可以使用 bvp4c 函数。 bvp4c 函数需要三个输入：
- 一个函数句柄，用于待求解的 ode 系统。
- 另一个函数的句柄，用于根据已知边界条件和这些点的预测值求解函数的残差值（误差）。
- 一组初始条件的猜测值。

第一个函数句柄与 ode 求解器使用的完全相同，它应包含感兴趣的导数方程，结果必须是列向量。

为了解决问题，先对所有因变量的初始值进行猜测，然后在求解过程中，程序通过比较计算出的边界值与实际值来检查求解效果。例如，若 $t = -1$ 时 $y = 0$，$t = 1$ 时 $y = 3$，程序将基于 $y$ 和 $\frac{dy}{dt}$ 的初始猜测值求解方程组，然后检查 $t = -1$ 和 $t = +1$ 时的结果是否接近实际值。这通过一个边界条件函数来实现，该函数的方程被安排为：若计算出正确的边界条件，函数值为零。在我们的例子中：

function residual=bc(y_initial, y_final)
    residual(1) = y_initial(1) + 0;
    residual(2) = y_final(1) - 3;
    residual = residual';
end

bc 函数的输入是两个数组（这里称为 y_initial 和 y_final ），每个数组由问题中使用的因变量组成（通常是 $y$，$\frac{dy}{dt}$ 等）。因此， y_initial(1) 是 $y$ 的初始值， y_initial(2) 是 $\frac{dy}{dt}$ 的初始值。若该函数对 y_initial(1)= 0 和 y_final(1)= 3 执行，结果将是一列零。任何其他结果意味着程序计算出的 y_initial 和 y_final 值错误，必须根据函数的算法（有限差分策略）更新初始条件的猜测值。

bvp4c 函数的最后一个输入是问题解的猜测网格，用作求解的起点。MATLAB 提供了一个辅助函数 bvpinit 来帮助创建这个网格，它存储为结构数组。它需要两个输入：一个对应自变量（这里是 $t$）的值数组，以及 ode 方程组中定义的每个变量的初始猜测值。在我们的例子中，有两个方程，所以需要对 $y$ 和 $\frac{dy}{dt}$ 进行猜测。网格不需要特别精细，初始猜测值也不需要非常准确。例如：

initial_guess = bvpinit(-1:.5:1, [0, -1]);

指定了从 $-1$ 到 $1$ 的五个 $t$ 值（$-1$，$-0.5$，$0$，$0.5$，$1$），并在所有 $t$ 值处对 $y$ 和 $\frac{dy}{dt}$ 进行了初始猜测。

一旦创建了描述 ode 系统的函数、定义残差的函数以及使用 bvpinit 创建的初始猜测值，就可以执行 bvp4c 函数：

bvp4c(@twoeq, @bc, initial_guess);

这将返回一个结构数组：

ans =
    x: [1x9 double]
    y: [2x9 double]
    yp: [2x9 double]
    solver: 'bvp4c'

结果是一个结构数组，其中 $x$ 是自变量的值（在这个问题中表示为 $t$），$y$ 值数组对应于 ode 系统的解（这里是 $y$ 和 $\frac{dy}{dt}$）。 bvp4c 函数默认假设自变量的名称是 $x$，但如果将自变量称为 $t$ 也同样适用。默认使用 $x$ 是因为边值问题通常基于位置变化，而不是时间变化。

要访问 $x$ 值数组，只需使用结构语法 ans.x 。若将结果命名为 solution 而不是默认的 ans ，则结构将称为 solution ，$x$ 值将存储在 solution.x 中。最感兴趣的值是 $y$ 值，也可以使用结构语法访问，如 solution.y 。要以类似于 ode 求解器显示的方式绘制结果，使用以下代码：

plot(ans.x,ans.y, '-o');

或者，如果结果命名为 solution ：

plot(solution.x, solution.y, '-o');

6. 偏微分方程

MATLAB 还包括一个有限的偏微分方程求解器 pdepe 。如需更多信息，请查阅 MATLAB 帮助功能：

doc pdepe

7. 插值与曲线拟合

• 插值 ：表格数据可用于总结技术信息，但如果需要表格中未包含的值，必须使用某种插值技术进行近似。MATLAB 中的 interp1 函数可以实现这一点，它需要三个输入：一组 $x$ 值、一组对应的 $y$ 值，以及一组想要估计 $y$ 值的 $x$ 值。该函数默认使用线性插值技术，假设可以将这些中间 $y$ 值近似为 $x$ 的线性函数，即 $y = f(x) = ax + b$。对于每对数据点，会找到一个不同的线性函数，确保近似数据的直线始终通过表格中的点。
  interp1 函数还可以使用高阶近似来建模数据，最常见的是三次样条。近似技术在 interp1 函数的第四个可选字段中以字符串形式指定。如果未指定，函数默认使用线性插值。语法示例如下：

new_y = interp1(tabulated_x, tabulated_y, new_x, 'spline');

除了 interp1 函数，MATLAB 还包括二维插值函数 interp2 、三维插值函数 interp3 和多维插值函数 interpn 。
- 曲线拟合 ：曲线拟合例程与插值技术类似，但它不是连接数据点，而是寻找一个尽可能准确地建模数据的方程。一旦有了方程，就可以计算相应的 $y$ 值。建模的曲线不一定通过测量数据点。MATLAB 的曲线拟合函数是 polyfit ，它使用最小二乘回归技术将数据建模为多项式。该函数返回形如 $y = a_1x^n + a_2x^{n - 1} + a_3x^{n - 2} + \cdots + a_nx + a_{n + 1}$ 的多项式方程的系数。

这些系数可用于在 MATLAB 中创建适当的表达式，也可作为 polyval 函数的输入，以计算任何 $x$ 值对应的 $y$ 值。例如，以下语句找到二阶多项式的系数以拟合输入的 $x - y$ 数据，然后使用第一个语句中确定的多项式计算新的 $y$ 值：

coef = polyfit(x,y,2);
y_first_order_fit = polyval(coef,x);

这两行代码可以通过嵌套函数缩短为一行：

y_first_order_fit = polyval(polyfit(x,y,1),x);

要将数据拟合到通过零点的直线，使用左除来计算模型 $y = a*x$ 的系数：

a = x'\y';

用作输入的 $x$ 和 $y$ 数组必须是列数组。

MATLAB 还包括交互式曲线拟合功能，允许用户不仅使用多项式，还可以使用更复杂的数学函数来建模数据。基本的曲线拟合工具可以从图形窗口的工具菜单中访问，更广泛的工具可在曲线拟合工具箱中找到。

8. 数值微分与积分

数值技术在工程中广泛用于近似导数和积分。求导数和积分的函数也可以在符号工具箱中找到。
- 数值微分 ：MATLAB 的 diff 函数用于求向量中相邻元素值之间的差值。通过将 diff 函数与 $x$ 和 $y$ 值的向量一起使用，可以使用以下命令近似导数：

slope = diff(y)./diff(x);

$x$ 和 $y$ 数据越接近，导数的近似值就越接近真实值。

gradient 函数使用前向差分方法近似数组中第一个点的导数，使用后向差分方法近似数组中最后一个值的导数，使用中心差分方法近似其余点的导数。一般来说，中心差分方法比其他两种技术更能准确地近似导数。
- 数值积分 ：有序数据对的积分使用梯形法则通过 trapz 函数完成。这种方法也可用于函数，通过基于一组 $x$ 值和对应的 $y$ 值创建一组有序对。

函数的积分可以更直接地使用 integral 函数完成。该函数要求用户输入一个函数及其积分限。函数可以表示为匿名函数，也可以表示为存储在单独文件中的函数。 integral 函数尝试返回精确到 $1 * 10^{-6}$ 的答案。

综上所述，MATLAB 提供了丰富的数值技术工具，可用于解决各种微分方程、插值、曲线拟合、数值微分和积分等问题，为工程和科学领域的数值计算提供了强大的支持。

MATLAB 数值求解微分方程全解析

9. 不同微分方程求解器的选择策略

在实际应用中，选择合适的微分方程求解器至关重要，它直接影响到求解的效率和精度。以下是一些选择求解器的策略：
- 非刚性问题 ：
- 若对问题特性了解较少，优先选择 ode45 。它是一种通用的求解器，对于大多数非刚性微分方程能提供较好的精度和效率。例如，在模拟简单的物理系统，如自由落体运动的速度变化时， ode45 可以快速准确地给出结果。
- 如果函数具有 “轻度刚性” 特征，可尝试使用 ode23 。它在某些情况下比 ode45 更高效，尤其是当计算资源有限或者对计算速度有较高要求时。
- 对于需要多步求解的非刚性问题， ode113 是一个不错的选择。它在处理一些复杂的非刚性系统时，可能会比单步求解器更稳定。
- 刚性问题 ：
- 对于刚性微分方程和微分代数方程， ode15s 是常用的求解器。虽然难以预测哪种技术对刚性方程最有效，但 ode15s 在很多情况下都能给出可靠的结果。
- ode23s 适用于刚性问题，特别是当需要使用改进的二阶罗森布罗克公式时。它在处理某些特定类型的刚性系统时可能具有更好的性能。
- ode23t 对于中度刚性微分方程和微分代数方程很有用，当需要无数值阻尼的解时，它能发挥很好的作用。
- ode23tb 也是用于刚性问题的求解器，它使用带梯形法则和二阶向后微分公式的隐式龙格 - 库塔公式，在一些刚性系统中能提供稳定的解。
- 隐式问题 ：对于完全隐式微分方程，使用 ode15i 求解。它专门用于求解形如 $f (y, y′, t ) = 0$ 的隐式微分方程。

问题类型	推荐求解器
非刚性且特性不明	ode45
轻度刚性非刚性	ode23
非刚性多步求解	ode113
刚性及微分代数方程	ode15s
特定刚性问题	ode23s
中度刚性及微分代数方程（无数值阻尼需求）	ode23t
刚性问题	ode23tb
完全隐式微分方程	ode15i

10. 求解微分方程的常见错误及解决方法

在使用 MATLAB 求解微分方程时，可能会遇到一些常见的错误，以下是一些常见错误及相应的解决方法：

graph TD;
    A[输入参数错误] --> B[检查函数句柄是否正确定义];
    A --> C[检查时间跨度和初始条件是否合理];
    D[函数输出格式错误] --> E[确保函数输出为列向量];
    F[求解器不收敛] --> G[尝试更换求解器];
    F --> H[调整初始条件的猜测值];
    I[边界条件不匹配] --> J[检查边界条件函数是否正确设置];
    I --> K[重新评估已知边界条件的准确性];

• 输入参数错误 ：
  • 函数句柄错误 ：确保函数句柄正确指向所需的函数，特别是在使用 M 文件时，要检查文件名和函数名是否一致。
  • 时间跨度和初始条件错误 ：检查时间跨度是否合理，初始条件是否符合实际问题的物理意义。例如，在求解物理系统的运动方程时，初始位置和速度不能为不合理的值。
• 函数输出格式错误 ：对于 ode 求解器，函数输出必须是列向量。如果输出格式不正确，求解器将无法正常工作。在编写函数时，要确保输出格式符合要求。
• 求解器不收敛 ：
  • 更换求解器 ：如果当前求解器无法收敛，可以尝试更换其他适合的求解器。例如，从 ode45 更换为 ode23 或 ode15s 。
  • 调整初始条件 ：初始条件的猜测值可能会影响求解器的收敛性。可以尝试调整初始条件，使其更接近真实解。
• 边界条件不匹配 ：
  • 检查边界条件函数 ：确保边界条件函数正确计算残差值，并且在满足边界条件时函数值为零。
  • 重新评估边界条件 ：检查已知边界条件的准确性，可能是边界条件本身存在错误导致无法求解。

11. 实际应用案例分析

• 物理系统模拟 ：考虑一个简单的弹簧 - 质量 - 阻尼系统，其运动方程可以表示为二阶微分方程。设质量为 $m$，弹簧刚度为 $k$，阻尼系数为 $c$，位移为 $y$，则运动方程为 $m\frac{d^2y}{dt^2} + c\frac{dy}{dt} + ky = 0$。
  首先，将二阶微分方程转化为一阶方程组。令 $z = \frac{dy}{dt}$，则有 $\frac{dy}{dt} = z$，$\frac{dz}{dt} = -\frac{k}{m}y - \frac{c}{m}z$。
  以下是实现该系统求解的 MATLAB 代码：

function dydt = spring_mass_damper(t,y)
    m = 1; % 质量
    k = 10; % 弹簧刚度
    c = 1; % 阻尼系数
    dydt(1) = y(2);
    dydt(2) = -(k/m)*y(1) - (c/m)*y(2);
    dydt = dydt';
end

[t,y] = ode45(@spring_mass_damper,[0, 10], [1, 0]); % 时间跨度从 0 到 10，初始位移为 1，初始速度为 0
plot(t,y(:,1));
xlabel('时间 (s)');
ylabel('位移 (m)');
title('弹簧 - 质量 - 阻尼系统位移随时间变化');

通过这个案例，可以看到如何将实际的物理问题转化为微分方程，并使用 MATLAB 进行求解和可视化。
- 电路分析 ：在一个简单的 RC 电路中，电容的充电和放电过程可以用一阶微分方程描述。设电容为 $C$，电阻为 $R$，电压源为 $V$，电容电压为 $v_c$，则有 $RC\frac{dv_c}{dt} + v_c = V$。
将其转化为标准的一阶微分方程形式 $\frac{dv_c}{dt} = \frac{V - v_c}{RC}$。
以下是求解该电路问题的 MATLAB 代码：

function dvc_dt = rc_circuit(t,vc)
    R = 1000; % 电阻
    C = 0.001; % 电容
    V = 5; % 电压源
    dvc_dt = (V - vc)/(R*C);
end

[t,vc] = ode45(@rc_circuit,[0, 5], 0); % 时间跨度从 0 到 5，初始电容电压为 0
plot(t,vc);
xlabel('时间 (s)');
ylabel('电容电压 (V)');
title('RC 电路电容电压随时间变化');

这个案例展示了如何利用 MATLAB 求解电路中的微分方程，从而分析电路的动态特性。

12. 总结与展望

MATLAB 提供了丰富的数值技术工具，涵盖了微分方程求解、插值、曲线拟合、数值微分和积分等多个方面。通过合理选择求解器和正确设置输入参数，可以高效地解决各种工程和科学领域的数值计算问题。

在未来，随着科学研究和工程应用的不断发展，对数值计算的精度和效率要求也会越来越高。MATLAB 可能会进一步优化其数值求解算法，提供更多更强大的功能。同时，与其他软件和工具的集成也将更加紧密，为用户提供更便捷的计算环境。此外，对于复杂的多物理场问题和大规模系统的求解，可能会出现新的数值技术和求解策略，MATLAB 有望在这些领域发挥更大的作用。

总之，掌握 MATLAB 的数值技术对于从事工程和科学研究的人员来说是非常重要的，它可以帮助我们更好地理解和解决实际问题。