# [SciPy-user] Integrate.odeint

Lorenzo Isella lorenzo.isella@gmail....
Wed Apr 11 11:48:27 CDT 2007

Dear All,
First of all I want to make clear that I am not looking for someone
taking care of my homework.
I am struggling to solve a system of nonlinear ODE's (mentioned in a
previous mail of mine).
I am resorting to integrate.odeint which worked wonderfully for me in the past.
What I am finding puzzling is that I can understand getting the wrong
result, but not somehow ending up with a solution which does not
respect the initial condition.
I cut and paste the code.
Does anyone have a clue at what is going on?
Kind Regards

Lorenzo

#! /usr/bin/env python
from scipy import *
import pylab  # used to read the .csv file

x=linspace(1.,300.,20) #set of initial particle diameters
lensum=len(x) # number of particle bins
myvar1=1.5 #standard deviation of the log-normal distribution
mu1=70. # mean of the distribution
A1=1000. #amplitude of the distribution
print 'x is', x

def distr(A1,mu1,myvar1,x): # function representing the initial
log-normal distribution
z=log(10.)*A1/sqrt(2.*pi)/log(myvar1)*exp(-((log(x/mu1))**2.) \
/2./log(myvar1)/log(myvar1))
return z

vec_distr=vectorize(distr) # I vectorized the previous functio

y0=vec_distr(A1,mu1,myvar1,x) # initial condition of log-normally
distributed particles

#print 'n_ini is', n_ini

#I plot the initial state condition

pylab.plot(x,y0)
pylab.xlabel('D_p')
pylab.ylabel('n_ini(0)')
#pylab.legend(('prey population','predator population'))
pylab.title('initial distribution')
pylab.grid(True)
pylab.savefig('N_initial')

pylab.hold(False)

# this is the system of 1st order ODE's I want to solve:
# \frac{dy[k]}{dt}=0.5*sum_{i+j=k}kernel[i,j]*y[i]*y[j] (creator) +
#     -y[k]sum_{i=1}^infty kernel[i,k]*y[i] (destructor)

# NB: careful since in the formula both i and j start from 1!

# In the following, I will be using a trivial constant kernel to test the code

kern=1e-6 # kernel

def coupling(y,t,kernel,lensum):
out=zeros(lensum) # array which I'll use to write down the
differential equations
creation=zeros(lensum)
destruction=zeros(lensum)

for k in range(0,lensum):
for i in range(0,lensum):
destruction[k]=destruction[k]-kern*y[i]
for j in range(0,lensum):
if (i+1+j+1==k+1): #I add 1 to correct the array indexing
creation[k]=creation[k]+kern*y[i]*y[j]

destruction[k]=y[k]*destruction[k]
creation[k]=0.5*creation[k]
if (creation[k]+destruction[k]>=0.):
out[k]=creation[k]+destruction[k] # output array

return out

t=arange(0.,100.,1.)
print 't is', t

#t_fix=0.
#test=coupling(y0,t_fix,kern,lensum)

#print 'test is', test

ysol = integrate.odeint(coupling, y0, t,args=(kern,lensum),printmessg=1)

print 'the shape of y is', shape(ysol)

pylab.plot(x,ysol[0,:]/max(ysol[0,:]),x,ysol[90,:]/max(ysol[90,:]))
pylab.xlabel('D_p')
pylab.ylabel('population')
pylab.title('evolution')
pylab.grid(True)
pylab.savefig('N_evolved')

print 'the solu at t=0 is', ysol[0,:]

print 'and the initial condition is', y0

print 'So far so good'