[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'


More information about the SciPy-user mailing list