In [1]:
# Libraries are in parent directory
import sys
sys.path.append('../')
from dtrw import *
import numpy as np
import matplotlib.pyplot as plt
In [43]:
# First we have a series of estimates use in proving that the DTRW scheme is stable in the sense of von Neumann
alpha = 0.75
N = 500
d = DTRW_subdiffusive([0.,0.], N, alpha)
V = 0.5
eps = np.zeros(N)
eps_bnd = np.zeros(N)
hopeful_K_telescope_bound = np.zeros(N)
real_eps_bnd = np.array([ alpha / pow(float(i), 1.-alpha) for i in range(1,N+1)])
better_real_eps_bnd = np.array([ alpha / pow(float(i), alpha) for i in range(1,N+1)])
K_eps_bnd = np.array([ d.K[:i].sum() for i in range(2,N+1)])
#eps_inv_bnd = np.array([ 1. / (1. - alpha / pow(float(i), 1.-alpha)) for i in range(1,N+1)])
eps_inv_bnd = 1. / (1. - V * K_eps_bnd)
eps_inv_bnd = 1. / (1. - V * better_real_eps_bnd)
div_bnd = np.array([ pow(float(i+1)/float(i), 1.-alpha) for i in range(1,N+1)])
eps[0] = alpha
#plt.plot(eps_bnd)
#plt.plot(div_bnd)
#plt.plot(div_bnd * div_bnd)
#plt.show()
# For n = 6 (i.e we are calculating eps[6] = \eps_{7})
#
# i eps K eps_prod eps_inv_bnd bnd_prod K_eps_bnd K_inv_prod
# [ 0 1 0 1 nil 1 nil (=1.0)
# 1 2 1 nil 2 6 1+2 5 (=1+2+3+4+5+6)
# 2 3 2 6 3 6*5 1+2+3 5*4 (=1+2+3+4+5+6 * 1+2+3+4+5)
# 3 4 3 6*5 4 6*5*4 1+2+3+4 5*4*3
# 4 5 4 6*5*4 5 6*5*4*3 1+2+3+4+5
# 5 6 5 6*5*4*3 6 6*5*4*3*2 1+2+3+4+5+6
# 6 - 6 6*5*4*3*2 7 6*5*4*3*2*1 1+2+3+4+5+6+7
# 7 - 7 6*5*4*3*2*1 8 -
#print 'eps_inv_bnd', eps_inv_bnd
for n in range(1,N):
eps_inv = 1. / (1. - V * eps)
# n means we're calculating eps_{n+1}, that is eps[n] = eps_{n+1}
# also eps_prod[n] = 1 / (1-eps_{1})...(1-eps_{n})
eps_inv_prod = np.array([eps_inv[n-i:n].prod() for i in range(n+1)])
eps[n] = (eps_inv_prod * d.K[1:n+2]).sum()
eps_inv_bnd_prod = np.array([eps_inv_bnd[n-i:n].prod() for i in range(n+1)])
eps_bnd[n] = (eps_inv_bnd_prod * d.K[1:n+2]).sum()
# TODO: CHECK THIS ONE OUT!
hopeful_K_telescope_bound[n] = np.array([d.K[i+1] * pow(float(n)/float(n-i), 1.-alpha) for i in range(n)]).sum()
#ratio_bnd = np.array([pow(float(n) / float(n-i), 1. - alpha) for i in range(n+1)])
#print 'eps_inv_bnd_prod', eps_inv_bnd_prod
#print 'K[1:n+2]', d.K[1:n+2]
#print 'eps_bnd[n]', eps_bnd[n]
#print 'eps', eps[n], '\t bnd', eps_bnd[n], "\t hope", hopeful_K_telescope_bound[n]
#print eps_bnd[:n-1], bnd_prod, d.K[1:n+1]
#print (bnd_prod * d.K[1:n+1]).sum(), '\t', (div * d.K[1:n+1]).sum()
print 'eps', eps[:10]
print 'eps_inv', eps_inv[:10]
print 'eps_inv_bnd', eps_inv_bnd[:10]
print 'K_eps_bnd', K_eps_bnd[:10]
print 'K_eps_bnd_inv', 1./(1.-K_eps_bnd[:10])
print 'eps_bnd_K_sum', eps_bnd[:10]
plt.plot(eps)
#plt.plot(hopeful_K_telescope_bound,'^')
#plt.plot(eps_bnd, '^')
#plt.plot(real_eps_bnd,'x')
#plt.plot(better_real_eps_bnd,'o')
#plt.plot(K_eps_bnd, 'v')
plt.show()
In [7]:
# Now we look at stability tests: is it really r < 1/2 that makes everything tick?
alpha = 0.5
N = 100
L = 40
X_init = np.zeros(L)
X_init[L/2] = float(L)
omega = 0.5
r = 1.
dif = DTRW_diffusive(X_init, N, omega, r=r)
sub = DTRW_subdiffusive(X_init, N, alpha, r=r)
dif.solve_all_steps()
sub.solve_all_steps()
# Lets look at the first few solutions to see if there's the characteristic oscillations...
print sub.Xs[0][:,:,0]
print sub.Xs[0][:,:,1]
print sub.Xs[0][:,:,2]
print sub.Xs[0][:,:,3]
print sub.Xs[0][:,:,4]
In [8]:
alpha = 0.75
N = 50
d = DTRW_subdiffusive([0.,0.], N, alpha)
b = [(1.-alpha) / pow(i, alpha) for i in range(1,N+1)]
b2 = [alpha / pow(i, alpha) for i in range(1,N+1)]
b = [(1.-alpha) / pow(i, 1.-alpha) for i in range(1,N+1)]
b2 = [alpha / pow(i, 1.-alpha) for i in range(1,N+1)]
plt.plot(-d.K)
plt.plot(b)
plt.plot(b2)
plt.show()
In [47]:
# Double checking the convergence of that ratio thing.
N=100
xi = np.zeros(N)
V = 4.0
xi[0] = 1. - V * alpha
for n in range(1,N):
xi_tel_prod = np.array([ 1. / xi[n-i:n].prod() for i in range(n+1) ])
xi[n] = 1. - V * (xi_tel_prod * d.K[1:n+2]).sum()
plt.plot(xi)
plt.show()
# This calculation is based on the Z-transform - should get the following
print 'According to z-transform, we should get V->', -1. / (pow(1.-1./xi[-1], -alpha) - 1. ), 'actual V =', V
Out[47]:
In [ ]:
# Plotting the relationship between fixed xi and V...
alpha = 0.7
xi = np.arange(-2., 2., 0.005)
# From the paper we get the following relationship:
V = 1. / (1.-pow(1.-1./xi, -alpha))
plt.plot(xi, V)
plt.show()
print alpha
print 1. / (1 - pow(2., -alpha))
In [12]:
alpha = 0.75
N = 500
L = 100
X_init = np.zeros(L)
X_init[L/2] = float(L)
xs = np.arange(0., 1., 1./float(L))
dtrw = DTRW_subdiffusive_with_death(X_init, N, alpha, 0.9, r = 0.1)
dtrw.solve_all_steps()
In [13]:
print xs.shape
print dtrw.Xs[0].shape
plt.plot(xs, dtrw.Xs[0][0,:,100])
plt.show()
In [16]:
n=4
print dtrw.Xs[0][:,:,n].sum()
print np.exp(-0.9*n)*100.
In [33]:
t = X_init > 1.
t.sum() == True
Out[33]:
In [4]:
alpha = 0.75
N = 500
d = DTRW_subdiffusive([0.,0.], N, alpha)
V = 0.5
K = np.zeros(N)
K[0] = 0
K[1] = alpha
for n in range(2,N):
K[n] = alpha * (alpha-1) / (n * (n-1))
for k in range(n-2):
K[n] *= (1.0 + alpha / float(k+1))
print d.K[:10]
print K[:10]
In [ ]: