[SciPy-User] weired results with ode solver

[Warren Weckesser]

Oz,

In rate_func(), try changing this:

	co=c1[0]#; %concentration of Oxygen
	cs=c1[1]#; %concentration of Substrate
	X=c1[2]#; %concentration of bacteria


to this

	co = c[0] #; %concentration of Oxygen
	cs = c[1] #; %concentration of Substrate
	X  = c[2] #; %concentration of bacteria


Also, you have all your parameters hard-coded in rate_func().   That's 
fine if you want to work that way, but you could also add additional 
arguments to rate_func(), and use the 'args' keyword argument of odeint 
to tell the solver what the additional arguments are.  There is an 
example in the SciPy cookbook:

    http://www.scipy.org/Cookbook/CoupledSpringMassSystem


Regards,

Warren


Oz Nahum wrote:
> 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
>