In [564]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from scipy import stats


### (a) - The first three cumulants of a distributuion¶

We can use the equivelence of the power series expansion of the characteristic function with $e^{log\langle e^{ikx} \rangle}$

$exp\left( \sum_{n=1}^\infty \frac{ik^n}{n!}C_n \right) = \sum_{n=1}^\infty \frac{ik^n}{n!} \langle x^n \rangle$ using the series expansion of $e^x$ $e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} ...$ we have

$$1 + x + \sum_{n=1}^\infty \frac{ik^n}{n!}C_n + \frac{( \sum_{n=1}^\infty \frac{ik^n}{n!}C_n)^2}{2!} + \frac{( \sum_{n=1}^\infty \frac{ik^n}{n!}C_n)^3}{3!} ... = \sum_{n=1}^\infty \frac{ik^n}{n!} \langle x^n \rangle$$

seperating the first order terms $\frac{ik}{1}C_1 = \frac{ik}{1} \langle x \rangle \\ C_1 = \langle x \rangle$

seperating the second order terms $-k^2 C_2+\frac{-k^2\:C_1^2}{2!} = \frac{-k^2}{2!} \langle x^2 \rangle \\ \frac{-k^2 C_2}{2!}+\frac{-k^2\:\langle x \rangle^2}{2!} = \frac{-k^2}{2!} \langle x^2 \rangle \\ C_2+\langle x \rangle^2 = \langle x^2 \rangle\\ C_2 = \langle x^2 \rangle - \langle x \rangle^2$

seperating the third order terms $\frac{-ik^3 C_3}{3!}+\frac{-ik^3\:C_1^3}{3!} = \frac{-ik^3}{3!} \langle x^3 \rangle \\ \frac{-ik^3 C_3}{3!}+\frac{-ik^3\:\langle x \rangle^3}{3!} = \frac{-ik^3}{3!} \langle x^3 \rangle \\ C_3+\langle x \rangle^3 = \langle x^3 \rangle\\ C_3 = \langle x^3 \rangle - \langle x \rangle^3$

### (b) Cumulants of a guassian¶

$F(t) = \int_{-\infty}^\infty P(x) \; e^{xt} dt = \langle e^{tx} \rangle \\ F(t) = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi{\sigma^2}}} e^{\frac{-(x^2)}{2{\sigma}^2}} \; e^{xt} dx \\ F(t) = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi{\sigma^2}}} e^{\frac{-(x^2)}{2{\sigma}^2}+xt} dx$

using the guassian integral

$\int_{-\infty}^\infty e^{-\alpha x^2 + \beta x}\;dx = \sqrt{\frac{\pi}{\alpha}}e^{\frac{\beta^2}{4\alpha}}\\ \alpha = \frac{-1}{2{\sigma}^2} \quad \text{and} \quad \beta = t \\ F(t) = e^{\frac{t^2{\sigma}^2}{2}} \\ log(F(t)) = \frac{t^2\sigma^2}{2}$

We see that there are no terms of a higher order than 2

### 6.2 - Transforming a uniform random number generator¶

In [ ]:



### 6.3 - Linear Feedback Shift Register¶

In [586]:
def LFSR(x):
#outbit = (np.sum([x[1],x[4]])%2)
outbit = (x[1]+x[4])%2
y = ([outbit,x[0],x[1],x[2],x[3]])
return(y)

In [587]:
inbits = ([0,1,1,1,1])
bits = ([0,1,1,1,1])
for i in range(0,100):
inbits = LFSR(inbits)
bits = np.vstack([bits,inbits])

fig = plt.imshow(bits.T)
fig.set_cmap('hot')
plt.ylabel('register values')
plt.xlabel('timesteps')
plt.title('bit values for a 4th order LFSR')
plt.show()

In [588]:
def autocorr(x):
result = np.correlate(x, x, mode='full')
return result[result.size/2:]

xx = bits[:,1]
x = autocorr(xx)
plt.plot(x,'k')
plt.show()

/Users/andrewbahle/anaconda/lib/python3.6/site-packages/ipykernel/__main__.py:3: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
app.launch_new_instance()


### (b) - If an LFSR has a clock rat of 1 GHz how many registers must if have for the repeat sequence length to be the age of the universe?¶

The repeat frequency for a LFSR is $\:2^M-1\:$ where M is the length of the register (excluding the register where bits are read out). If and LFSR has a clock rate of 1 GHz how long must the length of the register for the repeat length to be $\: 10^{10} \:$

$\bigg(\frac{10^9\:\text{GHz}}{1\:\text{Hz}}\bigg)\bigg(\frac{60\:\text{Hz}}{1\:\text{min}}\bigg) \bigg(\frac{60\:\text{min}}{1\:\text{hr}}\bigg)\bigg(\frac{24\:\text{hrs}}{1\:\text{day}}\bigg)\bigg(\frac{365\:\text{days}}{1\:\text{year}}\bigg)*10^{10}\:\text{years} = 2^M-1\\ M = log_2{(3.1536*10^{26}+1)} \\ M = 88.0271...\\ M = 89$

In [561]:
x = float(60*60*24*365*10**10*10**9+1)
np.ceil(np.log2(x))

Out[561]:
89.0

### (a) - solve the diffusion equation using a Fourier transform¶

$ft [f(x)] = \int_{-\infty}^\infty f(x) e^{-ikx} dk \qquad ift [f(k)] = \frac{1}{2\pi}\int_{-\infty}^\infty f(x)e^{-ikx} dx \qquad ft[\frac{d^nf(x)}{dx^n}] = (-ik)^nf(k) \qquad ft(\delta) = 1 \qquad ft(1) = \delta$

$\text{We begin with the diffusion equation}$

$\frac{\partial \hat{p}}{\partial{t}} = \frac{\langle \delta^2 \rangle}{2\tau} \frac{\partial^2p}{\partial{x^2}} \qquad \frac{\partial \hat{p}}{\partial{t}} = D \frac{\partial^2p}{\partial{x^2}}$

$\text{Taking the fourier transform of the diffusion equation and using the derivative rule for fourier transforms we have}$

$\frac{\partial \hat{p}}{\partial t} = (-ik)^2\hat{p} = -Dk^2\hat{p} \qquad \text{which is in the form of an ODE} \qquad \frac{\partial \hat{p}}{\partial t}+Dk^2\hat{p} = 0$

$\text{Using the integrating factor} \quad e^{\int{Dk^2dt}} = e^{Dk^2t} \quad \text{we have}$

$e^{Dk^2t}\frac{\partial \hat{p}}{\partial t}+e^{Dk^2t}Dk^2\hat{p} = 0$

$\text{which is derivative of} \quad e^{Dk^2t} \hat{p}$

$\text{integrating both sides} \quad \int {\frac{\partial}{\partial t} e^{Dk^2t} \hat{p} \:dt} = \int{0}\: dt\\ \text{since}\: \hat{p}\: \text{is a function of both k and t, zero becomes a function of k after integrating}\\ e^{Dk^2t} \hat{p} = g(k) \quad \text{or since dividing by}\: e^x \text{is the same as multiplying by}\: e^{-x} \quad \: \hat{p} = g(k)e^{-Dk^2t}$

$\text{our initial conditions stipulate that at} \: t=0 \: \text{we have a delta function} \: \text{i.e.} \: p(x,0) = \delta \: \text{therefore the fourier transform at} \: t = 0 \: \text{is 1}$ $\text{this simplifies our equation to} \; \hat{p}(x,0) = 1*e^{-Dk^2t}$

So our inverse Fourier transform is:

$ift [f(k)] = \frac{1}{2\pi}\int_{-\infty}^\infty f(x)e^{-ikx} dx \\ p(x) = \frac{1}{2\pi}\int_{-\infty}^\infty e^{-Dk^2t} e^{-ikx} dx \\ p(x) = \frac{1}{2\pi}\int_{-\infty}^\infty e^{-Dk^2t-ikx} dx \\$

This is easy to solve using the guassian integral used above

$\int_{-\infty}^\infty e^{-\alpha x^2 + \beta x}\;dx = \sqrt{\frac{\pi}{\alpha}}e^{\frac{\beta^2}{4\alpha}}\\ \text{where for our ift equation the parameters are} \quad \alpha = Dt \quad \text{and} \quad \beta = ix \\ p(x) = \frac{1}{2\pi}\sqrt{\frac{\pi}{Dt}}e^{\frac{ix^2}{4Dt}}\\ p(x) = \frac{1}{2\pi}\frac{1}{\sqrt{Dt}}e^{\frac{-x^2}{4Dt}}\\ p(x) = \frac{1}{\sqrt{4 \pi Dt}}e^{\frac{-x^2}{4Dt}}$

### (b) - What is the variance as a function of time?¶

$\sigma^2_{t} = \frac{\delta^2}{\tau}*t$

$\sigma^2_{t} = {D}*t$

(The diffusion coefficient for Brownian motion is related to the viscosity of a fluid using Stokes-Einstein Equation).

$D = \frac{k_B T}{6 \pi \eta r}$

Where $\;\eta\;$is the viscosity of the medium $\;r\;$ is the radius of the particle $\;T\;$ is the temperature and $\;k_B\;$ is Boltzmann's constant

### (d) - diffusion simulation¶

In [565]:
n = 1000;
nsteps = 1000;
walkers = np.zeros([n,1])
temp = np.zeros([n,1])
yax = np.arange(1.,n+1.)
x = 1 #initialize random num generator

In [566]:
for steps in range(0,nsteps):
for i in range(0, temp.size):
#x = np.random.rand(1,1)
x = (((8121*x) + 28411)%134456)
if x/134456 > 0.5:
# walker moves right
temp[i] += 1
else:
# walker moves left
temp[i] += -1

walkers = np.append(walkers, temp,1)

In [583]:
fig, ax1 = plt.subplots()

ax1.plot(temp,yax,'k.',ms=2)
ax1.set_xlabel('position')
ax1.axis([-100,100,0,n+1])

ax1.set_ylabel('walkers')
ax1.tick_params('y')

ax2 = ax1.twinx()
ax2.set_xlim([-100,100])
ax2.hist(temp,alpha=0.1,color='blue')
ax2.set_ylabel('frequency',color = 'blue')
ax2.tick_params('y')
plt.title('distribution of 1000 walkers after 1000 steps')

ax3 = ax1.twinx()
xt = plt.xticks()[0]
xmin, xmax = min(xt), max(xt)
xx = np.linspace(xmin, xmax, len(temp))

ax3.axis('off')
mu, sigma = stats.norm.fit(temp) # get mean and standard deviation
pdf_g = stats.norm.pdf(xx, mu, sigma) # now get theoretical values in our interval
ax3.plot(xx, pdf_g, label="Norm") # plot it

plt.show()

In [214]:
xx = np.arange(0.,1000.)
error1 = 1*np.sqrt(np.arange(0.,1000.))
error2 = 2*np.sqrt(np.arange(0.,1000.))
error3 = 3*np.sqrt(np.arange(0.,1000.))

plt.fill_between(xx, error1, -error1, color='blue', alpha='0.3')
plt.fill_between(xx, error2, -error2, color='blue', alpha='0.15')
plt.fill_between(xx, error3, -error3, color='grey', alpha='0.1')

for i in range(0,9):
plt.plot(walkers[i,:],'k',alpha = 0.7)

plt.ylabel('position')
plt.xlabel('timesteps')
plt.title('walker trajectories')
plt.show()


### (e) - What fraction of the trajactories should be contained in the error bars¶

$\text{The standard deviation grows as the square root of time and ~99.7% of the trajectories should lie within the first 3 standard deviations}$

$\sigma_{t} = \sqrt{\frac{\delta^2}{\tau}*t}\\ \sigma_{n} = \sqrt{\frac{1^2}{1}*n}\\ \sigma_{n} = \sqrt{n}$