GNU Octave求解含自定义R(T)的ODE代码错误原因咨询

你的Octave求解代码存在两个核心致命错误，会直接导致求解失败或结果完全错误：

• 导数函数参数顺序不符合lsode接口要求
  Octave内置的lsode求解器规定，用户传入的导数计算函数必须遵循xdot = f(状态变量, 时间变量)的签名格式，也就是第一个入参是当前步的状态值（本问题中就是温度T），第二个入参才是当前时间t。你写的函数签名是function xdot = f(t, T)，把时间和状态的位置写反了，求解器运行时会把温度值传给参数t、把时间值传给参数T，所有计算逻辑都会完全错乱。

  注意：这个接口规则和Matlab的ode系列求解器（如ode45）正好相反，Matlab要求时间在前、状态在后，是很多用户容易混淆的坑点。

• 错误对标量状态值做数组索引操作
  在求解器每一步调用导数函数计算时，传入的状态参数T是当前求解步对应的温度标量值，不是存储了全时间序列温度结果的数组。你代码里写的T(t)是试图用浮点数格式的时间值t作为数组下标去索引T，一方面T本身不是数组，另一方面浮点数无法作为合法的数组下标，运行时会直接触发索引报错。计算时直接使用T本身即可，不需要加(t)下标。

修正后的可运行代码 #

1;

function xdot = f(T, t)
  % 生热项
  heat = 987 * 0.0000696 * (1 + 0.0038 * (T - 25));
  % 对流散热项
  conv_loss = 0.0168 * (T - 25);
  % 辐射散热项
  rad_loss = 3.25e-13 * (T^4 - 25^4);
  xdot = heat - conv_loss - rad_loss;
endfunction

% 生成时间序列
t = linspace(0, 3600, 1000)';
% 调用lsode求解，初始温度100
[T_res, istate, msg] = lsode("f", 100, t);

补充说明：你原代码里把t=linspace(...)直接写在lsode入参位置的写法语法上可以运行，但可读性很差，建议单独赋值后传入。

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