James Pelletier
MAS.864 Nature of Mathematical Modeling
Random Systems

(6.1a) We begin with the series expansions for moments and cumulants $$\begin{eqnarray} \left< e^{ikx} \right> &=& \sum_{n = 0}^\infty \frac{(ik)^n}{n!} \left< x^n \right> \nonumber \\ \ln \left< e^{ikx} \right> &=& \sum_{n_m = 0}^\infty \frac{(ik)^{n_m}}{n_m!} C_{n_m} \end{eqnarray}$$
To express cumulants in terms of moments we should equate the above expressions, so we first exponentiate the cumulants series expansion $$\exp \left( \sum_{n_m = 1}^\infty \frac{(ik)^{n_m}}{n_m!} C_{n_m} \right) = \sum_{m = 0}^\infty \frac{1}{m!} \left( \sum_{n_m = 1}^\infty \frac{(ik)^{n_m}}{n_m!} C_{n_m} \right)^m.$$
We next distribute the products. Then, rather than write the above expression as a sum over powers $m$, we can change our perspective and write instead as a sum over index sets $\{ n_m \}$ of products over $m$. In other words, this expression describes numerous sums multiplied together - imagine an aisle at a grocery store represents each sum, and we walk from aisle to aisle, adding one item from each aisle to our bag, then at checkout we multiply all the items together (apologies for the stretched metaphor). The entire result equals all possible combinations of items, so long as we take one item from each aisle. $$\sum_{\{ n_m \}} \prod_m \frac{1}{m!} \left( \frac{(ik)^{m n_m}}{(n_m!)^m} C_{n_m}^m \right).$$
To compare powers of $ik$, we must constrain the sum over index sets such that $\sum m n_m = n$, to find the moments in terms of the cumulants $$\left< x^n \right> = \sum_{\{ n_m \}}^{\prime} \prod_m \frac{n!}{m! (n_m!)^m} C_{n_m}^m.$$
We can interpret this equation graphically. Each moment is a sum of products of cumulants, with each term multiplied by the number of groupings of $n$ points into connected clusters of size $m$ with redundancy $n_m$. The factor $n!$ describes the number of orderings of $n$ points, $m!$ accounts for orderings of connected clusters of the same size, and $(n_m!)^m$ accounts for orderings within each connected cluster. For example, we can write the first three moments as $$\begin{eqnarray} \left< x \right> &=& C_1 \nonumber \\ \left< x^2 \right> &=& C_2 + C_1^2 \nonumber \\ \left< x^3 \right> &=& C_3 + 3 C_2 C_1 + C_1^3 \end{eqnarray}$$
We solve this system of equations for the cumulants in terms of the moments to find $$\begin{eqnarray} C_1 &=& \left< x \right> \nonumber \\ C_2 &=& \left< x^2 \right> - \left< x \right>^2 \nonumber \\ C_3 &=& \left< x^3 \right> - 3 \left< x^2 \right> \left< x \right> + 2 \left< x \right>^3 \end{eqnarray}$$
The first cumulant is called the mean, the second is called the variance, and the third is called the skewness.

(6.1b) For a Gaussian distribution, we complete the square to calculate the characteristic function $$\begin{eqnarray} \left< e^{ikx} \right> &=& \frac{1}{\sqrt{2 \pi \sigma^2}} \int_{- \infty}^\infty dx \, \exp \left( - \frac{(x - \bar x)^2}{2 \sigma^2} + ikx \right) \nonumber \\ &=& \frac{1}{\sqrt{2 \pi \sigma^2}} \int_{- \infty}^\infty dx \, \exp \left( - \frac{(x - \bar x - \sigma^2 ik)^2}{2 \sigma^2} + ik \bar x + \frac{(ik)^2}{2} \sigma^2 \right) \nonumber \end{eqnarray}$$
We integrate over the translated Gaussian to yield a factor of $\sqrt{2 \pi \sigma^2}$ which cancels the normalization, then we take the natural logarithm to find $$\ln \left< e^{ikx} \right> = ik \bar x + \frac{(ik)^2}{2} \sigma^2.$$
Looking at the coefficients in the series expansion with respect to $ik$, we find the mean of the Gaussian distribution is $\bar x$, the variance is $\sigma^2$, and all higher cumulants are zero.

(6.2a) For small displacements, we use derivatives as linear approximations. For example, we find the change in $y_1$ as $$dy_1 = \frac{\partial y_1}{\partial x_1} dx_1 + \frac{\partial y_1}{\partial x_2} dx_2.$$
We can express this linear transformation in terms of the Jacobian matrix $J$, in which the $i^{\rm th}$ row and $j^{\rm th}$ column stores the derivative of $y_i$ with respect to $x_j$. The determinant of the Jacobian matrix gives the area of the mapped differential element.

(6.2b) To find the transformed probability distribution, we first calculate the matrix of partial derivatives $$\begin{eqnarray} \frac{\partial y_1}{\partial x_1} &=& - \frac{\sin x_2}{x_1 \sqrt{-2 \ln x_1}} \nonumber \\ \frac{\partial y_2}{\partial x_1} &=& - \frac{\cos x_2}{x_1 \sqrt{-2 \ln x_1}} \nonumber \\ \frac{\partial y_1}{\partial x_2} &=& \sqrt{-2 \ln x_1} \cos x_2 \nonumber \\ \frac{\partial y_1}{\partial x_1} &=& - \sqrt{-2 \ln x_1} \sin x_2 \nonumber \end{eqnarray}$$
The determinant of the Jacobian is $1 / x_1$, and $p(y_1, y_2) = p(x_1, x_2) / \det(J)$. If $p(x_1, x_2)$ is uniform, then $p(y_1, y_2) \propto x_1$. Solving for $x_1$ in terms of $y_1$ and $y_2$, with normalization we find $$p(y_1, y_2) = \frac{1}{2 \pi} \exp \left( - \frac{y_1^2 + y_2^2}{2} \right).$$
The transformation maps a uniform probability distribution to a two-dimensional Gaussian distribution with mean $(0, 0)$ and variance $(1, 1)$.

(6.2c) The uniform random number generator written in Python (txt, py) uses a linear feedback shift register (LFSR) to generate random bits. The code reads this table of lag values. For a given order $M$, the program loads ideal taps from a table, to make a "maximal" LFSR that exhibits a flat power spectrum, output uniformity, and other good properties. The program computes each new bit based on the current set of bits, and then the new bit shifts all other bits one step to the right and knocks off the last bit. From one step to the next, all but one bits remain and lead to correlated successive LFSR outputs (see figure, left subplot). For example, if the value of the LFSR at time $t$ is $v$, then the value at time $t+1$ will be either $v/2$ (if new bit 0) or $2^{M-1} + v/2$ (if new bit 1). Thus, to decrease correlations, per random number we iterate $M$ times to populate the register with all new bits (see figure, right subplot).
We then map these uniform random numbers to Gaussian random numbers, as prescribed in part (b). $x_1$ corresponds to some radius; we generate $x_1$ between 0 and 1 to keep the argument of the square root positive. $x_2$ corresponds to some angle; we generate $x_2$ between 0 and $2 \pi$. We calculate the first, second, and third cumulants of the distribution. As expected the computed mean is $(0.05, -0.02) \approx (0, 0)$, the computed variance is $(1.02, 1.01) \approx (1, 1)$, and the computed skewness is $(0.05, -0.09) \approx (0, 0)$.

(6.3a) For an order 4 maximal LFSR, the bit sequence is

(6.3b) If an LFSR has a clock rate of 1 GHz, then the register generates 1 number per nanosecond. There are approximately $3.2 \times 10^7$ seconds in a year, so the LFSR should generate about $3.2 \times 10^{26}$ numbers per universe. Since the repeat time of an LFSR is $2^M - 1$, the register must have $M = \log_2 \left( 3.2 \times 10^{26} \right) \approx 88$ bits.

(6.4a) Consider a probability distribution $p(x, t)$ governed by the differential equation $$\frac{\partial p}{\partial t} = D \frac{\partial^2 p}{\partial x^2}.$$
Due to the translational symmetry, we can use a Fourier transform then integrate by parts twice to solve this differential equation. When we integrate by parts, the boundary term equals zero due to the boundedness of the probability distribution. $$\begin{eqnarray} \int dx \, \frac{\partial p}{\partial t} e^{- i k x} &=& \int dx \, D \frac{\partial^2 p}{\partial x^2} e^{- i k x} \nonumber \\ \frac{\partial}{\partial t} \tilde p(k, t) &=& - D k^2 \tilde p(k, t) \end{eqnarray}$$
Therefore, $\tilde p(k, t) \propto \tilde p (k, 0) e^{- D k^2 t}$. We decompose the initial distribution into Fourier modes, and then each mode exponentially decays with a time scale proportional to the square of the wavelength (longer wavelength modes decay slower than shorter wavelength modes) and inversely proportional to the diffusion coefficient $D$. We note our solution looks similar to the solution of the Schrödinger equation via separation of variables, but exponential decay rather than a complex phase factor describes the time evolution. We next Fourier transform to find the initial probability distribution in frequency space in terms of the initial distribution in position space. We assume the initial distribution is a normalized delta function at the origin. The Fourier transform of a Dirac delta function is a uniform function with magnitude proportional to the system size. In order to work with normalizable functions, we model the delta function initial distribution as a Gaussian with narrow width. As in (6.1b), we complete the square to find $$\begin{eqnarray} \tilde p(k, 0) &\propto& \frac{1}{\sqrt{2 \pi}} \cdot \frac{1}{\sqrt{2 \pi \sigma^2}} \int_{- \infty}^\infty dx \, \exp \left( - \frac{x^2}{2 \sigma^2} - ikx \right) \nonumber \\ &=& \sqrt{ \frac{\sigma^2}{2 \pi}} \exp \left( - \frac{k^2 \sigma^2}{2} \right) \end{eqnarray}$$
A narrow distribution in position space (small $\sigma^2$) corresponds to a wide distribution in frequency space (large $\sigma^{-2}$).

(6.4b) We combine this result with the previous to find the full probability distribution in position space. $$p(x, t) \propto \frac{1}{\sqrt{2 \pi}} \int_{- \infty}^\infty dk \, \tilde p(k, t) e^{ikx}$$
We complete the square to obtain the normalized Gaussian distribution. $$p(x, t) = \frac{1}{\sqrt{2 \pi \left( \sigma^2 + 2 D t \right)}} \exp \left[ - \frac{x^2}{2 \left( \sigma^2 + 2 D t \right)} \right]$$
In the limit the initial distribution approaches a Dirac delta and the initial variance $\sigma^2 \rightarrow 0$, the time evolved variance equals $2 D t$, the same as an individual random walker in one dimension.

(6.4c) We use the Einstein relation to calculate the diffusion coefficient of a sphere of radius $a$ is inversely proportional to the fluid viscosity $\mu$. $$D = \frac{kT}{F / v} = \frac{kT}{6 \pi \mu a}$$

(6.4d) We simulated 10 random walkers over 1000 steps. All the walkers start at 0, and at each time step each walker has equal probability of making a step of $\pm 1$. The simulation written in Python (txt, py) uses the linear feedback shift register (LFSR) from (6.2c) to generate random steps. We plot the trajectories of the walkers and superimpose $3 \sigma(t)$ error curves ($1.5 \, \sigma$ above, $1.5 \, \sigma$ below), based on the analytical expression for the variance $\sigma(t) = \sqrt{2 D t} = \sqrt{t}$, where $t$ represents the time step, since $D = \left< \delta^2 \right> / 2 \tau$ and $\delta = \tau = 1$.

(6.4e) We integrate the Gaussian probability over $\pm 1.5 \, \sigma$ and find the error bars should enclose about 87% of the trajectories. Seems consistent with simulation!