weired results with ode solver

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
_______________________________________________
SciPy-User mailing list
SciPy-User <at> scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user