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

6.1 - Cumulants

(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

6.4

(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$

(c) - How is the diffusion coefficient of Brownian motion related to the viscosity of a fluid?

(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}$