MATLAB: Fixed point(equilibrium points)

equibrium pointsfixed points-

How can I write a script that uses a numerical method to identify any equilibrium values(fixed points) for a system of three ordinary differential equation for a choosen set of parameters.

function Y_m3a=m3a(t,X)

%#ok<*INUSL>

global lwd n a_cca1 b_cca a_cca2 a_prr1 b_prr a_prr2 g_cca g_prr b_toc …

g_toc k_cp k_ct k_tc k_cl k_pl k_pc k_pt

Y=zeros(1,3);

Y(1)=(((a_cca1*b_cca)+(a_cca2*b_cca*((lwd^n)/((k_cl^n) + (lwd^n)))))…

*((k_cp^n)/((k_cp^n) + (X(2)^n)))*((k_ct^n)/((k_ct^n) + (X(3)^n))))…

– (g_cca*X(1));

Y(2)=(((a_prr1*b_prr)+(a_prr2*b_prr*((lwd^n)/((k_pl^n) + (lwd^n)))))…

*((k_pc^n)/((k_pc^n) + (X(3)^n)))*((k_pt^n)/((k_pt^n) + (X(1)^n))))…

– (g_prr*X(2));

Y(3)=(b_toc*((k_tc^n)/((k_tc^n) + (X(1)^n)))) – (g_toc*X(3));

Y_m3a=Y';

end

Implementation of function

%Implementation of the m3 model from Joanito et al. 2018

% Chosen parameter see pdf note for reasons.

global lwd n a_cca1 b_cca a_cca2 a_prr1 b_prr a_prr2 g_cca g_prr b_toc …

g_toc k_cp k_ct k_tc k_cl k_pl k_pc k_pt

lwd = 0.08324;

n = 5 ; % Hill function coefficient.

a_cca1 = 0.2042; % Fraction of cca production rate

a_cca2 = 1-a_cca1; % fraction of cca production rate associated with lwd

a_prr1 = 0.6623; % fraction of prr production rate.

a_prr2 = 1-a_prr1; % fraction of prr production rate associated with lwd

g_cca = 5.5491; % degradation rate of cca.

g_prr = 0.0105; % degradation rate of prr.

g_toc = 1.7349; % degradation rate of toc.

b_cca = g_cca; %production rate of cca

b_prr = g_prr; %production rate of prr

b_toc = g_toc; % production rate of toc

%activators and repressors coefficients used in the Hill function

k_cl = 0.1886; % lwd activating cca genes

k_cp = 0.1261; % prr genes repressing the cca genes.

k_ct = 0.1617; % toc genes repressing the cca genes.

k_pc = 0.0661; % cca genes repressing the prr genes.

k_pl = 0.0264; % lwd activating the prr genes.

k_pt = 0.1327; % toc genes repressing the prr genes.

k_tc = 0.0867; % cca genes repressing the toc genes.

y0=[0.1, 0.1, 0.1];% holds the inital condition

tspan=[0:0.01:20]; % time range with a jump of 0.01

[t,y]=ode45(@m3a,tspan,y0);% odesolver to evaluate the function and return

%t and a matrix y . The vector t will contain the

% time range

plot(t,y)

axis([0 20 0 1])

xlabel('Time')

ylabel('Concentration')

legend({'CCA1/LHY','PRR9/7','PRR5/TOC1'},'Location','east')

fun = @m3a;

x0 = [0,0,0];

Equilibrium_values = fsolve(fun,x0);

error

Not enough input arguments.

Error in m3a (line 9)

*((k_cp^n)/((k_cp^n) + (X(2)^n)))*((k_ct^n)/((k_ct^n) + (X(3)^n))))…

Error in fsolve (line 242)

fuser = feval(funfcn{3},x,varargin{:});

Best Answer

Why would you use an ODE solver to determine equilibrium points for an ODE? Yes, I know it somehow seems logical. But this is not what ODE solvers are designed to do.

An equilibrium point is a point where the function does not change. Just set the derivatives to zero, then solve. That is, an equilibrium point is a point where Y' = 0. The same applies where you have a system of equations. As such, you use fsolve or solve or vpasolve to find that point or points, NOT an ODE solver.

Related Solutions

MATLAB: How to generate chromosome

Is maybe this what are you looking for?

chromo=40;
gene=20;
for i=1:1:chromo
    for j=1:1:gene
    pop=randi([0,1]);
    A(:,j)=[pop]
    end
    B(i,:)=A(:,j)
end
[m,n]=size(B) %just for test

B is matrix with m=40 and n=20, so you have 40 choromosomes with each having 20 genes.

Best of luck.

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

Maybe you mean

function main
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
F_FEED = 100;           %TOTAL FLOWRATE OF THE FEED (mol/s)
T0 = 650;          %INLET TEMPERATURE OF THE REACTOR (K)
P0 = 303975;    %INLET PRESSURE OF THE REACTOR (Pa)
WEIGHT_cat = 100;       %CATALYST WEIGHT (kg)
R = ToBeSpecified;       
Ua = ToBeSpecified;
Feed = zeros(8,1);
Feed(1) = 0.014*F_FEED;           %INLET FLOWRATE OF G (mol/s)
Feed(2) = 0.634*F_FEED;           %INLET FLOWRATE OF H2O (mol/s)
Feed(3) = 0.352*F_FEED;           %INLET FLOWRATE OF N2 (mol/s)
Feed(8) = T0; 
F_G0 = Feed(1);                 
F_H2O0 = Feed(2);              
F_N20 = Feed(3);                
F_A0 = Feed(4);               
F_AC0 = Feed(5);                
F_AA0 = Feed(6);                
F_FA0= Feed(7);                
F_T0 = F_FEED;
[W,F] = ode45(@Reactors , [0 WEIGHT_cat], Feed);
end
function dFdW = Reactors(t,F)
global F_FEED T0 P0 WEIGHT_cat R Ua F_G0 F_H2O0 F_N20 F_A0 F_AC0 F_AA0 F_FA0 F_T0;
dFdW = zeros(8,1);
%dFdW = [r_G]
%       [r_H2O]
%       [r_N2]
%       [r_O2]
%       [r_A]
%       [r_AC]
%       [r_AA]
%       [r_FA]
%       [T]     Temperature
%reaction
Hf_G = -577.9;                    %Heat of formation: G (kJ/mol)
Hf_H2O = -241.818;                %Heat of formation: H2O (kJ/mol)
Hf_N2 = 0;                        %Heat of formation: N2 (kJ/mol)
Hf_A = -65;                       %Heat of formation: A (kJ/mol)
Hf_AC = -370.06;                  %Heat of formation: AC (kJ/mol)
Hf_AA = -166.1;                   %Heat of formation: AA (kJ/mol)
Hf_FA = -108.6;                   %Heat of formation: FA (kJ/mol)
Cp_G = -35.3267443 + 0.645106658*F(8) - 0.00061741*F(8)^2 + 0.000000239229*F(8)^3;               %provide integrated form
Cp_H2O = 32.95266198 - 0.00391977*F(8) + 0.0000222314*F(8)^2 - 0.0000000104953*F(8)^3;
Cp_N2 = 28.883 -0.00157*F(8)+ 0.00000808*F(8)^2 -0.000000002871*F(8)^3;
Cp_A = -14.9617445 + 0.373240705*F(8) - 0.00034559*F(8)^2 + 0.000000129748*F(8)^3;
Cp_AC = -11.7479865 + 0.402868987*F(8) - 0.00031782*F(8)^2 + 0.000000109142*F(8)^3;
Cp_AA = -5.82005966 + 0.23787576*F(8) - 0.00017362*F(8)^2 + 0.0000000543415*F(8)^3;
Cp_FA = 19.10519375 + 0.048680586*F(8) + 0.00000878405*F(8)^2 - 0.0000000150776*F(8)^3;
integraldeltaCp_rxn_1 = 2*(32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)) + (-14.9617445*(F(8)-298) + 0.5*0.373240705*(F(8)^2-298^2) - (1/3)*0.00034559*(F(8)^3-298^3) + (1/4)*0.000000129748*(F(8)^4-298^4)) - ((-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4))); 
integraldeltaCp_rxn_2 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (-11.7479865*(F(8)-298) + 0.5*0.402868987*(F(8)^2-298^2) - (1/3)*0.00031782*(F(8)^3-298^3) + (1/4)*0.000000109142*(F(8)^4-298^4)) - (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4)));
integraldeltaCp_rxn_3 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
integraldeltaCp_rxn_4 = (32.95266198*(F(8)-298) - 0.5*0.00391977*(F(8)^2-298^2) + (1/3)*0.0000222314*(F(8)^3-298^3) - (1/4)*0.0000000104953*(F(8)^4-298^4) + (19.10519375*(F(8)-298) + 0.5*0.048680586*(F(8)^2-298^2) + (1/3)*0.00000878405*(F(8)^3-298^3) - (1/4)*0.0000000150776*(F(8)^4-298^4)) + (-5.82005966*(F(8)-298) + 0.5*0.23787576*(F(8)^2-298^2) - (1/3)*0.00017362*(F(8)^3-298^3) + (1/4)*0.0000000543415*(F(8)^4-298^4)) - (-35.3267443*(F(8)-298) + 0.5*0.645106658*(F(8)^2-298^2) - (1/3)*0.00061741*(F(8)^3-298^3) + (1/4)*0.000000239229*(F(8)^4-298^4)));
Hrxn_1 =  2*Hf_H2O + Hf_A - Hf_G + (integraldeltaCp_rxn_1);               %Heat of reaction: reaction 1
Hrxn_2 =  2*Hf_H2O + Hf_AC - Hf_G + (integraldeltaCp_rxn_2);              %Heat of reaction: reaction 2
Hrxn_3 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_3);        %Heat of reaction: reaction 3
Hrxn_4 =  Hf_H2O + Hf_AA + Hf_FA - Hf_G + (integraldeltaCp_rxn_4);        %Heat of reaction: reaction 4
%INITIAL CONCENTRATIONS  
VOLUMETRICFLOWRATEin = F_FEED*R*T0;                         %Initial volumetric flowrate into the reactor (m^3/s) (assuming ideal gas in)   NEGLECTING P/P0 to start...
C_G0 =      F_G0/VOLUMETRICFLOWRATEin;                      %Inital concentration of G
C_H2O0 =    F_H2O0/VOLUMETRICFLOWRATEin;                    %Inital concentration of H2O
C_N20 =     F_N20/VOLUMETRICFLOWRATEin;                     %Inital concentration of N2
C_A0 =  0;                                                  %Outlet concentration of A

C_AC0 = 0;                                                  %Outlet concentration of AC

C_AA0 = 0;                                                  %Outlet concentration of AA

C_FA0 = 0;                                                  %Outlet concentration of FA

C_T0 =      F_FEED/VOLUMETRICFLOWRATEin;                    %Initial total concentration
%FINAL Concentrations
F_T = F(1) + F(2) + F(3) + F(4) + F(5) + F(6) + F(7);       
C_G =   C_T0*(F(1)/F_T)*(T0/F(8));    %Outlet concentration of G              
C_H2O = C_T0*(F(2)/F_T)*(T0/F(8));      %Outlet concentration of H2O
C_N2 =  C_T0*(F(3)/F_T)*(T0/F(8));      %Outlet concentration of N2              NEGLECTING P/P0 to start...
C_A =   C_T0*(F(4)/F_T)*(T0/F(8));    %Outlet concentration of A
C_AC =  C_T0*(F(5)/F_T)*(T0/F(8));      %Outlet concentration of AC
C_AA =  C_T0*(F(6)/F_T)*(T0/F(8));      %Outlet concentration of AA
C_FA =  C_T0*(F(7)/F_T)*(T0/F(8));      %Outlet concentration of FA
%k-VALUES 
k1 = 23.3244*exp(-45007.1/(R*F(8)));           %Reaction rate constant k1 (m^3/kgcat*s)
k2 = 2528.8*exp(-79244.89/(R*F(8)));           %Reaction rate constant k2 (m^3/kgcat*s)
k3 = 2150.595*exp(-73702.78/(R*F(8)));         %Reaction rate constant k3 (m^3/kgcat*s)
k4 = 0.146768*exp(-24596.97/(R*F(8)));         %Reaction rate constant k4 (m^3/kgcat*s)
%reaction rates
r_1 = -k1*C_G;                                     %Rate of reaction 1
r_2 = -k2*C_G;                                     %Rate of reaction 2
r_3 = -k3*C_G;                                     %Rate of reaction 3
r_4 = -k4*C_G;                                     %Rate of reaction 4
r_G = -(k1+k2+k3+k4)*C_G;                  %Rate of formation of G
r_N2 = 0;                                  %Rate of formation of N2
r_H2O = k1*(C_G)^2+(k3+k4)*G_G;            %Rate of formation of H2O
r_A = k1*C_G;                              %Rate of formation of A
r_AC = k2*C_G;                             %Rate of formation of AC
r_AA = (k3+k4)*C_G;                        %Rate of formation of AA
r_FA = (k3+k4)*C_G;                        %Rate of formation of FA
%Differential Equations to Solve
dFdW(1) = (r_G);
dFdW(2) = (r_H2O);
dFdW(3) = (r_N2);
dFdW(4) = (r_A);
dFdW(5) = (r_AC);
dFdW(6) = (r_AA);
dFdW(7) = (r_FA);
dFdW(8) = (Ua(Ta-T0)+(-r_1)*(-Hrxn_1)+(-r_2)*(-Hrxn_2)+(-r_3)*(-Hrxn_3)+(-r_4)*(-Hrxn_4))/(F(1)*Cp_G + F(2)*Cp_H2O + F(3)*Cp_N2 + F(4)*Cp_A + F(5)*Cp_AC + F(6)*Cp_AA + F(7)*Cp_FA); %dTdW from pg 545 in Fogler;   
end

Best Answer

Related Solutions

MATLAB: How to generate chromosome

MATLAB: What could be the issue if a global variable gets an error “Undefined function or variable ‘F_FEED’ ” when called from a function. The global variable is defined in the file from which the function is called.

Related Question