weired results with ode solver
Oz Nahum <nahumoz <at> gmail.com>
2010-02-06 17:32:16 GMT
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
RSS Feed