\documentclass[12pt]{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{equations}
\usepackage{own}
%\setlength{\textwidth}{6.5in}
%\renewcommand{\baselinestretch}{1.1}
%\renewcommand{\theequation}{\thesection.\arabic{equation}}
% see page 19 of Advanced LaTeX course by Michel Goossens for nifty picture
%\voffset -1in
%\hoffset -1in
\oddsidemargin 0.55625in % Left margin on odd-numbered pages.
\evensidemargin 0in % Left margin on even-numbered pages.
\reversemarginpar % put notes in opposite margin than defaults
\textheight 8.70in
\addtolength{\textheight}{-\baselineskip}
\textwidth 5.98in
\topmargin -.28in % distance from top of page to header - 1in
\headheight 0.5in % height of box containing the header
\headsep \baselineskip % distance between head and body
%\footheight 0.5in % height of foot
\footskip 0.5in % bottom of last line to bottom of foot
\topsep=0in % get rid of extra vertical space in list modes
\partopsep=0in %
\parsep=0.25in %
\parskip 0.1 in
\title{Stochastic Modeling in Physics and Microbiology\\
MATH 6791--1 -- Spring 2013 \\
Homework 4}
\date{Due Wednesday, May 8 at 5:00 PM}
\newcommand{\oh}{\frac{1}{2}}
\newcommand{\grad}{\boldsymbol{\nabla}}
\newcommand{\dive}{\grad \cdot}
\newcommand{\divev}{\grad_{\bv} \cdot}
\newcommand{\divex}{\grad_{\bx} \cdot}
\renewcommand{\theenumi}{\alph{enumi}}
\newcommand{\mathi}{\mathrm{i}}
\newcommand{\expe}{\mathrm{e}}
\newcommand{\difd}{\mathrm{d}}
\newcommand{\Prob}{\,\mathrm{Prob}\,}
\newcommand{\bbE}{\mathbb{E}}
\newcommand{\Ito}{It\^{o}}
\newcommand{\pT}{p_T}
\newcommand{\pTeps}{p_T^{(\epsilon)}}
\newcommand{\xp}{x^{\prime}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bXavg}{\mathbf{X}^{\sharp}}
\newcommand{\bY}{\mathbf{Y}}
\newcommand{\bW}{\mathbf{W}}
\newcommand{\bV}{\mathbf{V}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bv}{\mathbf{v}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\br}{\mathbf{r}}
\newcommand{\bxp}{\mathbf{x}^{\prime}}
\newcommand{\bxprobe}{\mathbf{x}^{(p)}}
\newcommand{\bk}{\mathbf{k}}
\newcommand{\hbuk}{\hat{\mathbf{u}}_{\bk}}
\newcommand{\hbumk}{\hat{\mathbf{u}}_{\bk}}
\newcommand{\hpk}{\hat{p}_{\bk}}
\newcommand{\hqk}{\hat{q}_{\bk}}
\newcommand{\hgk}{\hat{g}_{\bk}}
\newcommand{\bFtherm}{\mathbf{F}_T}
\newcommand{\bfthermk}{\hat{\mathbf{f}}_{T,\bk}}
\newcommand{\bxinit}{\mathbf{x}_{\mathrm{in}}}
\newcommand{\bftherm}{\mathbf{f}_T}
\newcommand{\tp}{t^{\prime}}
\newcommand{\Id}{\mathit{\mathsf{I}}}
\newcommand{\bumeas}{\mathbf{U}_{m}}
\newcommand{\lprobe}{\ell_p}
\newcommand{\lpart}{a}
\newcommand{\bVpart}{\mathbf{V}}
\newcommand{\XXtens}{\mathit{\mathsf{C}}^{(\bX \bX)}}
\newcommand{\XVtens}{\mathit{\mathsf{C}}^{(\bX \bV)}}
\newcommand{\lF}{\ell_F}
\newcommand{\Fmag}{\bar{F}}
\newcommand{\bF}{\mathbf{F}}
\newcommand{\pxv}{p_{\bX,\bV}}
\newcommand{\px}{p_{\bX}}
\newcommand{\kB}{k_B}
\newcommand{\lkk}{\ell_{KK}}
\newcommand{\tkk}{\tau_{KK}}
\newcommand{\lobs}{\ell_{o}}
\newcommand{\tobs}{\tau_{o}}
\newcommand{\bxua}{\bx^{(1)}}
\newcommand{\bxub}{\bx^{(2)}}
\newcommand{\bXua}{\bX^{(1)}}
\newcommand{\bXub}{\bX^{(2)}}
\newcommand{\bVua}{\bV^{(1)}}
\newcommand{\bVub}{\bV^{(2)}}
\newcommand{\ma}{m^{(1)}}
\newcommand{\mb}{m^{(2)}}
\newcommand{\Gmat}{\mathit{\mathsf{\Gamma}}}
\newcommand{\Gaa}{\Gmat^{(11)}}
\newcommand{\Gab}{\Gmat^{(12)}}
\newcommand{\Gba}{\Gmat^{(21)}}
\newcommand{\Gbb}{\Gmat^{(22)}}
\newcommand{\Gij}{\Gmat^{(ij)}}
\newcommand{\sigmat}{\mathit{\mathsf{S}}}
\newcommand{\sigaa}{\sigmat^{(11)}}
\newcommand{\sigab}{\sigmat^{(12)}}
\newcommand{\sigba}{\sigmat^{(21)}}
\newcommand{\sigbb}{\sigmat^{(22)}}
\newcommand{\sigij}{\sigmat^{(ij)}}
\newcommand{\bWua}{\bW^{(1)}}
\newcommand{\bWub}{\bW^{(2)}}
\newcommand{\bWx}{\bW_x}
\newcommand{\bWy}{\bW_y}
\newcommand{\rada}{a_1}
\newcommand{\radb}{a_2}
\newcommand{\bphi}{\bar{\phi}}
\newcommand{\tphi}{\tilde{\phi}}
\newcommand{\tphip}{\tilde{\phi}^{\prime}}
\newcommand{\tphipp}{\tilde{\phi}^{\prime\prime}}
\newcommand{\aavg}{a^{\sharp}}
\newcommand{\pycsx}{p^{(s)}_{\bY|\bX}}
\newcommand{\bYx}{\bY_x}
\newcommand{\Imean}{\bar{J}}
\newcommand{\bExp}{\mathbb{E}}
\begin{document}
\maketitle
This homework has 150 points plus 30 bonus points available but, as always, homeworks are graded out of 100 points. Full credit will generally be awarded for a solution only if it is both correctly and efficiently presented using the techniques covered in the lecture and readings, and if the reasoning is properly explained. If you used software or simulations in solving a problem, be sure to include your code, simulation results, and/or worksheets documenting your work. If you score more
than 100 points, the extra points do count toward your homework total.
\section{Method of Averaging and Multiscale Simulation (80 points plus 30 bonus points)}
In class, we discussed the case of an overdamped physical system, which in suitably nondimensional form, reads
\begin{align}
\difd \bX (t) &= - \grad_{\bx} \phi (\bX (t),\bY (t)) \, \difd t + \sqrt{2 \sigma} \, \difd \bWx (t), \label{eq:fullxy}\\
\difd \bY (t) &= - \epsilon^{-1} \grad_{\by} \phi (\bX (t),\bY (t)) \, \difd t + \sqrt{\frac{2 \theta}{\epsilon}} \, \difd \bWy (t), \nonumber
\end{align}
where $ \phi (\bx,\by) $ is a nondimensionalized potential energy, $ \sigma $ is a noise coefficient, $ \theta $ is a nondimensionalized temperature, $ \epsilon \ll 1 $ is a small parameter describing the separation of time scales of the slow variables $ \bX (t) $ and fast variables $ \bY (t) $, and $ \bWx (t) $ and $ \bWy (t) $ are independent vector-valued Wiener processes. (The example in class had more general noise on the slow variables $ \bX (t) $ but to keep focus on the main issues in this problem, we make the simplest choice.)
We stated that the method of averaging indicates that the dynamics of $ \bX (t) $ can be well approximated for small $ \epsilon $ by the solution $ \bXavg (t) $ of the equation
\begin{equation}
\difd \bXavg (t) = - \grad_{\bx} \bphi (\bXavg (t)) \, \difd t + \sqrt{2 \sigma} \, \difd \bWx (t) \label{eq:xavg}
\end{equation}
where
\begin{equation*}
\bphi (\bx) \equiv - \theta \ln \int \expe^{-\phi (\bx,\by)/\theta} \, \difd \by
\end{equation*}
is a nondimensional free energy.
\begin{enumerate}
\item (\textbf{25 points}) Verify this approximation through choosing a potential $ \phi (\bx,\by) $ and simulating both the full system (\ref{eq:fullxy}) and the reduced system (\ref{eq:xavg}). It suffices to take each variable $ \bX $ and $ \bY $ to be one-dimensional if you like. Plot and compare some sample trajectories, then simulate an ensemble of trajectories for both the full and reduced system, and compare the resulting mean and covariance matrix of $ \bX (t) $. Choose a few values of $ \epsilon $ and comment on how well the reduced system approximates the full system as $ \epsilon $ is varied. \label{part:analytical}
\item (\textbf{10 points}) Comment on the relative computational costs of simulating the full system (\ref{eq:fullxy}) and reduced system (\ref{eq:xavg}). In particular, choose a relatively small value of $ \epsilon $ and see how large a time step you can choose for each system before your numerical method (which can just be the straightforward Euler-Marayama discussed in class) behaves badly (becomes very inaccurate or goes unstable).
\item (\textbf{25 points}) You may recall that the method of averaging more generally states that the reduced system (\ref{eq:xavg}) takes the form
\begin{equation}
\difd \bXavg = \aavg (\bXavg (t)) \, \difd t+ \sqrt{2 \theta} \, \difd \bWx (t) \label{eq:moa}
\end{equation}
where
\begin{equation*}
\aavg (\bx) \equiv \int - \grad_{\bx} \phi (\bx,\by) \pycsx (\by|\bx) \, \difd \by
\end{equation*}
and $ \pycsx (\by|\bx) $ is the stationary distribution for the stochastic differential system
\begin{equation*}
\difd \bYx = - \grad_{\by} \phi (\bx,\bYx (t)) \, \difd t + \sqrt{2 \theta} \, \difd \bWy (t).
\end{equation*}
with $ \bx $ treated as a fixed parameter.
We showed that
\begin{equation}
\pycsx (\by|\bx) = \frac{\expe^{-\phi (\bx,\by)/\theta}}{\int \expe^{-\phi (\bx,\by)/\theta} \, \difd \by}, \label{eq:pyboltz}
\end{equation}
from which the analytical reduction (\ref{eq:xavg}) follows after some elementary calculus manipulations. But suppose for some reason we don't realize that we have the exact solution (\ref{eq:pyboltz}) for $ \pycsx (\by|\bx) $.
Instead, as in multiscale computation, we estimate $ \pycsx (\by|\bx) $ from data obtained by
simulating the trajectories $ \bYx (t) $. Write a computational code that could simulate the reduced system (\ref{eq:xavg}) by estimating $ \pycsx (\by|\bx) $ in this spirit of multiscale simulations. \label{part:sample}
\item (\textbf{20 points}) Compare the simulations of $ \bXavg (t) $ from part (\ref{part:sample}) where you estimate $ \pycsx (\by|\bx) $ from simulations of $ \bYx (t) $ to the simulations from part (\ref{part:analytical}) where you are essentially using the analytical formula (\ref{eq:pyboltz}).
\item (\textbf{30 bonus points}) Consult Section 4 in the paper ``Hasselmann's Program Revisited: The Analysis of Stochasticity in Deterministic Climate Models'' by Ludwig Arnold (in the class reserves) for how to write an improved reduced system that also includes a stochastic correction term that accounts for fluctuations. (In that paper, no noise term appears explicitly in the original equations, but one can show that the only way their results need to be changed is that the noise term $ \sqrt{2 \sigma} \, \difd \bWx (t) $ simply is added on to their resulting equation for the reduced slow variable, as we did in Eq.~(\ref{eq:moa}) for the averaged equation). Write a code that incorporates this stochastic correction term, and compare the results of these simulations with those of Eq.~(\ref{eq:moa}) based on a direct method of averaging.
\end{enumerate}
\section{Brownian Motion in Tilted Periodic Potential (70 points)}
The dynamics for an overdamped particle moving through a periodic potential with an applied force and thermal fluctuations is given, in nondimensionalized form in one spatial dimension, by
\begin{equation}
\difd X (t) = \left[- \tphi^{\prime} (X (t)) + f\right] \, \difd t+ \sqrt{2 \theta} \difd W (t) \label{eq:periodsde}
\end{equation}
where $ \tphi (x) $ is the nondimensionalized potential with period $ 1 $, $ f $ is the nondimensionalized applied external force, and $ \theta $ is the nondimensionalized temperature. We stated in class that over a sufficiently long time, the effective drift of the particle is given by
\begin{equation}
\left\langle \frac{\difd X}{\difd t}\right\rangle = \lim_{t \rightarrow \infty} \left\langle \frac{X(t)}{t} \right\rangle
= \theta \frac{1 - \expe^{-f/\theta}}
{\int_0^1 \difd x \int_{x}^{x+1} \difd \xp \expe^{(\tphi (\xp) - \tphi (x) + f (x-\xp))/\theta}}
\label{eq:driftformula}
\end{equation}
\begin{enumerate}
\item (\textbf{20 points}) Show how to obtain this formula from the partial differential equation for the reduced flux and reduced probability density developed in class. Don't just verify the solution satisfies the partial differential equation -- show how you could have obtained this solution by systematically solving the partial differential equation. (It's fine to derive the solution in a different form and then show your solution is equivalent to (\ref{eq:driftformula})).
\item (\textbf{20 points}) Choose a particular nondimensionalized periodic potential $ \tphi (x) $ and simulate the stochastic differential equation (\ref{eq:periodsde}) for this choice of potential. Generate several trajectories and plot the average of $ (X(t)-X(0))/t $ over this statistical ensemble. How does this compare with the predicted rate of transport (\ref{eq:driftformula})? \label{part:specpot}
\item (\textbf{20 points}) Develop a simplified expression for the effective drift (\ref{eq:driftformula}) valid in the limit in which the thermal energy is small compared to the fluctuations in the potential $ \theta \ll 1 $ (and the applied external force is comparable or weaker than that of the potential, so that $ f $ may be taken as order unity.) You may either work with the general formula (\ref{eq:driftformula}) or with your particular choice of the nondimensionalized periodic potential from part~(\ref{part:specpot}).
\item (\textbf{10 points}) Compare your result for the effective drift from small $ \theta $ with the expression for the rate of escape from a potential at small temperature developed in class. If they are related, explain the connection; if not, explain the conceptual distinction between drifting across a periodic potential and escaping from a potential barrier, when in both cases the thermal energy is small compared to the barrier in potential energy.
\end{enumerate}
\section{Small Noise Asymptotics for Escape from Potential Well (95 points)}
In the (corrected) lecture notes, we derived an approximate expression for the rate of escape of a particle (undergoing Smoluchowski) dynamics from a potential well, when the ratio $ \epsilon = (k_B T)/A $ of the thermal energy to the potential barrier height $ A$ is small. In suitably nondimensional variables, the stochastic dynamics can be expressed as:
\begin{equation}
\difd X (t) = - \tphip (X(t)) \, \difd t + \sqrt{2 \epsilon} \difd W (t). \label{eq:escde}
\end{equation}
where $ \tphi (x) $ is the nondimensional potential energy, with minimum at $ 0 < m < 1 $ and local maxima at $ x= 0 $ and $ x=1 $.
The mean escape time $ f(x) = \bExp[T_1| X(0) = x]$, from the potential well given a starting location at $ x \in (0,1) $, is given by solving the boundary value problem:
\begin{equation}
\begin{aligned}
\epsilon \frac{\difd^2 f}{\difd x^2} - \tphip (x) \frac{\difd f}{\difd x} &= -1, \\
f(0) &= 0, \\
f(1) &= 0.
\end{aligned}
\label{eq:exitpde}
\end{equation}
Under the assumption that $ \tphi (0) > \tphi (1) $ (so the left peak is higher than the right peak), we obtained the asymptotic expression:
\begin{equation}
f(x) \sim \frac{\pi}{\sqrt{\left|\tphipp (1)\tphipp (m)\right|}} \exp \left(\frac{\tphi (1) - \tphi(m)}{\epsilon}\right) \text{ for }
\epsilon \ll 1, \epsilon \ll x \ll 1 - \epsilon. \label{eq:kramers}
\end{equation}
\begin{enumerate}
\item (\textbf{20 points}) Verify this formula (for some suitable choice of potential energy $ \tphi (x) $) through direct numerical simulations of Eq.~(\ref{eq:escde}). Consider a few choices of $ \epsilon $, and comment on what difficulties you find with the numerical simulations as $ \epsilon $ is made smaller.
\item (\textbf{15 points}) The mean escape time approximation~(\ref{eq:kramers}) is an exponentially large number (relative to the small parameter $ \epsilon $). This raises the question of whether the particle escapes the potential well by gradually crawling up the well, or whether the particle just suffers some fortuitous noise after some long time which moves it quickly from the bottom of the well to a peak.
Explore this question numerically by simulating trajectories until they escape the potential well, and examining/analyzing how they escaped. You can use a graphical and/or statistical approach, with more credit being given for more precise quantification and observations.
\item (\textbf{15 points}) Consider a sequence initial values of the particle position $ x $ that start near one of the boundaries. When the initial value is near the boundary, the asymptotic formula~(\ref{eq:kramers}) will begin to fail. Explore numerically on your example at what distance of the initial value from the boundary the asymptotic formula~(\ref{eq:kramers}) begins to break down, and examine whether there is a qualitative change in how trajectories are escaping from the potential well as the initial values cross this distance from the boundary.
\item (\textbf{20 points}) In the lecture, our derivation actually also excluded the case where $ |x - m| \sim O(\epsilon) $, but we said this restriction actually could be dropped. Offer some precise analytical arguments for why the formula~(\ref{eq:kramers}) still remains valid for the whole range of $ x $ stated in~(\ref{eq:kramers}), despite the exclusion of $ x $ near the minimum we made in the arguments in lecture. (That is, argue, with precision, why the restriction isn't really necessary.)
\item (\textbf{25 points}) If the right peak is higher than the left peak, then by a simple symmetry argument, the same formula~(\ref{eq:kramers}) applies with $ \tphi (1) $ and $ \tphipp (1) $ replaced by $ \tphi (0) $ and $ \tphipp (0) $. But when the peaks are the same height ($ \tphi (0) = \tphi (1) $), the formula changes more substantially, though the general argument still applies. Rework the asymptotic analysis involving Laplace's method to obtain the appropriate approximation for the escape time when $ \epsilon \ll 1 $ for the case where $ \tphi(0) = \tphi (1) $. You may either modify the approach taken in lecture, or apply Laplace's method directly to the exact solution of Eq.~(\ref{eq:exitpde}). But I am not interested in seeing a recapitulation of the various heuristic arguments you can find easily in various physical references. The credit for this part is for conducting an asymptotic analysis that is at least as systematic as the one adoped in the lecture.
\end{enumerate}
\end{document}