我为 Octave 编写的程序有问题。它是一个七个微分方程组。
下面我展示代码:
% 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
我已经根据2个微分方程组的问题改编了代码,所以也许它与我的问题有关,尽管我没有看到它。
我得到的错误如下:
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
错误消息意味着您正在尝试使用两个不同大小的向量进行数学运算(op1 的大小为 7x1 op2 的大小为 5x1)。所以octave在op2内部找不到op1的元素6和7的对应元素(因为op2在元素5之后结束)。
尝试修改代码以获得相同大小的向量。根据您尝试执行的数学运算,您可能还需要转置(旋转)其中一个向量。
建议阅读:第 16 页: http://www-mdp.eng.cam.ac.uk/web/CD/engapps/octave/octavetut.pdf