FRAP with a box shaped pulse.
The assignment was
Work out the solution of the diffusion equation for initial condition $$ P(x,0) = \cases{ c & x < -a \cr 0 & -a < x < a \cr c & x > a } $$
In the lecture, we have worked out the solution to one-sided case where the initial condition is \(c\) for \(x<0\) and zero otherwise. Since the diffusion equation is linear, solutions add up and we can use the previous solution to construct one that satisfies these boundary conditions
$$ P(x,t) = \frac{c}{\sqrt{4\pi Dt}} \left[\int_{-\infty}^{-a} dx_0 e^{-\frac{(x_0-x)^2}{4Dt}} + \int_{a}^{\infty} dx_0 e^{-\frac{(x_0-x)^2}{4Dt}}\right] = \frac{c}{\sqrt{4\pi Dt}} \left[\int_{-\infty}^{-a-x} dx_0 e^{-\frac{x_0^2}{4Dt}} + \int_{a-x}^{\infty} dx_0 e^{-\frac{x_0^2}{4Dt}}\right] $$ This can be expressed as error functions as before.
Diffusion equation in a closed box
Find a general solution to the diffusion equation in a box with no-flux boundaries at \(x=\pm L\). (Hint: Separate variables, trigonometric functions fulfill the boundary conditions.)
This problem is most readily solved by separation of variables. With this method, one is looking for a solution of the form \(P(x,t) = g(x)f(t)\). Plugging this into the diffusion equation, we find
$$g(x)\frac{d f(t)}{dt} = f(t) \left[D\frac{d^2g(x)}{dx^2} - v\frac{dg(x)}{dx}\right] $$
We now divide by \(g(x)f(t)\) to find
$$\frac{1}{f(t)}\frac{d f(t)}{dt} = -\lambda = \frac{1}{g(x)}D\frac{d^2g(x)}{dx^2}$$
Note that the left side depends only on \(t\) while the right side depends only on \(x\). The only way to make this happen is for both sides to be constant. In this case, the two equations can be solved and we find
$$f(t) = e^{-\lambda t} \quad g(x) = a \cos(\sqrt{\lambda/D}x) + b \sin(\sqrt{\lambda/D}x) $$
The boundary conditions now determine what kind of \(\lambda\) values are admissible. For no-flux boundaries at \(x=-L\) and \(x=L\), we require either \(b=0\) and
$$ \sqrt{\lambda/D}L = n \pi \Rightarrow \lambda_n = \frac{n^2\pi^2 D}{L^2} $$
or \(a=0\) and (lets call the eigenvalue \(\kappa\) instead of \(\lambda\))
$$ \sqrt{\kappa/D}L = (2n+1)\pi/2 \Rightarrow \kappa_n = \frac{(n+1/2)^2\pi^2 D}{L^2} $$
We can now write the general solution as
$$P(x,t) = \sum_{n=0}^\infty a_n e^{-\lambda_n t} \cos(\sqrt{\lambda_n/D}x) + \sum_{n=0}^\infty b_n e^{-\kappa_n t} \sin(\sqrt{\kappa_n/D}x) $$
We can now express any initial condition as a superposition of these functions with different \(a_n\) and \(b_n\).
$$P(x,t) = \sum_{n=0}^\infty a_n e^{-\lambda_n t} \cos(\frac{n\pi x}{L}) + \sum_{n=0}^\infty b_n e^{-\kappa_n t} \sin(\frac{(2n+1)\pi x}{2L}) $$
The trigonometric functions with these eigen values from a complete orthogonal basis such that any arbitrary initial condition can be represented by a suitable choice of \(a_n\) and \(b_n\). The different \(n\) are often called modes that decay at different rates. The slowest mode, \(\lambda_0\) doesn't decay at all. It corresponds to a flat stationary solution to which the system will decay. All other modes fade away over time, and the rate at which they fade is inversely proportional to the wave length of their spatial oscillation.
Return to the origin of a random walker.
from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
from collections import defaultdict
plt.ion()
def return_intervals(pos):
'''
determine the intervals between successive returns to the origin
'''
intervals = {}
intervals[1] = np.concatenate([np.diff(np.where(pos[:,i]==0)[0]) for i in range(pos.shape[1])])
intervals[2] = np.concatenate([np.diff(np.where((pos[:,i]==0)&(pos[:,j]==0))[0]) for i,j in [(0,1), (0,2), (1,2)]])
intervals[3] = np.concatenate(np.diff(np.where(np.all(pos==0, axis=1))))
return intervals
def return_times(pos):
'''
determine the times at which the walker returns to the origin
'''
intervals = {}
intervals[1] = np.concatenate([np.where(pos[:,i]==0)[0] for i in range(pos.shape[1])])
intervals[2] = np.concatenate([np.where((pos[:,i]==0)&(pos[:,j]==0))[0] for i,j in [(0,1), (0,2), (1,2)]])
intervals[3] = np.concatenate(np.where(np.all(pos==0, axis=1)))
return intervals
# generate multiple random walks and analyze them for their return properties.
t_max = 10000
d = 3
all_intervals = defaultdict(list)
all_returns = defaultdict(list)
for rep in range(100):
position = np.cumsum(1-2*np.random.randint(2, size=(t_max, d)), axis=0)
for start in range(0,t_max//10, 30):
tmp = return_intervals(position[start:] - position[start])
for i in tmp:
all_intervals[i].extend(tmp[i])
tmp = return_times(position[start:] - position[start])
for i in tmp:
all_returns[i].extend(tmp[i])
plt.figure()
x=np.concatenate([np.arange(2,10,2), np.arange(10,100,8), np.logspace(2,4,20)])
dx = np.diff(x)
for i in all_returns:
y, _x = np.histogram(all_returns[i], bins=x)
n = len(all_returns[i])
plt.plot(0.5*(x[1:]+x[:-1]), y/dx/n, label = "dim %d, n=%d"%(i, n))
plt.plot([1,1000], [1, 0.001**0.5])
plt.plot([1,1000], [1, 0.001])
plt.plot([1,1000], [1, 0.001**1.5])
plt.ylim([0.00001,2])
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.ylabel('probability of return')
plt.ylabel('time')
plt.figure()
x=np.concatenate([np.arange(2,10,2), np.arange(10,100,8), np.logspace(2,4,20)])
dx = np.diff(x)
for i in all_intervals:
plt.plot(sorted(all_intervals[i]), np.linspace(1,0,len(all_intervals[i])), label = "dim %d, n=%d"%(i, len(all_intervals[i])))
plt.plot([1,1000], [1, 0.001**0.5])
plt.plot([1,1000], [1, 0.001])
plt.ylim([0.001,2])
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.ylabel('fraction shorter')
plt.ylabel('interval')