[SciPy-User] weired results with ode solver

Oz Nahum nahumoz at gmail.com
Sat Feb 6 13:33:33 EST 2010


Hi Warren,
Thanks for your answer,

here is my updated code:
#% This file is a matlab script to simulate a batch problem as described
#% in problem 1a.

from scipy.integrate import odeint
from numpy import *
#%clear all;
#clc;

#%Define variables
Coi=3              # initial concentration of oxigen [mg/l]
Csi=10             # initial concentration of substrat [mg/l]
Cbi=0.2            # initial concentration of biomass [mg/l]

#define vector of concentrations
#c1 = array((Coi,Csi,Cbi))
c1=(Coi,Csi,Cbi)
y0=c1
ko=0.1             #% Monod coefficients Oxigen [mg/l]
ks=0.1             #% Monod coefficients Substract [mg/l]

umax=0.125/86400 #  % Maximun growth rate [1/d]

Yo=0.125#;           % Oxigen Yield [mg X / mg O]
Ys=0.25#;            % Substract Yield [mg X / mg S]

Kdec=0.01/86400#;    % Decay rate
kb=1#;              % Biomass Inhibition constant [mg/L]

def rate_func(c,time,ko,ks,umax,Kdec,Yo,Ys,kb):
	co=c[0]#; %concentration of Oxygen
	cs=c[1]#; %concentration of Substrate
	X=c[2]#; %concentration of bacteria
	kgr=umax*(co/(ko+co))*(cs/(ks+cs))#; % specific uptake rate [mol/m3/s]
	Ib=1+X/kb#;  				 		% inhibiton caused by biomass
	foxi=((-kgr/Yo/Ib)*X)#;  %dc/dt of Oxygen
	fsub=((-kgr/Ys/Ib)*X)#;  %dc/dt of Substrate
	fbio=((kgr/Ib-Kdec)*X)#; %inhibition for growth only
	f=array((foxi,fsub,fbio))
	return  f

#Solve ODE
timerange=arange(1,43*86400,3600)
c=odeint(rate_func,c1,timerange,args=(ko,ks,umax,Kdec,Yo,Ys,kb))

from pylab import *
plot(timerange/86400, c[:,0], 'bo')
plot(timerange/86400, c[:,1], 'go')
plot(timerange/86400, c[:,2], 'ko')
show()


I seems that it runs much faster than matlab (though I have not really
checked). It's syntacticly clearer to read. Much more fun.
There is one issue that bothers me still. I don't want to show the
concentrations every second. So I divide by 86400 sec perday. However,
this results in smeared concentrations, where every day has  a few of
them. Why is that ? how can I remove it ?
A line plot would be nice.

Thanks again,

-- 
Oz Nahum
Graduate Student
Zentrum für Angewandte Geologie
Universität Tübingen

---

Imagine there's no countries
it isn't hard to do
Nothing to kill or die for
And no religion too
Imagine all the people
Living life in peace



More information about the SciPy-User mailing list