Octave求解7阶微分方程组时出现维度不匹配错误求助

错误信息 #

error: operator +: nonconformant arguments (op1 is 7x1, op2 is 5x1)
error: called from
    starting_stepsize at line 66 column 6
    ode23 at line 200 column 25
    main_vinos at line 59 column 12

主程序代码 #

% Limpiar pantalla
close all;
clear all;

% Declaración de variables que usaremos al definir el valor de las derivadas
global mu_X0 Y_EA mu_S0 mu_L mu_DT mu_E0 mu_SD0 k_E k_S k_X mu_DY mu_AB S_0

% Concentración de tipos de bacterias; recordar que el porcentaje de bacterias
% vivas, latentes y muertas debe sumar 100%
X_inic=4;          % Concentración inicial de bacterias en g/L
X_A0=0.02.*X_inic; % Concentración inicial de bacterias activas en g/L
X_L0=0.48.*X_inic; % Concentración inicial de bacterias latentes en g/L
X_D0=0.5.*X_inic;  % Concentración inicial de bacterias muertas en g/L

% Otras condiciones iniciales de las variables de estado del proceso
S_0=130;           % Concentración inicial de sustrato (azúcares) en g/L
EtOH_0=0;          % Concentración inicial de etanol en g/L
DY_0=0;            % Concentración inicial de diacetilo en mg/L
EA_0=0;            % Concentración inicial de acetato de etilo en mg/L

T=13+273.15;       % Temperatura en °C + 273.15 (o sea Kelvin)

% Valores de parámetros calculados, ajustados a datos experimentales en función
% de la temperatura
mu_X0=exp(108.31-31934.0./T);                % Tasa máxima de crecimiento
Y_EA=exp(89.92-26589./T);                    % Coeficiente estequiomético de proporción EA/S
mu_S0=exp(-41.92+11654.64./T);               % Tasa específica máxima de consumo de S
mu_L=exp(30.72-9501.54./T);                  % Tasa específica de activación celular
mu_DT=exp(130.16-38313./T);                  % Tasa de muerte celular
mu_E0=exp(3.27-1267.24./T);                  % Tasa máxima de producción de EtOH_0
mu_SD0=exp(33.82-10033.28./T);               % Tasa específica de sedimentación máxima
k_E=exp(-119.63+34203.95./T);                % Función de afinidad de ecuación M-M en producción de EtOH
k_S=k_E;                                     % Constante de afinidad para consumo de S (M-M)
k_X=0.5;                                     % Constante de crecimiento celular
%mu_DY=-6.1344E-8.*T.^2+8.4266E-6.*T-1.7672E-2; % Tasa de producción de diacetilo
mu_DY=exp(1853.5-538541.2./T);               % Tasa de producción de diacetilo
%mu_AB=-9.1384E-7*T^2+6.7071E-5*T-0.1251E-3; % Tasa efectiva de reducción de diacetilo
mu_AB=exp(1304.2-379083.1./T);               % Tasa efectiva de reducción de diacetilo

% Límites de tiempo
t0=0;       % Tiempo inicial
tf=16;     % Tiempo final
Dtplot=0.1; % Intervalo de tiempo

% Condiciones iniciales
X_A=X_A0;
X_L=X_L0;
X_D=X_D0;
S=S_0;
EtOH=EtOH_0;
DY=DY_0;
EA=EA_0;
x_inicial=[X_A;X_L;X_D;S;EtOH;DY;EA];

% Resolución de ecuación diferencial
options=odeset('RelTol',1E-5,'AbsTol',1E-5);
options=odeset(options,'Stats','on','Events',@events);
t=[t0:Dtplot:tf];
[tout,xout]=ode45(@vinos_odes,t,x_inicial,options);

函数代码 #

function dxdt=vinos_odes(x,t)

  % Declaración de variables que usaremos al definir el valor de las derivadas
  global mu_x0 Y_EA mu_S0 mu_L mu_DT mu_E0 mu_SD0 k_E k_S k_X mu_DY mu_AB S_0

  % Asignación del valor de cada variable dependiente al vector x, que contienen
  % todas las variables dependientes cuya derivada se toma en las ecuaciones
  % diferenciales
  x=zeros(7,1);
  X_A=x(1,1);   % Células activas en g/L
  X_L=x(2,1);   % Células latentes en g/L
  X_D=x(3,1);   % Células muertas en g/L
  S=x(4,1);     % Concentración de azúcares en g/L
  EtOH=x(5,1);  % Concentración de etanol en g/L
  DY=x(6,1);    % Concentración de diacetilo en mg/L
  EA=x(7,1);    % Concentración de acetato de etilo en mg/L

  % Funciones auxiliares del proceso de fermentación
  mu_X=(mu_x0*S)./(k_X+EtOH);                  % Tasa específica de crecimiento de células activas
  mu_SD=(0.5.*mu_SD0.*S_0)./(0.5.*S_0+EtOH);   % Tasa de sedimentación de células muertas
  mu_S=(mu_S0.*S)./(k_S+S);                    % Tasa específica de consumo de sustrato
  mu_E=(mu_E0.*S)./(k_E+S);                    % Tasa de producción de etanol
  f=1-(EtOH./(0.5.*S_0));                      % Tasa de inhibición para producción de etanol

  % Ecuaciones dinámicas asociadas a la producción de cerveza
  dX_Adt=X_A.*(mu_X-mu_DT)+mu_L.*X_L;      % Células activas
  dX_Ldt=-mu_L.*X_L;                       % Células latentes
  dX_Ddt=-mu_SD.*X_D+mu_DT.*X_A;           % Células muertas
  dSdt=-mu_S.*X_A;                         % Sustrato fermentable
  dEtOHdt=f*mu_E.*X_A;                     % Concentración de etanol
  dDYdt=mu_DY.*S.*X_A-mu_AB.*DY.*EtOH;     % Concentración de diacetilo
  dEAdt=Y_EA.*mu_X.*X_A;                   % Concentración de acetato de etilo

  dxdt=[dX_Adt;dX_Ldt;dX_Ddt;dSdt;dEtOHdt;dDYdt;dEAdt];
endfunction

function [value, isterminal, direction]=events(t,x)
  X_A=x(1,1);   % Células activas en g/L
  X_L=x(2,1);   % Células latentes en g/L
  X_D=x(3,1);   % Células muertas en g/L
  S=x(4,1);     % Concentración de azúcares en g/L
  EtOH=x(5,1);  % Concentración de etanol en g/L
  DY=x(6,1);    % Concentración de diacetilo en mg/L
  EA=x(7,1);    % Concentración de acetato de etilo en mg/L

  value=S;
  isterminal=1;
  direction=-1;
endfunction

错误原因分析 #

• 覆盖输入状态向量：在vinos_odes函数中，x=zeros(7,1);语句直接将ODE求解器传入的7维状态向量覆盖为全零向量，导致后续计算使用的都是无效的零值。
• 变量名大小写不匹配：主程序中定义的参数是mu_X0（大写X），但函数中引用的是mu_x0（小写x），Octave变量名大小写敏感，导致该参数未被正确加载，后续计算mu_X时出现异常，最终使得dxdt的实际维度与预期的7维不符，引发维度不匹配错误。

修复方案 #

1. 删除覆盖输入向量的代码：移除vinos_odes函数中的x=zeros(7,1);语句，保留ODE求解器传入的原始状态向量x。
2. 修正变量名拼写：将函数中的mu_x0统一改为mu_X0，与主程序中的参数定义保持一致。

修改后的vinos_odes函数核心片段 #

function dxdt=vinos_odes(x,t)

  % Declaración de variables que usaremos al definir el valor de las derivadas
  global mu_X0 Y_EA mu_S0 mu_L mu_DT mu_E0 mu_SD0 k_E k_S k_X mu_DY mu_AB S_0

  % Asignación del valor de cada variable dependiente al vector x
  X_A=x(1,1);   % Células activas en g/L
  X_L=x(2,1);   % Células latentes en g/L
  X_D=x(3,1);   % Células muertas en g/L
  S=x(4,1);     % Concentración de azúcares en g/L
  EtOH=x(5,1);  % Concentración de etanol en g/L
  DY=x(6,1);    % Concentración de diacetilo en mg/L
  EA=x(7,1);    % Concentración de acetato de etilo en mg/L

  % Funciones auxiliares del proceso de fermentación
  mu_X=(mu_X0*S)./(k_X+EtOH);                  % Tasa específica de crecimiento de células activas
  % 其余代码保持不变...
endfunction

额外优化建议 #

• 替换全局变量：将参数打包为结构体，作为额外参数传入ODE函数，避免变量名冲突，提升代码可维护性。示例：

  % 主程序中打包参数
  params.mu_X0 = mu_X0;
  params.Y_EA = Y_EA;
  % ...其余参数
  [tout,xout]=ode45(@(t,x) vinos_odes(t,x,params),t,x_inicial,options);
  

  函数定义修改为function dxdt=vinos_odes(t,x,params)，通过params.mu_X0访问参数。
• 维度调试：在函数末尾添加disp(size(dxdt))，快速验证输出向量的维度是否符合预期。

内容的提问来源于stack exchange，提问作者Ben Duncan Málaga Espichán