Hi Adrien<br><br>i am new to SCILAB! I just want to say that.<br><br>Yes dx_1 is equal to dx1 but the only reason why i have programmed it like that is becasue the ODE's looks as follows - (see word file attached that can explain the ODE"s better. with regards to the HION and CO it actually refers to [H+] and CO2 as you said. the only reason why n multiplied HION and CO wit hsome ODE's is because some does have a influence on some ODE's.<br>
<br>This is the first time i am working with SCILAB thus i am struggling to understand how SCILAB wants the code so that all 9 ODE's are shown and so that the ODE's that is having a effect on other does happen. I though if you refer to dx(1) for example in a other ODE it means that SCILAb will know the dx(1) has a influence on the other ODE.<br>
<br>The main goal of my assesment is to deliver similar results obtained from MATLAB on SCILAB. all i got was how the grpahs should look like and the ODE's.<br><br><div class="gmail_quote">On Wed, Oct 5, 2011 at 2:43 PM, Adrien Vogt-Schilb <span dir="ltr"><<a href="mailto:vogt@centre-cired.fr">vogt@centre-cired.fr</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div bgcolor="#FFFFFF" text="#330000">
Hi<br>
<br>
Try to isalote your problem<br>
if i understood well, the following code<div class="im"><br>
<br>
<span style="color:rgb(51, 51, 255)">// initial values </span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">x0=[0.1, 1e-5, 0, 15, 1.16,
100,0,297,0.5]';</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">t=0:0.005:400;</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">y=ode(x0, 0, t, f);</span><span style="background-color:rgb(255, 255, 255);color:rgb(51, 51, 255)"><br>
<br>
</span></div>returns y such that sum(y(6:9,:)>x0) == 0 ?<br>
if this is true, we do not need the plots to solve the problem<br>
can you check that ?<br>
<br>
I believe the f function is erroneous.<br>
It seems that dx_1 should be equal to dx(1) at each time step, and
that HION should be equal to x(2) at each time step, etc.<br>
<br>
in other terms, some of your phisical variables seem to be
represented by to variables (i am guessing HION=[H+] and x(2)=[H+]
also) but scilab does not have any chance to know that.<br>
if my guess is right, you have to rewrite the f function in a way
that eliminates all references to HION, dx_1, dx_6 and so on<div><div></div><div class="h5"><br>
<br>
On 05/10/2011 14:27, Jaundre Venter wrote:
<blockquote type="cite">Hi all<br>
<br>
Can someone please explain to me the following:<br>
<br>
I am busy with a project of simulation the production of
penicillin in a bio reactor. Now i have 9 ODE's which i want to
simulate.<br>
<br>
now for some reason the last three graphs i am getting doesn't
show any response what so ever. i am using the following code.<br>
<br>
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(1)=(((mu)*(x(1)))-(((x(1))/(x(6)))*((dx_6)))*(CO)*(HION)),
//biomass concentration X</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(2)=((z*(((mu)*(x(1)))-(((F)*(x(1)))/(x(6)))))+(QQ)),
//hydrogen ion concentration H+</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(3)=((((mupp)*(x(1)))-((K)*(x(3)))-((x(3))/(x(6)))*(dx_6))*(HION)),
//Penicilin concentration P</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(4)=((-((mu)/(Yxs))*(x(1)))-(((mupp)/(Yps))*(x(1)))-((mx)*(x(1)))+((Fsf)/(x(6)))-((x(4)/(x(6)))*(dx_6))),
//Substrate concentration S</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(5)=(-(((mu)/(Yxo))*(x(1)))-(((mupp)/(Ypo))*(x(1)))-(((mo))*(x(1)))+((kla)*(cll-(x(5))))-(((x(5))/(x(6)))*(dx_6))),
//dissolved oxygen</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(6)=((F+Fab+Floss)*(HION)), // culture Volume V</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(7)=(((rq1)*(dx_1)*(x(6)))+(rq2)*(x(1))*(x(6))),
//Heat generation Qrxn</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(8)=((((F)/(sf))*(Tf-(x(8))))+(1/((x(6))*(pcp)))*(QT)),
// Temperature T</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">dx(9)=(((a1)*(dx_1))+((a2)*(x(1)))+(a3)), // CO2
evolution, CO2</span><br style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">
<span style="background-color:rgb(255, 255, 255);color:rgb(51, 102, 255)">endfunction</span><br style="background-color:rgb(255, 255, 255)">
<br>
now when i ask for plotting the graphs i am using the following.:<br>
<br>
<span style="color:rgb(51, 51, 255)">// initial values </span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">x0=[0.1, 1e-5, 0, 15, 1.16,
100,0,297,0.5]';</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">t=0:0.005:400;</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">y=ode(x0, 0, t, f);</span><span style="background-color:rgb(255, 255, 255);color:rgb(51, 51, 255)"></span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">// the plots of each
variable</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="BIOMASS
CONCENTRATION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.x_label.text="Time,
hours";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="X,g/l ";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(1);clf; //Opens and
clears figure 1</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(1,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="HYDROGEN ION
H+ CONCENTRATION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="H+,mol/l
";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(2);clf; //Opens and
clears figure 2</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(2,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="PENICILLIN
CONCENTRATION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="P,g/l ";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(3);clf; //Opens and
clears figure 3</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(3,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="SUBSTRATE
CONCENTRATION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="S,g/l ";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(4);clf; //Opens and
clears figure 4</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(4,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="DISSOLVED
OXYGEN CONCENTRATION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="C_l,g/l ";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(5);clf; //Opens and
clears figure 5</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(5,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="CULTURE
VOLUME"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="V,l";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(6);clf; //Opens and
clears figure 6</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(6,:))</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="HEAT OF
REACTION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="Qrxn,cal";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(7);</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">clf; //Opens and clears
figure 7</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(7),:)</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="TEMPERATURE"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="T,Kelvin";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(8);</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">clf; //Opens and clears
figure 8</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(8),:)</span><br style="color:rgb(51, 51, 255)">
<br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.title.text="CO2
EVOLUTION"</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">da.y_label.text="CO2,mmol/l/";</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">scf(9);</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">clf; //Opens and clears
figure 9</span><br style="color:rgb(51, 51, 255)">
<span style="color:rgb(51, 51, 255)">plot(t,y(9),:)</span><br>
<br>
Am i doing something wrong? before the ODE's i have just
programmed the initial values and constants :<br>
<br>
f<span style="color:rgb(51, 102, 255)">uncprot(0);</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">function dx = f(t,x)</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">K1=1.0e-10 //mol/l</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">K2=7.0e-05 //mol/l</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Kx=0.15 //
Contois saturation constant, g/l</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Kox=2e-02 //
oxygen limitation constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mux=0.092 //
maitenance coefficient on subsrate</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">p=3
//constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Kp=0.0002 //
inhibition constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Kop=2e-02 // oxygen
limitation constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">K=0.04 // Penicillin
hydrolysis constant, per h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Yxs=0.45 // yield
constant,g biomass/g glucose = dimensionless</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Yps=0.90 // yield
constant, g pinicillin/ g glucose = dimensionless</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mx= 0.014 //
Maintenance coefficient on substrate, per h </span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Yxo=0.04 // yield
constant, g biomass/g oxygen = dimensionless</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Ypo=0.20 // yield
constant, g penicillin/g oxygen= dimensionless</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mo= 0.467 //
maintenance coefficient of oxygen, per h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mup=0.0005 // specific
rate of penicilline production (per h)</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">sf=600 // Feed substrate
concentration, g/l</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">kla=23 // function of
agitation power input and oxugen flow rate, dimensional</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">cll=1.16 // dissolved
oxygen concentration, g/l</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Cab=3 //
concentrations in both solutions</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Fa=5 // acid flow
rate, l/h !! </span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Fb=5 // base flow
rate, l/h !! </span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">delta_t=0.01 // time
step in digital PID controller - arbitrary value!!!</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">z=10e-5 // constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">F=0.042 // feed
substrate flow rate l/h </span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">T0=273 //
temperature at freezing, K</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Tv=373 //
temperature at boiling</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">T=298 // feed temp of
substrate</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">h=(2.5e-4) //
constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Floss=(x(6)*(h)*(exp(5)*((T-T0)/(Tv-T0))))</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Fab=Fa+Fb // volume
increase due to influx of acid Fa and base Fb</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Fsf=((sf)*(F))</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">kg= 7e-3 // Arrhenius
constant for growth</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">kd=10e33 // Arrhenius
constant for cell death</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Eg= 5100 // Activation
energy for growth, cal/mol</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Ed= 50000 // Activation
energy for cell death, cal/mol</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">R= 1.987 // gas
constant, cal/mol k</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">T= 297 // Temperature</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">RT= R*T</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">alpa= 70 // constant in
Kla</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">betha= 0.4 // constant
in Kla</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Pw= 30 // Agitation
power input, W</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">fg= 8.6 // Flow rate of
oxygen</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">V=100 // Volume</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">QE=
((kg*exp(-(Eg/RT)))-(kd*exp(-(Ed/RT))))</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">kla=
alpa*((sqrt(fg)*(Pw/x(6)))^betha)</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mu
=(((mux)/(1+((K1)/(x(2)))+((x(2))/(K2))))*((x(3))/(((Kx)*(x(1)))+(x(3))))*((x(5))/(((Kox)*(x(1)))+(x(5))))*(QE))
// Specific growth rate</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">mupp =
((mup)*((x(4))/((Kp)+(x(4))+(x(4)^2)/(K1)))*((x(5)^p)/((Kop)*(x(1)))+(x(5)^p)))
// Specific penicillin production rate</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">B
=(((1e-14/x(2)-x(2))*x(6)-Cab*(Fa+Fb)*delta_t)/(x(6)+(Fa+Fb)*delta_t))</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">QQ
=((-B+sqrt(B^2+4e-14))/2-(x(2)))*(1/delta_t)</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">dx_6 = (F+Fab+Floss)
//Culture Volume V</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">dx_1 =
(((mu)*(x(1)))-((x(1))/(x(6)))*(dx_6)) //biomass concentration X</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">rq1 = 60 // yield of
heat generation, cal/g biomass</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">rq2 = 1.6783e-4 //
Constant, cal/g biomass h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Tf = 296 // substrate
feed temperature, Kelvin</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">a = 1000 // heat
transfer coefficient of cooling/heating liquid, cal/h degree C</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">b = 0.60 // constant</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">Fc=0.1 // Cooling water
flow rate, not sure about value, l/h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">pcCpc = 1/2000 //
Density times heat capacity of cooling liquid, per l degree C</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">pcp = 1/1500 // density
times heat capacity of medium</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">QT =
((x(7)-(((a)*(Fc^b+1))/((Fc)+((a)*(Fc^b))/2*pcCpc))))</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">a1=0.143 // constant
relating CO2 to growth, mmol CO2/g biomass</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">a2=4e-7 // Constant
relating CO2 to mainteneance energy, mmol CO2/g biomass h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">a3=1e-4 // Constant
relating CO2 to penicillin production, mmol CO2/l h</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">CO=
(((a1)*(dx_1))+((a2)*(x(1)))+(a3)), // CO2 evolution, CO2</span><br style="color:rgb(51, 102, 255)">
<span style="color:rgb(51, 102, 255)">HION=((z*(((mu)*(x(1)))-(((F)*(x(1)))/(x(6)))))+(QQ))</span><br>
<br>
Thanks.<br>
</blockquote>
<br>
<br>
</div></div><font color="#888888"><div>-- <br>
Adrien Vogt-Schilb (Cired) <br>
Tel: (+33) 1 43 94 <b>73 77</b></div>
</font></div>
</blockquote></div><br>