关于Matlab中ode_system函数y(2)、y(3)的ODE求解疑问

我遇到一段用于求解三阶常微分方程 y'''(t) + a y(t) = -b y''(t) + u(t) 的Matlab代码，但对ode_system函数定义中的y(2)、y(3)部分感到困惑，希望了解其含义及对整体求解的作用。代码如下：

% Define the parameters a and b
a = 1;
b = 2;
    
% Define the time horizon [0,1]
time_horizon = [0, 1];
    
% Define the initial conditions for y, y', and y''
initials = [0; 0; 0];
    
% Define the function handle for the input function u(t)
%sin(t) is a common example of a time-varying function.
% You can change the definition of u to any other function of time,
% such as a constant, a step function, or a more complex function, depending on your needs
    
u = @(t) sin(t);
    
% Define the function handle for the system of ODEs
odefunction = @(t, y) ode_system(t, y, a, b, u);
    
% Solve the ODEs using ode45
[t, y] = ode45(odefunction, time_horizon, initials);
    
% Plot the solution
plot(t, y(:,1), '-', 'LineWidth', 2);
xlabel('t');
ylabel('y');
    
function dydt = ode_system(t, y, a, b, u)
    
%Define the system of ODEs
    
dydt = [y(2); y(3); -b*y(3) + u(t) - a*y(1)];
    
end

核心解释：y(2)、y(3)的含义 #

Matlab的ode45这类求解器只能直接求解一阶常微分方程组，因此必须把高阶常微分方程拆解为多个一阶方程的组合，这里的y(2)和y(3)就是拆解过程中引入的状态变量：

• y(1)：对应原方程的未知函数 y(t)
• y(2)：对应y(t)的一阶导数 y'(t)
• y(3)：对应y(t)的二阶导数 y''(t)

拆解过程与作用 #

原方程整理后可以得到三阶导数的表达式：
y'''(t) = -b y''(t) + u(t) - a y(t)

基于状态变量的定义，我们可以构造出三个一阶常微分方程组成的系统：

1. 对y(t)求导：dy(t)/dt = y'(t) = y(2) → 对应dydt的第一个元素
2. 对y'(t)求导：dy'(t)/dt = y''(t) = y(3) → 对应dydt的第二个元素
3. 对y''(t)求导：dy''(t)/dt = y'''(t) = -b*y(3) + u(t) - a*y(1) → 对应dydt的第三个元素

通过这种方式，三阶方程被转化为ode45可处理的一阶方程组，y(2)和y(3)作为中间状态变量，让求解器能够在每一步迭代中跟踪原函数的各阶导数，最终计算出原函数y(t)的数值解。

内容的提问来源于stack exchange，提问作者mikewazowskigreeen