在使用 ODE15 求解时，如果使用函数句柄调用向量的索引，如何迭代向量的索引？

Question

我正在尝试求解每次 T 更新时都会更新的向量 k（反应常数），并且 k 的新值用于查找下一个 T 等等。但每次 k 迭代调用 T 时，应使用与 k 的同一索引对应的 T 值，第二次调用 T（即 T_prev）时应引用前一个索引的值。我不太确定我做错了什么。除非必须，我通常不使用函数句柄。

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_0oc, NB_0oc, 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

这是我当前遇到的错误：

无法识别的函数或变量“k”。

Monsanto Plant Disaster>myode 中的错误（第 61 行） xt = NA_0 - NB_0 * (1 - exp(-k(T, T_prev) * t)); % 假设一级动力学

孟山都工厂灾难中的错误（第 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）；

odenumjac 中的错误（第 131 行） Fdel(:,j) = F(Fargs{1:diffvar-1},ydel(:,j),Fargs{diffvar+1:end});

ode15s 中的错误（第 349 行） [dfdy,Joptions.fac,nF] = odenumjac(ode, {t,y,odeArgs{:}}, f0, Joptions); %#好的

Answer 1

您需要将 k 传递到函数中或在函数中定义它。