--- jupytext: formats: ipynb,md:myst text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.6 kernelspec: display_name: Python 3 (math-583) language: python name: math-583 --- ```{code-cell} :tags: [hide-cell] import mmf_setup; mmf_setup.nbinit() import os from pathlib import Path import logging logging.getLogger("matplotlib").setLevel(logging.CRITICAL) logging.getLogger('numba').setLevel(logging.CRITICAL) logging.getLogger('OMP').setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt ``` (sec:Solution11)= # Solution 11: Brownian Motion (Random Walks) Here we present a partial solution of {ref}`sec:Assignment11`. For details about Brownian motion and random walks, please see the following (part of my [Standard Model course][]): * [Renormalization Group Techniques](https://physics-581-the-standard-model.readthedocs.io/en/latest/ClassNotes/RenormalizationGroup.html) * [Renormalizing Random Walks][] The assignment asks you to consider Brownian motion where the steps $x_n$ are independent and identically and independently distributed ([iid][]) with some probability density function bracket ([PDF][]) $p(x)$. The distribution $p_N(X)$ after $N$ steps is simply the distribution obtained by summing $N$ random variables $x$ with [PDF][] $p(x)$. As discussed in [Renormalizing Random Walks][], the [PDF][] of the sum of random variables is simply the [convolution][] of the individual [PDF][]s. This means that the Fourier transform of $p_N(X)$ is just the Fourier transform of $p(x)$ raised to the $N$th power. :::{margin} This integral can be done by completing the square, changing variables, then using the standard result \begin{gather*} \int_{-\infty}^{\infty} e^{-x^2/2}\d{x} = \sqrt{2\pi}. \end{gather*} ::: For the gaussian distribution this is simple: \begin{gather*} \tilde{p}(k) = \int_{-\infty}^{\infty}\d{x}\;p(x)e^{-\I k x} = \frac{1}{\sqrt{2\pi}\sigma} \int_{-\infty}^{\infty}\d{x}\;e^{-\I k x - (x-\mu)^2/2\sigma^2}\\ = \exp\Bigl(-\I \mu k - \sigma^2 \frac{k^2}{2}\Bigr),\\ \tilde{p}_{N}(k) = \tilde{p}(k)^N = \exp\Bigl( - \I \underbrace{N \mu}_{\mu_N} k - \underbrace{N\sigma^2}_{\sigma_N^2} \frac{k^2}{2}\Bigr). \end{gather*} :::{margin} This property, that $\tilde{p}_{N} = \tilde{p}^N$ explains the central limit theorem. If we rescale our variable $k$ so that either the mean ($k \rightarrow k/N$) or the variance ($k \rightarrow k/\sqrt{N}$) stay fixed, then all higher order terms vanish in the limit of large $N$: most general distributions become gaussian. See [Renormalizing Random Walks][] for details. ::: Thus, the position after $N$ steps $X_{N}$ is normally distributed with mean $N\mu$ and standard deviation $\sigma_N = \sqrt{N}\sigma$: \begin{gather*} X_{N} \sim \mathcal{N}\bigl(\mu_N = N\mu, \sigma_N = \sqrt{N}\sigma\bigr). \end{gather*} For the binary walk, the number of ways of ending up at position $X_N = m-n$ is the number of ways of choosing $m$ steps to be $+1$ and $n = N-m$ steps to be $-1$. This gives the [binomial distribution][] scaled by 2 and shifted by $-N$ \begin{gather*} \binom{N}{m}p^{m}(1-p)^{N-m}\\ p_N(X) = \sum_{m=0}^{N}\binom{N}{m}p^{m}(1-p)^{N-m}\delta(X - \underbrace{(2m - N)}_{m-n}), \end{gather*} whose mean and standard deviation are: :::{margin} The standard [binomial distribution][] has \begin{gather*} \mu_N = Np, \qquad \sigma_N = \sqrt{N}\sqrt{p(1-p)}, \end{gather*} where we take step $+1$ with probability $p$ and no step with probability $1-p$, yielding values from $0$ to $N$. Here we must multiply by $2$ then shift by $-N$ to obtain values from $-N$ to $N$. ::: \begin{gather*} \mu_N = N(2p-1), \qquad \sigma_N = \sqrt{N} 2\sqrt{p(1-p)}. \end{gather*} This gives the answer about how to relate the binary random walk and the gaussian random walk: \begin{gather*} \mu = 2p-1, \qquad \sigma = 2\sqrt{p(1-p)}. \end{gather*} :::::{admonition} Another Approach :class: dropdown We can follow the same approach we used for a gaussian, considering the Fourier transform of the binary distribution: \begin{gather*} \tilde{p} = pe^{-\I k} + (1-p)e^{\I k}. \end{gather*} Expanding for small $k$ we have: \begin{gather*} \tilde{p} \approx 1 - (2p - 1)\I k - \frac{k^2}{2} + O(k^3)\\ \approx e^{-\mu\I k - \sigma^2 k^2/2 + O(k^3)} \approx 1 - \mu\I k - (\mu^2 + \sigma^2) k^2/2 + O(k^3)\\ \mu = 2p-1, \qquad \sigma = \sqrt{1-\mu^2} = 2\sqrt{p(1-p)}, \end{gather*} as stated above. ::::: ```{code-cell} :tags: [hide-input] def get_binary_walk(N, p=0.5, seed=2): """Return a binary random walk of N steps.""" rng = np.random.default_rng(seed=seed+10) # Slight trick: we get a binary distribution from the uniform # distribution on [0, 1]. dx = np.where(rng.uniform(size=N) < p, 1, -1) x = np.cumsum(dx) return x def get_normal_walk(N, mean=0, stdev=1, seed=2): """Return a gaussian random walk of N steps.""" rng = np.random.default_rng(seed=seed+10) dx = rng.normal(loc=mean, scale=stdev, size=N) x = np.cumsum(dx) return x ps = [0.2, 0.5, 0.9] fig, axs = plt.subplots(1, 2, figsize=(4*2, 4)) N = 500 Ns = 1000 step = 1+np.arange(N) for n, p in enumerate(ps): mu, sigma = 2*p-1, 2*np.sqrt(p*(1-p)) x_bs = np.array([get_binary_walk(N, p=p, seed=n) for n in range(Ns)]) x_ns = np.array([get_normal_walk(N, mean=mu, stdev=sigma, seed=n) for n in range(Ns)]) kw = dict(c=f"C{n}", alpha=0.5) label = r"$N \mu $" if n == 0 else None axs[0].plot(step, step*mu, '-', label=label, **kw) axs[0].plot(step, x_bs.mean(axis=0), '--', label=f"Binary ({p=})", **kw) axs[0].plot(step, x_ns.mean(axis=0), ':', **kw) label = r"$\sqrt{N}\sigma$" if n == 0 else None axs[1].plot(step, np.sqrt(step)*sigma, '-', label=label, **kw) axs[1].plot(step, x_bs.std(axis=0), '--', **kw) axs[1].plot(step, x_ns.std(axis=0), ':', label=f"Normal ({mu=:.2f}, {sigma=:.2f})", **kw) axs[0].set(xlabel="step $N$", ylabel=r"$\langle X_N \rangle$"); axs[1].set(xlabel="step $N$", ylabel=r"$\sigma_N$"); axs[0].legend() axs[1].legend() plt.suptitle(f"Averaged over {Ns=} samples"); ``` [Renormalizing Random Walks]: [Standard Model course]: [iid]: [PDF]: [convolution]: [binomial distribution]: