使用ODE15s时，如何通过函数句柄迭代向量索引求解反应常数k？

我想求解随温度T更新的反应常数向量k，每次迭代k时，首次调用的T用当前索引值，第二次调用的T_prev用前一索引值。但运行代码时报错Unrecognized function or variable 'k'，具体代码和错误信息如下：

NA_0 = 3170; % mols of ONCB
NB_0 = 43000; % mols of NH3
B = NB_0 / NA_0;
dHrx = -5.9E5;
Vn = 3.265; % Volume of reactant ONCB at normal operating conditions
T_a = 298; % ambient temperature (K)

C_pA = 40E-3; % heat capacity of ONCB
C_pB = 8.38E-3; % heat capacity of H2O
C_pW = 18E-3; % heat capacity of NH3

UA = 35.85;
xt = linspace(0, 120, 240)'; % time vector

% Initial rate constant k at T_a
k = 1.167E-4;

tspan = [0 120]; 
ic = 448; % initial temperature in K

% ODE function for temperature change
[t_s, T] = ode15s(@(t_s, T) myode(t_s, T, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA), tspan, ic);

% Plot the results
plot(t_s, T)
title("Excess Reagents")
xlabel('Time (minutes)')
ylabel('Temperature (K)')
title('Temperature vs Time')


function dTdt = myode(t, T, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA)
    % Calculate the current rate constant k based on temperature
    
    if T == ic
        T_prev = ic - 0.1;
        k = @(T, T_prev) 1.167E-4 * exp(11273 / (1.987 * ((1/T) - (1 ./ T_prev))));
    else
        for z = 2:length(T)
            T_prev(z) = T(z-1);
            k = @(T, T_prev) 1.167E-4 * exp(11273 / (1.987 * ((1/T) - (1 ./ T_prev))));
        end
    end
  
    % Calculate dX/dt
    xt = NA_0 - NB_0 * (1 - exp(-k(T, T_prev) * t)); % Assuming first-order kinetics
    dXdt = k(T, T_prev) * NA_0 .* (1 - xt) .* (B - 2 .* xt) ./ Vn;
    
    % Calculate the heat generated by the reaction
    Q_b = k(T, T_prev) * (NA_0^2) .* (1 - xt) .* (B - 2 .* xt) .* (-dHrx) / Vn;

    % Calculate the heat removed by the reactor cooling
    Q_r = UA * (T - T_a);

    % Calculate the total heat capacity
    NC_p = NA_0 * C_pA + NB_0 * C_pB + NB_0 * C_pW;

    % Check if any element of Q_b is infinite
    if any(isinf(Q_b))
        dTdt = 0;
    else
        % Calculate the rate of temperature change
        dTdt = (Q_b + Q_r) / NC_p; % Heat generated - Heat removed
    end
end

错误信息 #

Unrecognized function or variable 'k'.

Error in Monsanto Plant Disaster>myode (line 61)
xt = NA_0 - NB_0 * (1 - exp(-k(T, T_prev) * t)); % Assuming first-order kinetics

Error in Monsanto Plant Disaster (line 24) [t_s, T] = ode15s(@(t_s, T)
myode(t_s, T, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA),
tspan, ic);

Error in odenumjac (line 131)
Fdel(:,j) = F(Fargs{1:diffvar-1},ydel(:,j),Fargs{diffvar+1:end});

Error in ode15s (line 349)
[dfdy,Joptions.fac,nF] = odenumjac(ode, {t,y,odeArgs{:}}, f0, Joptions); %#ok

核心问题 #

1. 作用域与变量定义漏洞：主脚本里的k、ic无法被myode函数直接访问；函数内的k仅在if/else块定义，但ode15s调用myode时传入的T是单个标量，else块的for循环永远不会执行（length(T)=1），导致k未定义。
2. ODE逻辑误解：你试图用数组索引获取T_prev，但ODE求解器每次只传入当前时刻的温度值，不是整个温度数组。
3. 不必要的函数句柄：计算k不需要用函数句柄，直接计算标量值即可。

修复后的代码 #

NA_0 = 3170; % mols of ONCB
NB_0 = 43000; % mols of NH3
B = NB_0 / NA_0;
dHrx = -5.9E5;
Vn = 3.265; % Volume of reactant ONCB at normal operating conditions
T_a = 298; % ambient temperature (K)

C_pA = 40E-3; % heat capacity of ONCB
C_pB = 8.38E-3; % heat capacity of H2O
C_pW = 18E-3; % heat capacity of NH3

UA = 35.85;
xt = linspace(0, 120, 240)'; % time vector

% Initial rate constant k at T_a
k_initial = 1.167E-4;
T_initial = 448; % initial temperature in K

tspan = [0 120]; 
% 扩展状态变量：包含当前温度T和前一时刻温度T_prev
ic = [T_initial; T_initial - 0.1]; 

% ODE function for temperature change
[t_s, sol] = ode15s(@(t, y) myode(t, y, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA, k_initial), tspan, ic);

% 提取温度数据
T = sol(:,1);

% Plot the results
plot(t_s, T)
xlabel('Time (minutes)')
ylabel('Temperature (K)')
title('Temperature vs Time')


function dydt = myode(t, y, NA_0, NB_0, B, dHrx, Vn, T_a, C_pA, C_pB, C_pW, UA, k_initial)
    T_current = y(1);
    T_prev = y(2);
    
    % 直接计算当前k值，无需函数句柄
    k = k_initial * exp(11273 / (1.987 * ((1/T_current) - (1/T_prev))));
    
    % 修正转化率计算：假设反应为A + 2B -> 产物，一级动力学针对A
    xt = 1 - exp(-k * t); % 转化率xt = (NA0 - NA)/NA0
    NA = NA_0 * (1 - xt);
    NB = NB_0 - 2 * NA_0 * xt;
    
    % 修正反应速率计算
    dXdt = k * NA .* NB ./ (NA_0 * Vn);
    
    % 计算反应产热
    Q_b = k * NA .* NB .* (-dHrx) / Vn;

    % 计算冷却移除的热量
    Q_r = UA * (T_current - T_a);

    % 计算总热容（包含生成的水）
    NC_p = NA * C_pA + NB * C_pB + (2 * NA_0 * xt) * C_pW;
    
    % 处理无穷值情况
    if isinf(Q_b)
        dTdt = 0;
    else
        % 修正热量平衡：产热减去移除的热量
        dTdt = (Q_b - Q_r) / NC_p;
    end
    
    % 返回状态变量导数：当前温度变化率，以及前一温度的更新（下一时刻T_prev=T_current）
    dydt = [dTdt; dTdt];
end

关键修复点 #

1. 扩展状态向量：把T_prev加入ODE状态变量，确保每次迭代能直接获取前一时刻温度。
2. 移除函数句柄：直接计算k的标量值，彻底解决变量未定义问题。
3. 修正动力学与热量公式：原代码中转化率、反应速率、热量平衡的逻辑存在错误，重新推导后修正。
4. 参数传递修正：将初始k等必要参数传入函数，解决作用域问题。

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