[SciPy-User] weired results with ode solver

Oz Nahum nahumoz@gmail....
Sat Feb 6 11:32:16 CST 2010


Hi Everyone,

I'm trying to convert the following matlab code into python, the code
simulates reduction of organic matter by bacteria in a batch system.
It has to files in the matlab version, and 1 in the python version.
I'm not sure I did the conversion correctly, but I'm sure the result
I'm getting from python is wrong - probably because I did not really
understand how to use the solver. Any insights would be appreciated.

Here are the codes:
% This file is a matlab script to simulate a batch problem as described

% in problem 1a.

% This file calls batch1_rate.m


%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=[Coi;Csi;Cbi];


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]


%Solve ODE

timerange=[0 43*86400];

[t,c]=ode15s(@batch1_rate,timerange,c1,[],ko,ks,Ys,Yo,Kdec,umax,kb);


%Plot

hf=figure(1);

clf;

%plot(,c(:,1:3),'linewidth',2);

days=t/86400;

plot(days,c(:,1),'b-',days,c(:,2),'k--','linewidth',2);

hold on

plot(days,c(:,3),'-','Color',[0 0.498 0],'linewidth',1);


xlabel('T [days]','fontsize',14);

ylabel('C [mg/L]','fontsize',14);

ylim([0 10]);

xlim([0 22])


#end of 1st matlab file.
% This file is an aid function to the script batch1.m which describes

% a batch system relations of substrate, oxygen and one microbial specie

% over time.

% Oz Nahum, Anibal Perez, Georgia Kastanioti, Januar 2010.





function f=batch1_rate(t,c1,ko,ks,Ys,Yo,Kdec,umax,kb)

co=c1(1); %concentration of Oxygen

cs=c1(2); %concentration of Substrate

X=c1(3); %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=[foxi;fsub;fbio];

end

#end of second matlab file

#python script:
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=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]
	co=c1[0]#; %concentration of Oxygen
	cs=c1[1]#; %concentration of Substrate
	X=c1[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,60)
parms=(c1,ko,ks,Ys,Yo,Kdec,umax,kb)
#[t,c]=ode15s(@batch1_rate,timerange,c1,[],ko,ks,Ys,Yo,Kdec,umax,kb);
#c=odeint(rate_func,parms,timerange)
c=odeint(rate_func,c1,timerange)
#sprint c

from pylab import *
plot(timerange, c[:,0], 'b-')
plot(timerange, c[:,1], 'g-')
plot(timerange, c[:,2], 'k-')
show()



Many thanks in advance.


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