如何在ODE求解器中追踪无导数的显式中间状态？

方案1：循环内同步计算存储（贴合你现有思路） #

你当前的伪代码逻辑完全可行，只需调整状态的存储方式，保证每个时间点的状态都被记录下来：

% 手动划分时间序列
t = linspace(t_start, t_end, step_num);
rates = zeros(length(t), rate_dim);
rates(1, :) = initial_rate_val;
states = zeros(length(t), state_dim);

for i = 1:length(t)-1
    % 计算当前t(i)对应的状态值
    states(i, :) = 2 - rates(i, :); % 替换为你的实际状态计算逻辑
    % 用当前状态求解[t(i), t(i+1)]区间的ODE
    [~, rate_segment] = ode45(@(t, r) dydt(t, r, states(i, :)), [t(i), t(i+1)], rates(i, :));
    % 取区间终点值作为下一时间点的初始rate
    rates(i+1, :) = rate_segment(end, :);
end
% 补全最后一个时间点的状态计算
states(end, :) = 2 - rates(end, :);

plot(t, rates);
plot(t, states);

这种方式逻辑直观，不需要改动核心导数函数，状态与时间点严格对应，不会出现错位问题。

方案2：用ODE求解器的回调函数自动记录 #

如果不想手动控制时间步，多数主流ODE求解器（比如MATLAB的ode45、Python的solve_ivp）支持回调机制，积分过程中每一步都会触发回调，你可以在回调里实时计算并存储状态：
以MATLAB为例：

% 初始化存储变量
states_log = [];
t_log = [];

% 定义回调函数，用于记录状态
function status = log_states(t, r, ~)
    persistent states_log t_log
    % 计算当前状态
    current_state = 2 - r;
    states_log = [states_log; current_state];
    t_log = [t_log; t];
    status = 0; % 让求解器继续积分
end

% 调用ODE求解器并指定回调
[t, rates] = ode45(@(t, r) dydt(t, r), [t_start, t_end], initial_rate_val, ...
    'OutputFcn', @log_states);

% 补全终点状态（防止回调漏捕获）
if t(end) ~= t_log(end)
    states_log = [states_log; 2 - rates(end, :)];
    t_log = [t_log; t(end)];
end

plot(t, rates);
plot(t_log, states_log);

这个方法的优势是求解器会自动选择最优步长，同时能完整记录积分过程中所有时刻的状态，适合对精度要求较高的场景。

方案3：将状态作为“伪状态”加入ODE系统（不推荐） #

如果一定要把状态放进导数函数，可以给状态添加导数为0的方程，让求解器把它当作普通状态变量处理：
导数函数示例：

function dydt = ode_system(t, y)
    % y前半部分是rates，后半部分是states
    rates = y(1:rate_dim);
    states = y(rate_dim+1:end);
    
    % 计算rates的导数
    drdt = ... % 你的实际导数计算逻辑，可直接使用states
    % 状态的导数设为0，保持其值不变
    dstdt = zeros(state_dim, 1);
    
    dydt = [drdt; dstdt];
end

调用时初始值要包含初始状态：

initial_y = [initial_rate_val; initial_state_val];
[t, y] = ode45(@ode_system, [t_start, t_end], initial_y);
rates = y(:, 1:rate_dim);
states = y(:, rate_dim+1:end);

plot(t, rates);
plot(t, states);

这种方法虽然能把状态与rates一起求解，但会增加ODE系统的维度，大型系统下会降低求解效率，且状态导数为0属于冗余计算，仅在特殊场景下考虑使用。

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