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>