使用 MATLAB 的 ode45 求解常微分方程

使用 MATLAB 的 `ode45` 求解常微分方程

简介

ode45 是 MATLAB 中一个常用的数值解常微分方程 (ODE) 的函数。它基于经典的 Runge-Kutta 方法，适用于求解非刚性方程。本文将介绍 ode45 的基本使用方法和一些简单的实例。

非刚性方程是指那些在数值求解时，不需要特别小的时间步长来保证稳定性的常微分方程。这类方程在数值求解过程中一般没有明显的数值振荡和稳定性问题，可以使用常规的数值方法如 ode45 等进行求解。

“45” 代表的是该方法使用了一个四阶的 Runge-Kutta 方法和一个五阶的 Runge-Kutta 方法。

四阶和五阶 Runge-Kutta 方法

四阶 Runge-Kutta 方法 (RK4) ：这个方法以其简单性和广泛的精度而著名。它在每个时间步长 h 中，通过评价四次斜率，来估计当前点的增量。
五阶 Runge-Kutta 方法 (RK5) ：它是四阶方法的扩展，通过评价额外次数的斜率，进一步提高了精度。

ode45 的基本原理

ode45 组合使用了四阶和五阶 Runge-Kutta 方法，属于所谓的 Runge-Kutta-Fehlberg 方法。其基本工作原理如下：
每一步计算时，ode45 同时执行四阶和五阶的方法。
比较两者的结果，估计误差。
根据误差估计，动态调整步长：如果误差较大，减小步长；如果误差较小，增大步长。

`ode45` 的基本语法

在 MATLAB 中，调用 ode45 的基本语法如下：

matlab
[t, y] = ode45(odefun, tspan, y0, options)

其中：

odefun：定义微分方程的函数句柄。
tspan：时间区间向量，表示求解的时间范围。
y0：初始条件向量，在起始时间点的状态。
options：可选参数，设置求解的选项（如误差容限）。

函数返回：

t：时间点数组。
y：在这些时间点上求解得到的解。

定义微分方程

必须定义一个函数，该函数返回状态变量的导数。典型地，这个函数接收当前时间 t 和状态 y 作为输入，返回导数 dy/dt。

例如，对简单的微分方程 $\dot{y} = -2y$ ：

matlab
function dydt = odefun(t, y)
    dydt = -2 * y;
end

你可以将这个函数保存为一个 .m 文件，文件名为 odefun.m。

使用 `ode45` 求解

使用 ode45 求解上述微分方程，初始条件为 $y(0) = 1$ 和时间范围为 $[0, 5]$ ：

matlab
% 定义初始条件和时间范围
y0 = 1;
tspan = [0 5];

% 调用 ode45 求解
[t, y] = ode45(@odefun, tspan, y0);

% 绘制结果
plot(t, y);
xlabel('Time');
ylabel('Solution y');
title('Solution of dy/dt = -2y');
grid on;

示例：二阶微分方程

假设我们要求解一个二阶微分方程，比如： $\ddot{x} + 2\dot{x} + x = 0$

我们可以将其转化为一阶微分方程组：

\begin{cases} \dot{x_1} = x_2 \\ \dot{x_2} = -2x_2 - x_1 \end{cases}

分别用 $x_1$ 和 $x_2$ 来表示变量 $x$ 和 $\dot{x}$ ，然后定义状态变量向量 $y = [x_1; x_2]$ 。

定义函数：

matlab
function dydt = second_order_ode(t, y)
    dydt = zeros(2, 1);
    dydt(1) = y(2);
    dydt(2) = -2 * y(2) - y(1);
end

使用初始条件 $x(0) = 1$ 和 $\dot{x}(0) = 0$ ：

matlab
% 定义初始条件和时间范围
y0 = [1; 0];
tspan = [0 10];

% 调用 ode45 求解
[t, y] = ode45(@second_order_ode, tspan, y0);

% 提取解
x = y(:, 1);
x_dot = y(:, 2);

% 绘制结果
figure;
subplot(2,1,1);
plot(t, x);
xlabel('Time');
ylabel('x');
title('Solution of \ddot{x} + 2\dot{x} + x = 0');
grid on;

subplot(2,1,2);
plot(t, x_dot);
xlabel('Time');
ylabel('\dot{x}');
grid on;

设置选项参数

可以通过 odeset 定制求解选项，如改变容差：

matlab
options = odeset('RelTol', 1e-4, 'AbsTol', 1e-5);
[t, y] = ode45(@odefun, tspan, y0, options);

这可以控制求解的精度和性能，具体取决于问题的需求。

总结

MATLAB 的 ode45 提供了一种强大且简便的方法来数值解常微分方程。通过定义合适的微分方程函数（odefun），设定初始条件和时间范围，你可以在 MATLAB 中轻松求解各种 ODE 问题。

目录