I have a problem with a program I'm writing for Octave. It is a system of seven differential equations.
Below I show the code:
% 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);
If it would be useful, below I attach the code for the functions mentioned in the main program:
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
I have adapted the code from a problem of a system of 2 differential equations, so maybe it has something to do with my problem, although I don't see it.
The error I get is the following:
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
The error message means you are trying a mathematical operation with two vectors of different size (op1 has size 7x1 op2 has size 5x1). So octave doesn't find the corresponding elements for element 6 and 7 of op1 inside op2 (because op2 ends after element 5).
Try to modify the code to get the same sizes of the vectors. Depending on the mathematical operation you try to do you may additionally need to transpose (rotate) one of the vectors.
Suggested reading: page 16ff of: http://www-mdp.eng.cam.ac.uk/web/CD/engapps/octave/octavetut.pdf