Diffusion -- Solutions to problems

from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
from collections import defaultdict
plt.ion()


d = 3
L = 100
n = 1000

#
def make_trajectory(self_avoiding=True, steps = 500):
    box = np.zeros([L]*d, dtype=bool)
    x0 = np.array([L//2]*d, dtype=int)
    x = np.array(list(x0))
    box[tuple(x)]=True
    trials = np.zeros((n,d), dtype=int)
    trials[np.arange(n), np.random.randint(d,size=n)] = 1-2*np.random.randint(2,size=n)
    success=0
    xtraj = [x0]
    tp = [0]
    for step in trials:
        xtest = x+step
        if np.any(xtest<0) or np.any(xtest>=L):
            break

        if not box[tuple(xtest)] or self_avoiding==False:
            box[tuple(xtest)]=True
            x = xtest
            success+=1
            xtraj.append(x)
            tp.append(success)
            if success>=steps:
                break

    xtraj = np.array(xtraj)
    return np.array(tp), xtraj - xtraj[0]

# generate an ensemble of trajectories
end2end=[]
bins=30
counts = np.zeros(bins)
distance = np.zeros(bins)
final_distances = []
dt = n/bins
T = np.arange(bins)*dt
for ii in range(5000):
    tp, xtraj = make_trajectory(steps=n/2)
    if tp[-1]>=n/2:
        final_distances.append(xtraj[-1])
    for bi in range(bins):
        ind = tp//dt==bi
        counts[bi]+=np.sum(ind)
        distance[bi]+=np.sum(np.sqrt(np.sum(xtraj[ind]**2, axis=1)))

distance/=counts
plt.figure()
plt.plot(T,distance**2, label="self avoiding")


counts_nsw = np.zeros(bins)
distance_nsw = np.zeros(bins)
final_distances_nsw = []
for ii in range(5000):
    tp, xtraj = make_trajectory(self_avoiding=False)
    if tp[-1]>=n/2:
        final_distances_nsw.append(xtraj[-1])
    for bi in range(bins):
        ind = tp//dt==bi
        counts_nsw[bi]+=np.sum(ind)
        distance_nsw[bi]+=np.sum(np.sqrt(np.sum(xtraj[ind]**2, axis=1)))

distance_nsw/=counts_nsw
plt.plot(T,distance_nsw**2, label="random")
plt.legend()
plt.ylabel('distance squared')
plt.xlabel('time')

### plot end-to-end distribution

plt.figure()
y, x = np.histogram(np.sqrt(np.sum(np.array(final_distances)**2, axis=1)), bins=20)
bc = 0.5*(x[1:] + x[:-1])
plt.plot(bc, y/y.max(),label='self avoiding')

y, x = np.histogram(np.sqrt(np.sum(np.array(final_distances_nsw)**2, axis=1)), bins=20)
bc = 0.5*(x[1:] + x[:-1])
plt.plot(bc, y/y.max(), label='random/no self avoidance')

t=500
y= bc**2*np.exp(-3*bc**2/2/t)
plt.plot(bc,y/y.max(), label='radial Gaussian')

plt.ylabel('distribution [scaled]')

plt.xlabel('end to end distance')