提交: 0f4d28f35042a997ebf4bc6c416233b3ed133098 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修订版: 0f4d28f35042a997ebf4bc6c416233b3ed133098 (tree)
• 时间: 2007-10-17 07:03:09
• 作者: iselllo
• Commiter: iselllo

Log Message

I added a linear interpolation on the boundary to estimate where the
particle gets deposited and also wrote a routine which saves only 2
positions (present and immediately precedent one) of a Brownian particle
to save memory.

更改概述

• modified: Python-codes/sde.py

差异

diff -r f55b71d65be7 -r 0f4d28f35042 Python-codes/sde.py
--- a/Python-codes/sde.py	Mon Oct 15 09:55:12 2007 +0000
+++ b/Python-codes/sde.py	Tue Oct 16 22:03:09 2007 +0000
@@ -172,6 +172,257 @@
 		return final_pos
 
 
+
+
+
+
+
+class many_particles_2D_box_interpolation:
+	#I initialize the position of the particle
+	def __init__(self, x0,y0,l_boun,r_boun,low_boun,up_boun,N_part):
+		self.x0 = x0
+		self.y0 = y0
+		self.l_boun=l_boun
+		self.r_boun=r_boun
+		self.low_boun=low_boun
+		self.up_boun=up_boun
+		self.N_part=N_part
+		#self.y0 = y0
+	#then I calculate its area
+	def position(self, dt):
+		posx_temp=self.x0
+		posy_temp=self.y0
+		exit_pos=s.zeros(2)
+		#dt=t/N_step
+		#dWx=dWy=s.zeros(1)
+		#dWx=self.x0
+		#dWy=self.y0
+		Wx=s.zeros(1)+self.x0
+		Wy=s.zeros(1)+self.y0
+		randx=n.random.normal(0.,1.,1)
+		randy=n.random.normal(0.,1.,1)
+		while (posx_temp<self.r_boun and posx_temp>self.l_boun and \
+			  posy_temp<self.up_boun and posy_temp>self.low_boun):
+			randx=n.random.normal(0.,1.,1)
+			randy=n.random.normal(0.,1.,1)
+			dWx=s.sqrt(dt)*randx
+			dWy=s.sqrt(dt)*randy
+			Wx=s.append(Wx,(Wx[-1]+dWx))
+			Wy=s.append(Wy,(Wy[-1]+dWy))
+			
+			posx_temp=Wx[-1]
+			posy_temp=Wy[-1]
+		#exit_pos[0]=Wx[-1]
+		#exit_pos[1]=Wy[-1]
+		#Now I think about obtaining the exit position by interpolation
+		a=(Wy[-1]-Wy[-2])/(Wx[-1]-Wx[-2]) #linear coefficient of the straight line connecting
+		# the last value inside the boundary with the first value outside the boundary
+		if (posy_temp>self.up_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.up_boun
+		elif (posy_temp<self.low_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.low_boun
+		elif(posx_temp>self.r_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.r_boun+b
+			exit_pos[0]=self.r_boun
+			exit_pos[1]=y_inter
+		elif(posx_temp<self.l_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.l_boun+b
+			exit_pos[0]=self.l_boun
+			exit_pos[1]=y_inter
+		#results=s.column_stack((Wx,Wy))
+		return exit_pos
+	
+	def iteration_on_particles(self,dt):
+		final_pos=s.zeros((self.N_part,2))
+		for i in xrange(0,self.N_part):
+			final_pos[i]=self.position(dt)
+		return final_pos
+
+
+
+
+
+
+class many_particles_2D_box_interpolation_time:
+	#I initialize the position of the particle
+	def __init__(self, x0,y0,l_boun,r_boun,low_boun,up_boun,N_part):
+		self.x0 = x0
+		self.y0 = y0
+		self.l_boun=l_boun
+		self.r_boun=r_boun
+		self.low_boun=low_boun
+		self.up_boun=up_boun
+		self.N_part=N_part
+		#self.y0 = y0
+	#then I calculate its area
+	def position(self, dt):
+		posx_temp=self.x0
+		posy_temp=self.y0
+		exit_pos=s.zeros(3)
+		#dt=t/N_step
+		#dWx=dWy=s.zeros(1)
+		#dWx=self.x0
+		#dWy=self.y0
+		Wx=s.zeros(1)+self.x0
+		Wy=s.zeros(1)+self.y0
+		randx=n.random.normal(0.,1.,1)
+		randy=n.random.normal(0.,1.,1)
+		while (posx_temp<self.r_boun and posx_temp>self.l_boun and \
+			  posy_temp<self.up_boun and posy_temp>self.low_boun):
+			randx=n.random.normal(0.,1.,1)
+			randy=n.random.normal(0.,1.,1)
+			dWx=s.sqrt(dt)*randx
+			dWy=s.sqrt(dt)*randy
+			Wx=s.append(Wx,(Wx[-1]+dWx))
+			Wy=s.append(Wy,(Wy[-1]+dWy))
+			
+			posx_temp=Wx[-1]
+			posy_temp=Wy[-1]
+		#exit_pos[0]=Wx[-1]
+		#exit_pos[1]=Wy[-1]
+		#Now I think about obtaining the exit position by interpolation
+		a=(Wy[-1]-Wy[-2])/(Wx[-1]-Wx[-2]) #linear coefficient of the straight line connecting
+		# the last value inside the boundary with the first value outside the boundary
+		if (posy_temp>self.up_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.up_boun
+			exit_pos[2]=len(Wy)*dt
+		elif (posy_temp<self.low_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.low_boun
+			exit_pos[2]=len(Wy)*dt
+		elif(posx_temp>self.r_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.r_boun+b
+			exit_pos[0]=self.r_boun
+			exit_pos[1]=y_inter
+			exit_pos[2]=len(Wy)*dt
+		elif(posx_temp<self.l_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.l_boun+b
+			exit_pos[0]=self.l_boun
+			exit_pos[1]=y_inter
+			exit_pos[2]=len(Wy)*dt
+		#results=s.column_stack((Wx,Wy))
+		#print 'exit pos is,' , exit_pos
+		return exit_pos
+	
+	def iteration_on_particles(self,dt):
+		final_pos=s.zeros((self.N_part,3))
+		for i in xrange(0,self.N_part):
+			final_pos[i]=self.position(dt)
+		return final_pos
+
+
+
+
+class many_particles_2D_box_interpolation_time_short_array:
+	#I initialize the position of the particle
+	def __init__(self, x0,y0,l_boun,r_boun,low_boun,up_boun,N_part):
+		self.x0 = x0
+		self.y0 = y0
+		self.l_boun=l_boun
+		self.r_boun=r_boun
+		self.low_boun=low_boun
+		self.up_boun=up_boun
+		self.N_part=N_part
+		#self.y0 = y0
+	#then I calculate its area
+	def position(self, dt):
+		posx_temp=self.x0
+		posy_temp=self.y0
+		exit_pos=s.zeros(3)
+		#dt=t/N_step
+		#dWx=dWy=s.zeros(1)
+		#dWx=self.x0
+		#dWy=self.y0
+		Wx=s.zeros(2)  +self.x0
+		Wy=s.zeros(2)   +self.y0
+		#Wx=s.zeros[0]  +self.x0
+		#Wy=s.zeros[0]   +self.y0
+		randx=n.random.normal(0.,1.,1)
+		randy=n.random.normal(0.,1.,1)
+		time=0.0
+		while (posx_temp<self.r_boun and posx_temp>self.l_boun and \
+			  posy_temp<self.up_boun and posy_temp>self.low_boun):
+			randx=n.random.normal(0.,1.,1)
+			randy=n.random.normal(0.,1.,1)
+			dWx=s.sqrt(dt)*randx
+			dWy=s.sqrt(dt)*randy
+			#Wx=s.append(Wx,(Wx[-1]+dWx))
+			#Wy=s.append(Wy,(Wy[-1]+dWy))
+			Wx[0]=Wx[1]    #now Wx (and Wy) are arrays with only 2 elements and they contain the 
+			                #present position of the particle and the previous one (nothing else
+					# is needed for the evolution, since the particle is Brownian and has no 
+					# memory) 
+			Wx[1]=Wx[1]+dWx
+			Wy[0]=Wy[1]
+			Wy[1]=Wy[1]+dWy
+			
+			
+			posx_temp=Wx[-1]
+			posy_temp=Wy[-1]
+			time=time+dt # I update the time; now I have to calculate it this way since I no longer have
+			# the length of the array which I can use
+			
+			
+			
+		#exit_pos[0]=Wx[-1]
+		#exit_pos[1]=Wy[-1]
+		#Now I think about obtaining the exit position by interpolation
+		a=(Wy[-1]-Wy[-2])/(Wx[-1]-Wx[-2]) #linear coefficient of the straight line connecting
+		# the last value inside the boundary with the first value outside the boundary
+		if (posy_temp>self.up_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.up_boun
+			exit_pos[2]=time+dt
+		elif (posy_temp<self.low_boun):
+			b=Wy[-1]-a*Wx[-1]
+			x_inter=(Wy[-1]-b)/a
+			exit_pos[0]=x_inter
+			exit_pos[1]=self.low_boun
+			exit_pos[2]=time+dt
+		elif(posx_temp>self.r_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.r_boun+b
+			exit_pos[0]=self.r_boun
+			exit_pos[1]=y_inter
+			exit_pos[2]=time+dt
+		elif(posx_temp<self.l_boun):
+			b=Wy[-1]-a*Wx[-1]
+			y_inter=a*self.l_boun+b
+			exit_pos[0]=self.l_boun
+			exit_pos[1]=y_inter
+			exit_pos[2]=time+dt
+		#results=s.column_stack((Wx,Wy))
+		#print 'exit pos is,' , exit_pos
+		return exit_pos
+	
+	def iteration_on_particles(self,dt):
+		final_pos=s.zeros((self.N_part,3))
+		for i in xrange(0,self.N_part):
+			final_pos[i]=self.position(dt)
+		return final_pos
+
+
+
+
+
+
 #my_particle=particle_1D_new(2.) #initialize particle position
 #now I can send a message to the object
 
@@ -193,7 +444,7 @@
 
 
 twoD=1
-many=1
+many=4
 
 if (twoD ==1):
 	if (many==0):
@@ -222,7 +473,7 @@
 		p.savefig('Brownian_xy.pdf')
 	elif(many==1):
 		print "here I am!"
-		my_particles=many_particles_2D_box(0.,0.,-1.,1.,-1.,1.,1000) # inject 10 Brownian particles
+		my_particles=many_particles_2D_box(0.,0.,-1.,1.,-1.,1.,10000) # inject 10 Brownian particles
 		exit_positions=my_particles.iteration_on_particles(dt)
 		print "the exit positions are, ", exit_positions
 		p.save('exit_positions',exit_positions)
@@ -273,6 +524,215 @@
 		
 		print'my_count is, ', my_count
 		print 'the total number of escaped particles is, ', sum(my_count)
+	elif(many==2):
+		print "here I am!"
+		my_particles=many_particles_2D_box_interpolation(0.,0.,-1.,1.,-1.,1.,100) # inject 10 Brownian particles
+		exit_positions=my_particles.iteration_on_particles(dt)
+		print "the exit positions are, ", exit_positions
+		p.save('exit_positions',exit_positions)
+		#Now I record the exit positions along the 4 sides of the box
+		
+		
+		#the particle escapes from the top (y particle >1)
+		up=s.where(exit_positions[:,1]>=1.)
+		up_exit=exit_positions[up,0]
+		
+		#the particle escapes from the bottom (y particle < -1) 
+		low=s.where(exit_positions[:,1]<=-1.)
+		low_exit=exit_positions[low,0]
+		
+		
+		#the particle escapes from the left (x particle < -1) 
+		left=s.where(exit_positions[:,0]<=-1.)
+		left_exit=exit_positions[left,1]
+		
+		# the particle escapes from the right (x particle > 1)
+		right=s.where(exit_positions[:,0]>=1.)
+		right_exit=exit_positions[right,1]
+		
+		print "right_exit is ", right_exit
+		
+		p.hist(right_exit)
+		p.savefig('right_exit.pdf')
+		p.hold(False)
+		
+		p.hist(left_exit)
+		p.savefig('left_exit.pdf')
+		p.hold(False)
+		
+		p.hist(up_exit)
+		p.savefig('up_exit.pdf')
+		p.hold(False)
+		
+		p.hist(low_exit)
+		p.savefig('low_exit.pdf')
+		p.hold(False)
+		
+		#now I keep track of the # of particles leaving from each side
+		my_count=s.zeros(4)
+		my_count[0]=s.shape(right_exit)[1]
+		my_count[1]=s.shape(left_exit)[1]
+		my_count[2]=s.shape(up_exit)[1]
+		my_count[3]=s.shape(low_exit)[1]
+		
+		print'my_count is, ', my_count
+		print 'the total number of escaped particles is, ', sum(my_count)
+	
+	elif(many==3):
+		print "here I am!"
+		my_particles=many_particles_2D_box_interpolation_time(0.,0.,-1.,1.,-1.,1.,100) # inject 10 Brownian particles
+		exit_positions=my_particles.iteration_on_particles(dt)
+		#print "the exit positions are, ", exit_positions
+		p.save('exit_positions',exit_positions)
+		#Now I record the exit positions along the 4 sides of the box
+		
+		
+		#the particle escapes from the top (y particle >1)
+		up=s.where(exit_positions[:,1]>=1.)
+		up_exit=exit_positions[up,0]
+		time_up=exit_positions[up,2]
+		
+		#the particle escapes from the bottom (y particle < -1) 
+		low=s.where(exit_positions[:,1]<=-1.)
+		low_exit=exit_positions[low,0]
+		time_low=exit_positions[low,2]
+		
+		#the particle escapes from the left (x particle < -1) 
+		left=s.where(exit_positions[:,0]<=-1.)
+		left_exit=exit_positions[left,1]
+		time_left=exit_positions[left,2]
+		
+		# the particle escapes from the right (x particle > 1)
+		right=s.where(exit_positions[:,0]>=1.)
+		right_exit=exit_positions[right,1]
+		time_right=exit_positions[right,2]
+		
+		
+		
+		print "right_exit is ", right_exit
+		
+		p.hist(right_exit)
+		p.savefig('right_exit.pdf')
+		p.hold(False)
+		
+		p.hist(left_exit)
+		p.savefig('left_exit.pdf')
+		p.hold(False)
+		
+		p.hist(up_exit)
+		p.savefig('up_exit.pdf')
+		p.hold(False)
+		
+		p.hist(low_exit)
+		p.savefig('low_exit.pdf')
+		p.hold(False)
+		
+		#Now the plots of the exit times
+		
+		p.hist(time_right)
+		p.savefig('right_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_left)
+		p.savefig('left_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_up)
+		p.savefig('up_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_low)
+		p.savefig('low_exit_time.pdf')
+		p.hold(False)
+		
+		#now I keep track of the # of particles leaving from each side
+		my_count=s.zeros(4)
+		my_count[0]=s.shape(right_exit)[1]
+		my_count[1]=s.shape(left_exit)[1]
+		my_count[2]=s.shape(up_exit)[1]
+		my_count[3]=s.shape(low_exit)[1]
+		
+		print'my_count is, ', my_count
+		print 'the total number of escaped particles is, ', sum(my_count)
+		
+	elif(many==4):
+		print "here I am!"
+		my_particles=many_particles_2D_box_interpolation_time_short_array(0.,0.,-1.,1.,-1.,1.,100) # inject 10
+		                                                                              #Brownian particles
+		exit_positions=my_particles.iteration_on_particles(dt)
+		#print "the exit positions are, ", exit_positions
+		p.save('exit_positions',exit_positions)
+		#Now I record the exit positions along the 4 sides of the box
+		
+		
+		#the particle escapes from the top (y particle >1)
+		up=s.where(exit_positions[:,1]>=1.)
+		up_exit=exit_positions[up,0]
+		time_up=exit_positions[up,2]
+		
+		#the particle escapes from the bottom (y particle < -1) 
+		low=s.where(exit_positions[:,1]<=-1.)
+		low_exit=exit_positions[low,0]
+		time_low=exit_positions[low,2]
+		
+		#the particle escapes from the left (x particle < -1) 
+		left=s.where(exit_positions[:,0]<=-1.)
+		left_exit=exit_positions[left,1]
+		time_left=exit_positions[left,2]
+		
+		# the particle escapes from the right (x particle > 1)
+		right=s.where(exit_positions[:,0]>=1.)
+		right_exit=exit_positions[right,1]
+		time_right=exit_positions[right,2]
+		
+		
+		
+		print "right_exit is ", right_exit
+		
+		p.hist(right_exit)
+		p.savefig('right_exit.pdf')
+		p.hold(False)
+		
+		p.hist(left_exit)
+		p.savefig('left_exit.pdf')
+		p.hold(False)
+		
+		p.hist(up_exit)
+		p.savefig('up_exit.pdf')
+		p.hold(False)
+		
+		p.hist(low_exit)
+		p.savefig('low_exit.pdf')
+		p.hold(False)
+		
+		#Now the plots of the exit times
+		
+		p.hist(time_right)
+		p.savefig('right_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_left)
+		p.savefig('left_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_up)
+		p.savefig('up_exit_time.pdf')
+		p.hold(False)
+		
+		p.hist(time_low)
+		p.savefig('low_exit_time.pdf')
+		p.hold(False)
+		
+		#now I keep track of the # of particles leaving from each side
+		my_count=s.zeros(4)
+		my_count[0]=s.shape(right_exit)[1]
+		my_count[1]=s.shape(left_exit)[1]
+		my_count[2]=s.shape(up_exit)[1]
+		my_count[3]=s.shape(low_exit)[1]
+		
+		print'my_count is, ', my_count
+		print 'the total number of escaped particles is, ', sum(my_count)
+		
 
 elif (twoD==0):
 	if (box_1D==1):