import numpy as np
import scipy
import ipywidgets as widgets
import matplotlib.pyplot as plt
%matplotlib inline
d=3
t= np.linspace(0, 1, 10000)
for D in [0.01, 0.1, 1, 10]:
plt.plot(t, 2*d *D * t**0.5)
#plt.loglog(t, 2*d *D * t**0.5)
plt.ylabel('$MSD(t)$')
plt.xlabel('$t$')
Text(0.5, 0, '$t$')
:::{figure-md} markdown-fig
Animation of brownian motion :::
def brown(T, N, dt=1, D=1, d=3):
"""
Creates 3D brownian path given:
time T
N=1 trajecotires
dt=1 timestep
D=1 diffusion coeff
returns np.array with shape (N, T, 3)
"""
n = int(T/dt) # how many points to sample
dR = np.sqrt(2*D*dt) * np.random.randn(N, n, d) # 3D position of brownian particle
R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
return R
R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
def brownian_plot(t=10):
fig, ax = plt.subplots(ncols=2)
ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
ax[1].hist(R[:, 10, 0], density=True, color='red')
ax[1].hist(R[:, t, 0], density=True)
ax[1].set_ylim([0,0.1])
ax[0].set_ylim([-200, 200])
ax[0].set_xlim([-200, 200])
fig.tight_layout()
import holoviews as hv
hv.extension('plotly')
plots = []
for i in range(10):
plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
plots.append(plot)
hv.Overlay(plots)
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]
xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])
hv.Overlay(rw_curve) << ydist << xdist
Diffusion equation is a 2nd order PDE! Unlike Netwon's or Schrodinger equations diffusion equation shows irreversible behavior.The diffusion coefficient $D$ has of units $[L^2]/[T]$.
For the case of free diffusion the equation is exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically. Can solve PDE numerically by writing derivatives as finite differences!
An important special case solution is free diffusion for 1D case. Notice that it is gaussian with MSD that we derived in CLT and LLN sections $\sigma(t)=\sqrt{2{D}t}$. You may heck that this is indeed a solution by plugging into the diffusion equation!
def sigma(t, D = 1):
return np.sqrt(2*D*t)
def gaussian(x, t):
return 1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
def diffusion(t=1):
R = brown(T=101, N=1000)
x = np.linspace(-20, 20, 100)
plt.plot(x, gaussian(x, t), label=f'Diff eq solution at t={t}', color='black')
plt.hist(R[:,t,0], density=True, label='brownian diffusion sim', color='green', alpha=0.6)
plt.legend()
plt.ylabel('$p(x)$')
plt.xlabel('$x$')
plt.xlim([-25, 25])
plt.ylim([0, 0.5])
widgets.interact(diffusion, t=(1, 100))
interactive(children=(IntSlider(value=1, description='t', min=1), Output()), _dom_classes=('widget-interact',)…
<function __main__.diffusion(t=1)>