(6.1)

(a) Work out the first three cumulants $C_1$, $C_2$, and $C_3$.

Using Equations 6.24-6.26 from the textbook, we can expand and solve both sides, paying attention only to terms with $k$ terms of the first, second, and third orders for the first three cumulants.

$exp(\sum_{n=1}^{\infty} \frac{(ik)^n}{n!} C_n) = \sum_{n=0}^{\infty} \langle x^n \rangle$

$\implies e^{(ikC_1 + \frac{(ik)^2}{2}C_2 + \frac{(ik)^3}{6}C_3 + \ldots)} = \langle x_0 \rangle + ik \langle x^1 \rangle + \frac{(ik)^2}{2} \langle x^2 \rangle + \frac{(ik)^3}{6} \langle x^3 \rangle$

$e^x = 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \ldots$

$x = ikC_1 - \frac{1}{2}k^2C_2 - \frac{1}{6}ik^3C_3 + \ldots$

$\implies 1 + (ikC_1 - \frac{1}{2}k^2C_2 - \frac{1}{6}ik^3C_3) + \frac{1}{2}(ikC_1 - \frac{1}{2}k^2C_2 - \frac{1}{6}ik^3C_3)^2 + \frac{1}{6}(ikC_1 - \frac{1}{2}k^2C_2 - \frac{1}{6}ik^3C_3)^3 = \langle 1 \rangle + ik \langle x \rangle - \frac{k^2}{2} \langle x^2 \rangle - \frac{ik^3}{6} \langle x^3 \rangle$

Since we only care about $k^1$, $k^2$, and $k^3$, we will only extract those terms.

$\implies ikC_1 - \frac{1}{2}k^2C_2 - \frac{1}{6}ik^3C_3 - \frac{1}{2}k^2C_1^2 - \frac{1}{2}ik^3C_1C_2 - \frac{1}{6}ik^3C_1^3$

$ikC_1 = ik \langle x \rangle \implies C_1 = \langle x \rangle$

$-\frac{1}{2}k^2C_2 - \frac{1}{2}k^2C_1^2 = -\frac{1}{2}k^2 \langle x^2 \rangle \implies C_2 = \langle x^2 \rangle - C_1^2 \implies C_2 = \langle x^2 \rangle - \langle x \rangle^2$

$-\frac{1}{6}ik^3C_3 - \frac{1}{2}ik^3C_1C_2 - \frac{1}{6}ik^3C_1^3 = -\frac{1}{6}ik^3 \langle x^3 \rangle \implies C_3 = \langle x^3 \rangle - 3C_1C_2 - C_1^3 \implies C_3 = \langle x^3 \rangle - 3 \langle x \rangle \langle x^2 \rangle + 2 \langle x \rangle^3$

Note that the first cumulant is the mean and the second cumulant is the variance.

(b) Evaluate the first three cumulants for a Gaussian distribution

$p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(x - \bar{x})^2/2 \sigma^2}$

$\langle x \rangle = \int_{-\infty}^{\infty} x \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(x - \bar{x})^2/2 \sigma^2} dx = \bar{x}$

$\langle x^2 \rangle = \int_{-\infty}^{\infty} x^2 \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(x - \bar{x})^2/2 \sigma^2} dx = \bar{x}^2 + \sigma^2$

$\langle x^3 \rangle = \int_{-\infty}^{\infty} x^3 \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(x - \bar{x})^2/2 \sigma^2} dx = \bar{x}^3 + 3 \bar{x} \sigma^2$

$C_1 = \langle x \rangle = \bar{x}$

$C_2 = \langle x^2 \rangle - \langle x \rangle^2 = \bar{x}^2 + \sigma^2 - \bar{x}^2 = \sigma^2$

$C_3 = \langle x^3 \rangle - 3 \langle x \rangle \langle x^2 \rangle + 2 \langle x \rangle^3 = \bar{x}^3 + 3 \bar{x} \sigma^2 - 3 \bar{x} (\bar{x}^2 + \sigma^2) + 2 \bar{x}^3 = 0$

(6.2)

(a) If $\overrightarrow{y}(\overrightarrow{x}) = (y_1(x_1,x_2),y_2(x_1,x_2))$ is a coordinate transformation, what is the area of a differential element $dx_1 dx_2$ after it is mapped into the $\overrightarrow{y}$ plane? Recall that the area of a parallelogram is equal to the length of its base times its height.

With the notion of the vectors creating a parallelogram, we have to determine the base and height lengths so we can multiply them. We can easily represent these using the differential elements of the vectors:

$d\overrightarrow{y}(1) = (\frac{\partial y_1}{\partial x_1}dx_1, \frac{\partial y_2}{\partial x_1}dx_1)$

$d\overrightarrow{y}(2) = (\frac{\partial y_1}{\partial x_2}dx_2, \frac{\partial y_2}{\partial x_2}dx_2)$

Therefore, the base is $|d\overrightarrow{y}(1)|$, and the height is $|d\overrightarrow{y}(2)|sin\theta$ ($\theta$ is the angle between the vectors). Multiplying these two is equivalent to finding their cross product.

$|d\overrightarrow{y}(1)||d\overrightarrow{y}(2)|sin\theta = |d\overrightarrow{y}(1) \times d\overrightarrow{y}(2)|$

$ = |\frac{\partial y_1}{\partial x_1} \frac{\partial y_2}{\partial x_2} - \frac{\partial y_1}{\partial x_2} \frac{\partial y_2}{\partial x_1}| dx_1 dx_2$

$ = |\frac{\partial \overrightarrow{y}}{\partial \overrightarrow{x}}| dx_1 dx_2$

Which is the Jacobian of the vectors.

(b) Let

$y_1 = \sqrt{-2 ln x_1} sin(x_2) ,\; y_2 = \sqrt{-2 ln x_1} cos(x_2) \;$.

If $p(x_1,x_2)$ is uniform, what is $p(y_1,y_2)$?

First, we will find the partial derivatives of the two equations:

$\frac{\partial y_1}{\partial x_1} = \frac{-1}{x_1 \sqrt{-2 ln x_1}} sin(x_2)$

$\frac{\partial y_1}{\partial x_2} = \sqrt{-2 ln x_1} cos(x_2)$

$\frac{\partial y_2}{\partial x_1} = \frac{-1}{x_1 \sqrt{-2 ln x_1}} cos(x_2)$

$\frac{\partial y_2}{\partial x_2} = -\sqrt{-2 ln x_1} sin(x_2)$

We then find the magnitude of the determinant of the matrix assembled by these partial derivatives.

$|\frac{\partial y_1}{\partial x_1} \frac{\partial y_2}{\partial x_2} - \frac{\partial y_1}{\partial x_2} \frac{\partial y_2}{\partial x_1}| = |\frac{1}{x_1}sin^2(x_2) + \frac{1}{x_1}cos^2(x_2)|$

$ = \frac{1}{|x_1|} = \frac{1}{e^{-(y_1^2 + y_2^2)/2}}$

This gives us the new distribution

$p(y_1,y_2) = p(x_1,x_2)|\frac{\partial \overrightarrow{y}}{\partial \overrightarrow{x}}|^{-1}$

$ = p(x_1,x_2)e^{-(y_1^2 + y_2^2)/2}$

Therefore, $y_1$ and $y_2$ will have Gaussian distributions if $p(x_1,x_2)$ is constant.

(c) Write a uniform random number generator, and transform it by equation (6.78). Numerically evaluate the first three cumulants of its output.

In [3]:
# Adapted closely from Neil's gauss.py code

import numpy as np

n = 1
N = 100000

m1 = 0
m2 = 0
m3 = 0

a = 8121
b = 28411
c = 134456

for i in range(N):
    n = (a*n + b) % c
    x1 = (1 + n) / c
    n = (a*n + b) % c
    x2 = (1 + n) / c
    
    y1 = np.sqrt(-2*np.log(x1))*np.sin(2*np.pi*x2)   # From Eq. (6.78)
    y2 = np.sqrt(-2*np.log(x1))*np.cos(2*np.pi*x2)   # From Eq. (6.78)
    
    m1 = m1 + y1 + y2
    m2 = m2 + y1**2 + y2**2
    m3 = m3 + y1**3 + y2**3
    
m1 = m1/(2*N)
m2 = m2/(2*N)
m3 = m3/(2*N)

C1 = m1
C2 = m2 - m1**2
C3 = m3 - 3*m1*m2 + 2*m1**3

print("Cumulant 1 =", C1)
print("Cumulant 2 =", C2)
print("Cumulant 3 =", C3)
Cumulant 1 = 0.001229637077087833
Cumulant 2 = 0.99886425322971
Cumulant 3 = 0.005550206343330839

(6.3)

(a) For an order 4 maximal LFSR write down the bit sequence.

We can think of the LFSR as utilizing the following equation:

$x_n = \sum_{i=1}^M a_i x_{n-i} (mod 2)$

Because it's an order 4 maximal LFSR, we know that $M = 4 \implies i = 1,4$

So, essentially, we have a logical XOR statement of $x_n = x_{n-1} \oplus x_{n-4}$.

I used a seed value of "1111" which gives us the following:

step $x_n$ $x_{n-1}$ $x_{n-2}$ $x_{n-3}$ $x_{n-4}$
0 0 1 1 1 1
1 1 0 1 1 1
2 0 1 0 1 1
3 1 0 1 0 1
4 1 1 0 1 0
5 0 1 1 0 1
6 0 0 1 1 0
7 1 0 0 1 1
8 0 1 0 0 1
9 0 0 1 0 0
10 0 0 0 1 0
11 1 0 0 0 1
12 1 1 0 0 0
13 1 1 1 0 0
14 1 1 1 1 0
15 0 1 1 1 1

(b) If an LFSR has a clock rate of 1 GHz, how long must the register be for the time between repeats to be the age of the universe (∼ $10^{10}$ years)?

The number of distinct patterns is $2^N - 1$ where $N$ is the number of bits in the register. We'll use $f$ to signify the clock rate frequency, and $t$ to be time. We get the following equation:

$\frac{(2^N - 1)}{f} = t \implies N = log_2(ft + 1)$

$10^{10} \ \text{years}$ is equivalent to $10^{10} \times 365.2425 \times 24 \times 60 \times 60 = 315569520000000000 \ \text{seconds}$.

$N = log_2(31556952 \times 10^{10} \times 10^9) \approx 88 \ \text{bits}$

(6.4)

(a) Use a Fourier transform to solve the diffusion equation (6.57) (assume that the initial condition is a normalized delta function at the origin).

Equation (6.57):

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

The spatial Fourier transform pair is:

$F(k,t) = \int_{-\infty}^{\infty} e^{ikx} f(x,t) dx$

$f(x,t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} e^{-ikx} F(k,t) dk$

We will substitute this into the diffusion equation, which gives us

$\frac{1}{2 \pi} \int_{-\infty}^{\infty} e^{-ikx} [\frac{\partial F}{\partial t} + k^2 DF] dk = 0$

Which results in

$\frac{dF}{dt} + k^2 DF = 0$

We solve this to find $F = Ae^{-k^2 Dt}$. The inverse transform of $F$ is

$f = A'e^{-x^2/2\sigma^2} = A'e^{-x^2/4Dt}$

After comparing to a normalized Gaussian equation, we see our solution is

$f = \frac{1}{\sqrt{4 \pi Dt}} e^{-x^2/4Dt}$

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

We can find this from part (a): $\sigma^2 = 2Dt$

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

We can compare the variance to the Langevin equation to get

$2Dt = \frac{kTt}{3 \pi \mu d} \implies D = \frac{kT}{6 \pi \mu d}$

(d) Write a program (including the random number generator) to plot the position as a function of time of a random walker in 1D that at each time step has an equal probability of making a step of $\pm 1$. Plot an ensemble of 10 trajectories, each 1000 points long, and overlay error bars of width $3\sigma(t)$ on the plot.

First off, we need to find the value of $\sigma$ so we can properly determine the error width.

$\sigma^2 = 2Dt$

$D = \frac{\langle \delta^2 \rangle}{2 \tau}$

$\sigma^2 = 2 \frac{\langle \delta^2 \rangle}{2 \tau} t$

$\langle \delta^2 \rangle, \; \tau = 1$

$\implies = \sigma^2 = t$

Therefore, $3\sigma$ goes from $\pm 1.5 \sqrt{t}$

Note that for this problem, I used an LCG Random Number Generator with suggested values from Wikipedia:

$a = 1103515245, \; c = 12345, \; m = 2^{32}$

In [25]:
import matplotlib.pyplot as plt
import numpy as np
import time

(a, c, m) = (1103515245, 12345, 2**32)

def seed (x):
    global previous
    previous = (a*int(x) + c) % m;

def rand ():
    global previous;
    x = (a*previous + c) % m;
    previous = x;
    return x/m;

seed(time.time())

def stepper (n):
    x = [0]
    for i in range(n - 1):
        if rand() > 0.5:
            step = 1
        else:
            step = -1
        xn = x[-1] + step
        x.append(xn)
    return x

n = 1000;
for i in range(10):
    x = stepper(n)
    plt.plot(x)
    
plt.ylabel('position')
plt.xlabel('n (step number)')

t = np.arange(0, n, n/10)
y = np.zeros(t.size)
error = 1.5*np.sqrt(t)
plt.errorbar(t, y, yerr=error, capsize=5, fmt="none", color="black")

plt.show()

(e) What fraction of the trajectories should be contained in the error bars?

$\int_{-1.5\sigma}^{1.5\sigma} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-x^2/2\sigma^2} dx = \text{erf}(\frac{3}{2\sqrt{2}}) = 0.866$