Matlab求解ODE得到错误结果：dsolve返回常数值7119/125原因问询

问题根因 #

你得到固定常数解的核心原因是：代入初始值yo后，你定义的yef(t)计算结果为0，导致微分方程右侧的导数项diff(y(t),t)恒等于0，因此y(t)会一直保持初始值不变，不会出现随t变化的项。你得到的7119/125数值恰好等于你定义的初始值yo=0.40680*140，也可以佐证这个结论。

代码存在的错误 #

• 第一类错误（直接导致常数解）：yef(t)的表达式括号嵌套逻辑错误，完全不符合你预期的动力学模型形式。你可以运行以下代码验证：

subs(yef(t),y(t),yo)

运行后会看到返回值为0，证明初始点的导数为0，方程直接达到稳态，自然不会有随t变化的解。你需要核对你原本要写的动力学公式，调整括号的嵌套位置，保证yef(t)在合理的y取值范围内不会恒为0。

• 第二类错误（不导致常数解，但会让结果数值完全错误）：阿伦尼乌斯公式的运算顺序错误，你写的exp(-E1/R*T)按照Matlab运算优先级会先算E1/R再乘以T，实际应该是exp(-E1/(R*T))，漏掉了R*T的括号，会导致反应速率常数K1的计算结果偏差数个数量级。

修正参考 #

你首先需要核对原始动力学方程，修正yef(t)的括号写法，以下为修正通用问题后的示例代码（yef(t)的具体表达式请替换为你自己的正确公式）：

syms y(t) yef(t);
ymax=120*0.40680;
% 示例修正yef表达式括号，此处以饱和动力学为例，请替换为你的正确公式
yef(t)= 0.95/(1 + y(t)/ymax);
yo=0.40680*140;
k10=2.37;
R=0.00831447;
T=395;
HA=0.08;
n1=1.51;
E1=83.3;
% 修正阿伦尼乌斯公式的括号
K1=k10*(10^10)*(HA^n1)*exp(-E1/(R*T));
ode=diff(y(t),t)==-K1*yef(t);
cond=y(0)==yo;
ySol(t)=dsolve(ode,cond)

运行修正后的代码即可得到带t的解析解。

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