Demo entry 5367846

Test

   

Submitted by anonymous on Jun 21, 2016 at 01:11
Language: Python 3. Code size: 2.7 kB.

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages

# Set parameters
T_end = 10.0
steps= np.power(2,11)
dt= T_end/steps
samples = 2

# parameter for lspde:  dX_t = b*X_t dt + sigma*dW_t
sigma = 1
b = 1.5
x_0 = 0.0

alpha = 0.5

seed = np.random.randint(100000)
#or using a fixed seed
#seed = 95097
np.random.seed(seed)
print ('seed of random generator: ',seed)

#
## init pot -------------------------------------------------------------------------------------
#
#init the figure
fig = plt.figure(figsize=[25,20])
ax1 = plt.subplot(411,autoscale_on=False, xlim=(0, T_end), ylim=(-5, 1))
ax2 = plt.subplot(412)
ax3 = plt.subplot(413,autoscale_on=False, xlim=(0, T_end), ylim=(-5, 1))
ax4 = plt.subplot(414)

ax1.set_xlabel('time')
ax1.set_ylabel('X^1(t)')
ax1.set_title('sample of 1.particle drift '+str(b)+' diff '+str(sigma))
ax1.grid(True)
ax2.set_xlabel('time')
ax2.set_title('sample of counting proc. 1.particle')
ax2.grid(True)
ax3.set_xlabel('time')
ax3.set_ylabel('X(t)')
ax3.set_title('empiricale mean of particles')
ax3.grid(True)
ax4.set_xlabel('time')
ax4.set_title('sample of counting proc. des 2.Partikels')
ax4.grid(True)

#
## simulation -----------------------------------------------------------------------------------
#Euler-Maruyama method on linear SDE, i.e.
#  dX_t = b*X_t dt + sigma* dW_t
#   X_0 = x_0
# mit Ruecksprungregel

Phi = np.zeros( (samples,steps) )
Count = np.zeros( (samples,steps) )

Phi[:,0] = x_0*np.ones( samples )
for t in xrange(steps-1):
    Phi[:,t+1] = Phi[:,t] + dt*b+ sigma*np.sqrt(dt)*np.random.randn(samples) + dt*alpha/samples*np.sum(Count[:,t]) - dt*Count[:,t]
    Count[:,t+1] = Count[:,t]
    for j in range(samples):
	if Phi[j,t+1] >= 1:
	    Phi[j,t+1] = 0;
	    Count[j,t+1] = Count[j,t]+1;
	    for k in range(samples):
		if k != j and Phi[k,t+1] < (1-alpha/samples):
		    Phi[k,t+1] = Phi[k,t+1] + alpha/samples
	    

#fuer Konvergenze des empirischen Mittels untersuche Prozess ohne Ruecksprung !!!!


#
## pot ------------------------------------------------------------------------------------------
#
#ax1.plot(np.linspace(0,T_end,T_end),1*np.ones(T_end),'r')
ax1.plot(np.linspace(0,T_end,steps),Phi[0,:],np.linspace(0,T_end,T_end),(1-alpha/samples)*np.ones(T_end),'r')
ax2.plot(np.linspace(0,T_end,steps),Count[0,:])
ax3.plot(np.linspace(0,T_end,steps),Phi[1,:],np.linspace(0,T_end,T_end),(1-alpha/samples)*np.ones(T_end),'r')
ax4.plot(np.linspace(0,T_end,steps),Count[1,:])

##grafic output
name = 'simu_'+str(b)+'_'+str(sigma)+'_'+str(alpha)+'.pdf'
pp = PdfPages(name)
pp.savefig()
pp.close()

plt.show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).