CoCalc Shared Filesmiw / discretization1.tex
Author: Bill Page
Views : 13
Description: Discretization of the Quantum Dynamical Equations of Motion - One Dimension
1\documentclass[notitlepage]{article}
2\usepackage{titling}
3\usepackage[utf8]{inputenc}
4\usepackage{titling}
5\usepackage{fullpage}
6\usepackage{mathtools}
7\usepackage{etoolbox}
8\cslet{lrbox*}\lrbox
9\expandafter\patchcmd\csname lrbox*\endcsname{\setbox}{\global\setbox}{}{}
10%\expandafter\show\csname lrbox*\endcsname % uncomment to see if it has worked
11\cslet{endlrbox*}\endlrbox
12\usepackage{sagetex}
13\usepackage{url}
14%\usepackage{amsmath}
15\usepackage{cite}
16\usepackage{xcolor}
17\usepackage{breqn}
18\usepackage{collect}
19\setlength{\parskip}{1ex}
20\usepackage[titletoc,title]{appendix}
21\usepackage{float}
22\usepackage{hyperref}
23\hypersetup{
26    citecolor={blue!50!black},
27    urlcolor={blue!80!black}
28}
29
30\title{Discretization of the Quantum Dynamical Equations of Motion - One Dimension}
31\author{Bill Page \small{(\href{mailto:bill.page@newsynthesis.org}{bill.page@newsynthesis.org})}}
32\date{Draft \today}
33
34\begin{document}
35
36\maketitle
37\thispagestyle{empty}
38
39\begin{abstract}
40\textbf{From Poirier’s Bohmian Mechanics without Wavefunctions, to Hall’s Many Interacting Worlds
41in One Dimension}
42
43In this note we show that Hall's toy'' model of quantum mechanics as many interacting interacting classical worlds in one dimension follows directly from the approximation of Poirier's quantum mechanical equations of motion as difference equations with a judicious choice of difference operators. These calculations are a preliminary exercise for the case of more than one dimension. The main emphasis is methodological with the main applications to follow in later articles.
44
45We provide a numerical simulation that reproduces the results given previously by Poirier\cite{poirier2008}.
46\end{abstract}
47
48%\setcounter{section}{-1}
49
50\section{Introduction}
51
52In a general sense the quantum mechanical equations of motion originates with the work of Erwin Madelung\cite{madelung1927} on quantum hydrodynamics. David Bohm\cite{bohm1993} rediscovered and extended an interpretation of quantum mechanics originating with Louis de Broglie\cite{debroglie1960}, aka. Bohmian mechanics, which the notion of particles guided'' by the wavefunction. Peter Holland\cite{holland1995} and co-workers have continued to extend the de Broglie-Bohm theory in this direction, ultimately showing that it is possible to dispense with the wavefunction all together in favor of a continuous ensemble of quantum trajectories. Besides the conceptual metaphysical issues, this has lead to some new more powerful calculation methods that have become important for practical application of quantum mechanics to physical chemistry (e.g. Wyatt \cite{wyatt2005}). Poirier continues in this tradition while Hall, Deckert, and Wiseman\cite{hall2014quantum} emphasize the interpretation of the quantum trajectories as interaction between parallel worlds reminiscent but differing from Hugh Everett's many-worlds interpretation of quantum mechanics.
53
54In this article the Hamiltonian formulaion of these equations of motion is important. Conservation of energy is an important indicator of the quality of the numerical integration. The Störmer–Verlet method implemented in this article is well known and preserves the symplectic structure on the system phase space.
55
56From the point of view of Poirier\cite{schiff2012communication} this simulation is viewed as an approximation to the quantum dynamical field equations. One should think of this in hydrodynamic terms. But from the point of view of Hall, et al. these are discrete trajectories of the \textbf{same} particle in many different parallel'' worlds - but just one particle!
57
58In the expressions which follow $\mathbf{C}$ is the co-ordinate which indexes these worlds. One should think of this as the same particle in many different worlds interacting with the copy of itself in other worlds. It is a an important aspect of the equations of motion that trajectories of these particles never cross. The co-ordinate function $\mathbf{C}$ is assumed to be uniformized''\footnote{ Schiff and Poirier\cite{schiff2012communication}} so that discretization of $\mathbf{C}$  assigns an equal density'' to the particle in each world.
59
60In the example simulation included in this article the starting distribution of the particle in these different worlds is a Gaussian distribution. But the distribution evolves as it interacts with the classical potential. Some trajectories are scattered from the potential while other tunnel through. From the point of view of Poirier, an important goal of this simulation is to provide an efficient calculation of the scattering and tunneling amplitudes. But from the point of view of Hall this is a just purely "classical" mechanical problem with a rather unusual form of interaction. It is the form of this interaction between worlds that leads to statistics that approximate quantum mechanics in the limit of a very large number of interacting worlds.
61
62\newcommand{\J}{\begin{vmatrix}\scriptstyle{J}\end{vmatrix}}
63
64\definecollection{prelim}
65\begin{sagesilent}[prelim]
66import numpy
67from numpy import array,concatenate,vectorize,isnan
68rdf=numpy.array([RDF.pi()]).dtype # ensure compatible types
69import mpmath
70from mpmath import erfinv
71
72def fun(*y): # fun(x,y)(f)(x,y)=f
73    cache={}
74    def upd(f,x,y): # memoized
75        try: return cache[x]
77            try: cache[x] = f.subs(dict([[y[i],x[i]]
78                                        for i in range(len(y))]))
79            except AttributeError: # ask forgiveness
80                def sub(f): return f.subs(dict([[y[i],x[i]]
81                                        for i in range(len(y))]))
82                cache[x] = map(sub,f)
83        return cache[x]
84    return lambda f: lambda *x: upd(f,x,y)
85def argscript(self, *args):
86    return "%s_{%s}"%(self.name(),','.join(map(repr, args)))
87def functions(*names):
88    return map(lambda nam:function(nam, print_latex_func=argscript),names)
89def detJscript(self, *args):
90    return '\J_{%s}'% (','.join(map(repr, args)))
91J__= function('J_') # determinant of J
92J_ = function('J_', print_latex_func=detJscript)
93def squared(x): # x^2
94    try: return x.dot(x) # numpy array
95    except: return power(x,2)
96\end{sagesilent}
97
98\definecollection{shift}
99\begin{sagesilent}[shift]
100def shift(x,n):
101    try: # deduce slicing
102        start = (x.__array_interface__['data'][0] -
103                 x.base.__array_interface__['data'][0])/x.itemsize
104        stop = start + x.shape[0]
105        if start+n<0 or stop+n>x.base.shape[0]:
106            raise RuntimeError('bounds')
107        return x.base[start+n:stop+n]
108    except AttributeError:
109        raise RuntimeError('expecting view')
110\end{sagesilent}
111
112
113In this document we have used the computer algebra system
114Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of
115SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}
116and
117mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}
118for both symbolic and numeric calculations. This document and the source code for this work can be found
119online\footnote{\LaTeX\  \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.tex}{source} \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.pdf}{pdf} \href{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/miw/\jobname.sagews}{SageWorksheet}. The online versions may be more up to date.}. The inset bold text below and many of the formulas in this paper are input and output produced directly by the Sage system. For more details please refer to Appendix \ref{SageAppendix}.
120
121For example we define some special symbols for Planck's constant and mass:
122\begin{sageblock}
123hbar = var('hbar',latex_name=r'\hbar')
124m = var('m')
125\end{sageblock}
126
127\definecollection{fancy}
128\definecollection{prelim1}
129\begin{collect}{prelim1}{\noindent}{\par}
130The following is a work-round for a bug\footnote{Some commands for manipulating and simplifying expressions do not respect the \textbf{\rm print\_latex\_func} option. \url{https://trac.sagemath.org/ticket/19151}} in Sage.
131\end{collect}
132\definecollection{prelim2}
133\begin{sagecollect}{prelim2}
134def fix(it):
135    ret = it.subs(dict([[function('x')(n+i), X[0](n+i)] for i in range(-2,3)] +
136                       [[J__(n+i), J_(n+i)] for i in range(-2,3)]))
137    return ret
138def factor(x):
139    try: return fix(x.factor())
140    except: return vector(map(lambda y:fix(y.factor()),x))
141\end{sagecollect}
142
143\section{One Dimensional System}
144
145We show below that in one dimension Schiff and Poirier's quantum dynamical equation of motion,  \cite{schiff2012communication} eqs. (18-20),
146\begin{equation} \label{eq18}
147m \ddot{x^i}+\frac{\partial V(x)}{\partial x^i}- \frac{\hbar^2}{4 m}\frac{\partial}{\partial C^m}\left(K^k_i K^m_j \frac{\partial^2 K^l_j}{\partial C^k \partial C^l}\right)=0
148\end{equation}
149that follows from a quantum potential of the form
150\begin{equation} \label{eq20}
151Q = -\frac{\hbar^2}{4 m} \left(K^k_j \frac{\partial^2 K^l_j}{\partial C^l \partial C^k} + \frac{1}{2}\frac{\partial K^l_j}{\partial C^l}\frac{\partial K^k_j}{\partial C^k} \right)
152\end{equation}
153
154when discretized by replacing derivatives with difference operators, is equivalent to Hall’s \cite{hall2014quantum} toy'' expression eq. (24) for a finite number of classically interacting parallel'' worlds.
155
156\begin{align}
157r_N(x_n;X) &= \frac{\hbar^2}{4 m}[\sigma_{n+1}(X)-\sigma_n(X)] \label{eq24} \\
158\sigma_n(X) &= \frac{1}{(x_n-x_{n-1})^2} \left[ \frac{1}{x_{n+1}-x_n} - \frac{2}{x_n-x_{n-1}} + \frac{1}{x_{n-1}-x_{n-2}}\right]  \nonumber
159\end{align}
160
161By discretized we mean that the uniformizing co-ordinate'' $\mathbf{C}$ will be replaced with discrete values representing a finite number of worlds''. Derivatives with respect to $\mathbf{C}$ will be replaced with difference operators in the expressions for the quantum potential and quantum force. $\mathbf{C}$ becomes a discrete index $n$.
162
163In the following we write most expressions in their general vector form even though in one dimension considerable simplification is possible to connect with related articles where we use these expressions in more than one dimension.
164
165Let $E =[x_{n}]$ represent the configuration space of one particle in 1-D space with world co-ordinate $n$. In Sage this can be represented as:
166\begin{sageblock}
167X = functions('x')
168d = len(X)
169C = (var('n'),)
170E = tuple([x(*C) for x in X])
171\end{sageblock}
172\begin{equation} \sage{E} \end{equation}
173
174\section{Forward/Backward Difference Approximation}
175
176In this approximation to the derivative each world interacts only locally but with the
177\emph{two} nearest neighbors to the left and right.
178
179The forward and backward difference operators
180\begin{sageblock}
181def DM(x,n):return x-x.subs(n==n-1)
182def DP(x,n):return x.subs(n==n+1)-x
183\end{sageblock}
184 approximate derivatives.
185
186To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes a $1x1$ matrix
187\begin{sageblock}
188J = fun(n) (matrix([[DM(E[j],C[i]) for i in range(d)] for j in range(d)]))
189\end{sageblock}
190\begin{equation} \sage{J(n)} \end{equation}
191It will be convenient for our calculations to represent the inverse of the Jacobian in terms of its adjugate and delay evaluation of the determinant until needed. The adjugate matrix'' is just
192\begin{sageblock}
194\end{sageblock}
195\begin{equation} \sage{A(n)} \end{equation}
196so the inverse matrix'' is
197\begin{sageblock}
198K = fun(n) (A(n)/J_(n))
199\end{sageblock}
200\begin{equation} \sage{K(n)} \end{equation}
201where the symbol $\J_n$ represents the determinant. We will eventually need to replace the symbolic determinants with then explicit expressions
202\begin{sageblock}
203_J = fun(n) (factor(J(n).determinant()))
204dets = {J_(n+i):_J(n+i) for i in range(-2,3)}
205\end{sageblock}
206\begin{equation} \sage{_J(n)} \end{equation}
207It is $\sage{factor(J(n)^(-1))==K(n).subs(J_(n)==_J(n))}$ that $J(n)^{-1}=K(n) \text{ where }{\sage{J_(n)}=\operatorname{_J}(n)}$.
208
209\subsection{Quantum Potential}
210
211The Schiff and Poirier expression for the quantum potential \eqref{eq20} can be approximated in terms of difference operations. Let $Q(n)=Q_2(n)+Q_3(n)$ be the 2nd and 3rd order terms, respectively.
212\begin{sageblock}
213Q = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
214    (K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k] ) +
215        1/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]))
216            for k in range(d)) for j in range(d)) for l in range(d)) )
217\end{sageblock}
218\begin{dmath}
219\sage{ expand(factor(Q(n))) }
220\end{dmath}
221
222\begin{sageblock}
223Q2 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
224    (1/2 * DP(K(n)[l,j],C[l]) * DP(K(n)[k,j],C[k]) )
225        for k in range(d)) for j in range(d)) for l in range(d)) )
226\end{sageblock}
227\begin{dmath}
228\sage{ expand(factor(Q2(n))) }
229\end{dmath}
230
231\begin{sageblock}
232Q3 = fun(n) ( -(hbar^2/4/m) * sum( sum( sum(
233    (K(n)[k,j] * DM( DP( K(n)[l,j],C[l] ),C[k]) )
234        for k in range(d)) for j in range(d)) for l in range(d)) )
235\end{sageblock}
236\begin{dmath}
237\sage{ expand(factor(Q3(n))) }
238\end{dmath}
239
240In the discretized expressions it will be convenient to replace the $X$ co-ordinates with scalar quantities. The forward and backward difference approximations will then require $2x5x5 = 50$ scalar quantities.
241
242\begin{sageblock}
243xyk = [X[i](n+j)
244          for j in range(-2,3) for i in range(d)]
245xyv = [var(X[i].name()+str(2+j))
246          for j in range(-2,3) for i in range(d)]
247xykv = dict(zip(xyk,xyv))
248xyvk = dict(zip(xyv,xyk))
249xn = X[0](n).subs(xykv)
250\end{sageblock}
251\begin{dmath}
252\sage{(n==xn,)}
253\end{dmath}
254
255\subsection{Quantum Force}
256
257The Schiff and Poirier expression for the quantum force is the last term on the left of \eqref{eq18}. It can be rewritten as terms of difference operations as follows
258\begin{sageblock}
259R = fun(n) (vector([ (hbar^2/4/m) * sum( sum( sum( sum(
260   DP( K(n)[k,i] * K(n)[p,j] * DP( DM( K(n)[l,j],C[l] ),C[k] ),C[p] )
261       for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
262           for i in range(d)]) )
263\end{sageblock}
264\begin{dmath}
265\sage{ expand(R(n)[0]).subs(dets) }
266\end{dmath}
267
268We can compute the discretized quantum force as the negative gradient of the discretized quantum potential symbolically. We would like to compare this with the discretization of Poirier's expression for the force.
269
270To get the corresponding expression for the force we sum the potential symbolically over a sufficiently large area and then differentiate by $x(n)$. Since the potential is local, involving only nearby worlds, the result is independent of the size of the area so long as the area is large enough to include all the neighbors of a given world at $n$.
271
272\noindent
274\begin{sageblock}
275r = fun(n) (vector([ -sum(Q(n+i) for i in range(-2,3)
276                  ).subs(dets).subs(xykv).diff(xn).expand()
277                      for k in range(d)]))
278\end{sageblock}
279\begin{dmath}
280\sage{ expand(r(n).subs(xyvk)) }
281\end{dmath}
282with $r_2$ and $r_3$ the gradients of the 2nd and 3rd order terms of the potential, separately.
283\begin{sageblock}
284r2 = fun(n) (vector([ -sum(Q2(n+i) for i in range(-2,3)
285          ).subs(dets).subs(xykv).diff(xn).expand()
286              for k in range(d)]))
287\end{sageblock}
288\begin{dmath}
289\sage{ expand(r2(n).subs(xyvk)) }
290\end{dmath}
291
292\begin{sageblock}
293r3 = fun(n) (vector([ -sum([Q3(n+i) for i in range(-2,3)]
294         ).subs(dets).subs(xykv).diff(xn).expand()
295             for k in range(d)]))
296\end{sageblock}
297\begin{dmath}
298\sage{ expand(r3(n).subs(xyvk)) }
299\end{dmath}
300
301Of course $r2(n)+r3(n)=r(n)$ is $\sage{r2(n)+r3(n)==r(n)}$.
302
303Surprisingly, the force term from the 3rd order term is exactly -2 times the force term obtained from the 2nd order term so the difference between the force from just the 2nd order term so the force from the full expression for the quantum potential differs from the 2nd term by only a sign.
304
305It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.
306
307This result is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression eq. (18) for the quantum force. We will see later that this is not the case in more than one dimension.
308
309It is $\sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))}$ that $r(n)==R(n)$
310
311Hall's equation \eqref{eq24} can be written
312
313\begin{sageblock}
314sigma = fun(n) ( (1/(X[0](n)-X[0](n-1))^2) *
315    (1/(X[0](n+1)-X[0](n)) - 2/(X[0](n)-X[0](n-1)) + 1/(X[0](n-1)-X[0](n-2))) )
316r_hall = fun(n) (hbar^2/4/m * (sigma(n+1) - sigma(n)))
317\end{sageblock}
318\begin{dmath}
319\sage{ expand(r_hall(n)) }
320\end{dmath}
321
322It is $\sage{bool(r(n)[0]==r_hall(n).subs(xykv))}$ that $rr=r_{hall}$.
323
324Thus Hall's model of quantum mechanics as many interacting interacting classical worlds in one dimension follows directly from the approximation of Poirier's quantum mechanical equations of motion as difference equations with a certain choice of operators. Obviously this is not the only way to approximate Poirier's equations as difference equations.
325
326Save these calculations for later.
327\begin{sageblock}
328Q_hall = Q(n).subs(dets)
329R_hall = R(n).subs(dets)
330xyv_hall = xyv
331xykv_hall = xykv
332\end{sageblock}
333
334\section{Symmetric Difference Approximation}
335
336The central difference approximation to this derivative is
337\begin{sageblock}
338def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2
339\end{sageblock}
340
341In this approximation we will see that each world interacts with its four neighbors to the left and right.  In fact we may continue this process of approximation to higher orders with more neighbors but in all these Hall-like models this non-locality in phase space is always bounded. Perhaps most interesting form the point of view of efficient computation is the lowest order approximation given previously. But the surprise is that all of these approximations gives quantum mechanics in the limit of a very large number of worlds (i.e. a continuum) and are good approximations even for a very small number of worlds.
342
343Repeating the computations above using a symmetric difference approximation to the derivatives we again find that the quantum force computed as the gradient of the discretized quantum potential is identical to the discretization of Poirier's expression for the force.
344
345\begin{sageblock}
346J = fun(n) (matrix([[D(E[j],C[i]) for i in range(d)] for j in range(d)]))
347\end{sageblock}
348\begin{equation} \sage{J(n)} \end{equation}
349\begin{sageblock}
351K = fun(n) (A(n)/J_(n))
352\end{sageblock}
353\begin{equation} \sage{K(n)} \end{equation}
354\begin{sageblock}
355_J = fun(n) (fix(J(n).determinant().expand()))
356dets = {J_(n+i):_J(n+i) for i in range(-6,7)}
357\end{sageblock}
358\begin{equation} \sage{_J(n)} \end{equation}
359
360\begin{sageblock}
361Q = fun(n) ( -(hbar^2/4/m) *
362    sum( sum( sum(
363        (K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ) +
364            1/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))
365                for k in range(d)) for j in range(d)) for l in range(d)) )
366Q2 = fun(n) ( -(hbar^2/4/m) *
367    sum( sum( sum( (1/2 * D(K(n)[l,j],C[l]) * D(K(n)[k,j],C[k]))
368        for k in range(d)) for j in range(d)) for l in range(d)) )
369Q3 = fun(n) ( -(hbar^2/4/m) *
370    sum( sum( sum( (K(n)[k,j] * D( D( K(n)[l,j],C[l] ),C[k] ))
371        for k in range(d)) for j in range(d)) for l in range(d)) )
372\end{sageblock}
373\begin{dmath}
374\sage{ expand(factor(Q(n))) }
375\end{dmath}
376\begin{sageblock}
377R = fun(n) (vector([ (hbar^2/4/m)*sum( sum( sum( sum(
378    D( K(n)[k,i] * K(n)[p,j] * D( D( K(n)[l,j],C[l] ),C[k] ),C[p] )
379        for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
380            for i in range(d)]) )
381\end{sageblock}
382\begin{dmath}
383\sage{ expand(R(n)[0].subs(dets)) }
384\end{dmath}
385
386\begin{sageblock}
387xyk = [X[i](n+j) for j in range(-4,5) for i in range(d)]
388xyv = [var(X[i].name()+str(5+j)) for j in range(-4,5) for i in range(d)]
389xykv = dict(zip(xyk,xyv))
390xn = X[0](n).subs(xykv)
391\end{sageblock}
392
393\begin{dmath}
394\sage{(n==xn)}
395\end{dmath}
396
397The neighborhood for the symmetric difference operators is larger than previous case.
398\begin{sageblock}
399r = fun(n) (vector([-sum(Q(n+i) for i in range(-3,4)
400            ).subs(dets).subs(xykv).diff(xn).expand()
401                for k in range(d)]))
402r2 = fun(n) (vector([-sum(Q2(n+i) for i in range(-3,4)
403            ).subs(dets).subs(xykv).diff(xn).expand()
404                for k in range(d)]))
405r3 = fun(n) (vector([-sum([Q3(n+i) for i in range(-3,4)]
406            ).subs(dets).subs(xykv).diff(xn).expand()
407               for k in range(d)]))
408\end{sageblock}
409
410It is $\sage{r3(n)==-2*r2(n)}$ that $r3(n)= -2\ r2(n)$.
411
412It is $\sage{bool(factor(r(n))==factor(R(n).subs(dets).subs(xykv)))}$ that $r(n)=R(n)$
413
414This expression is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression \eqref{eq18} for the quantum force.
415
416Save these calculations for later.
417\begin{sageblock}
418Q_poirier = Q(n).subs(dets)
419R_poirier = R(n).subs(dets)
420xyv_poirer = xyv
421xykv_poirier = xykv
422\end{sageblock}
423
424
425\section{Simulation}
426
427Using the expressions derived above we can generate efficient numerical functions which we use later in determining the motion of a particle in many worlds simultaneously.
428
429We will use Hall's difference approximation to Poirier's equations above for a numerical simulation of quantum motion in a small number of parallel'' worlds.
430
431For efficiency we expression the following operations in  vectorized'' form. First we need the difference operator
432\begin{sageblock}
433def dM(x,n):return shift(x,n) - shift(x,n-1)
434\end{sageblock}
435
436
437\subsection{Classical Potential}
438Poirier's problem involves scattering a proton-mass particle from an Eckart potential slightly less than the particles incident energy.
439\begin{sageblock}
440alf = 2.5 ; x0=0 ; B = 0.0024
441V_(x) = B /cosh(alf *(x-x0))^2
442#V_(x) = 2/cosh(2*x)^2
443#V_(x) = 0
444V = vectorize(fast_callable( V_(x), vars=[x], domain=RDF))
445_ = plot(V,(-3,3),axes_labels=['distance','energy'],axes_labels_size=1,
446         legend_label='potential')
447\end{sageblock}
448\begin{figure}[H] \centering
449\sageplot[width=.9\textwidth]{_}
450\caption{Classical Potential}
451\end{figure}
452
453\subsection{Classical Force}
454\begin{sageblock}
455F = vectorize(fast_callable( -diff(V_(x),x), vars=[x], domain=RDF))
456_ =plot(F,(-3,3),axes_labels=['distance','force'],axes_labels_size=1,
457        legend_label='force')
458\end{sageblock}
459\begin{figure}[H] \centering
460\sageplot[width=.9\textwidth]{_}
461\caption{Classical Force}
462\end{figure}
463
464\subsection{Quantum Potential}
465
466\begin{sagesilent}[fancy]
467def qq(hbar,m,x):
468    dx0  = dM(x,0)
469    dxp1 = dM(x,+1)
470    dxm1 = dM(x,-1)
471    return hbar^2/4/m*(-1/2/dxp1^2 + 3/2/dx0^2 - 1/dxm1/dx0)
472\end{sagesilent}
473\includesagecollection{fancy}
474
475\subsection{Quantum Force}
476In general the force is a vector but in the following we treat it as a single scalar quantity.
477
478Vectorized Quantum Force
479\begin{sagesilent}[fancy]
480def ff(hbar,m,x):
481    dx0  = dM(x,0)
482    dxp1 = dM(x,+1)
483    dxp2 = dM(x,+2)
484    dxm1 = dM(x,-1)
485    return hbar^2/2/m*(-1/dxp1^3 + 1/2/dxp2/dxp1^2 + 1/dx0^3 - 1/2/dxp1/dx0^2 - 1/2/dxm1/dx0^2 + 1/2/dxp1^2/dx0)
486\end{sagesilent}]
487\includesagecollection{fancy}
488
489\subsection{Initial Data}
490
491\subsubsection{Spatial Distribution}
492
493In order to compute the quantum force on the particles in world $n$ we introduce four widely separated fictitious particles in worlds to the left and right to represent open boundary conditions.
494\begin{sagesilent}[fancy]
495def boundary(x): # returns x as a view over a base with boundary values
496    return concatenate((
497        array([-1e18,-1e14,-1e10,-1e6]),
498        x,
499        array([1e6,1e10,1e14,1e18])))[4:-4]
500\end{sagesilent}
501\includesagecollection{fancy}
502
503For example consider the spatial distribution of a Gaussian ensemble of one particle in  $N$ parallel'' worlds. Everything is in atomic units below, so $\hbar=1$.
504\begin{sagesilent}[fancy]
505N = 2000
506var('alf,x0,m,v,x1')
507model1={m:RDF(2000),alf:RDF(0.70),v:RDF(0.00164317),x0:RDF(-7.0)}
508W(x)=abs((alf/pi)^(1/4)*exp(-(alf/2)*(x-x0)^2)*exp(I*m*v*x))^2
509\end{sagesilent}
510\includesagecollection{fancy}
511\begin{equation}
512\sage{W(x)}
513\end{equation}
514We may compute a uniform mapping into this distribution by inverting the cummulative distribution.
515\begin{sagesilent}[fancy]
516assume(alf>0)
517erf0=integrate(W(x),x,-oo,x1).subs(model1)
518y1=var('y1');inverse_erf=function('inverse_erf')
519erfinv0=solve(
520    integrate(W(x),x,-oo,x1).subs(model1)==y1,x1)[0].rhs()
521def erfinv1(x): return erfinv0.subs(y1==x).substitute_function(inverse_erf,erfinv)
522\end{sagesilent}
523\includesagecollection{fancy}
524After adding the open boundary conditions, we can check the accuracy by comparing the discrete representation to the original distribution.
525\begin{sageblock}
526p0 = array([erfinv1((i+0.001)/(N-1+0.002)) for i in range(N)], dtype=rdf)
527b0 = boundary(p0)
528_ = (list_plot([[b0[i],2/N/(b0[i+1]-b0[i-1])] for i in range(4,len(b0)-4) ],
529              axes_labels=['location','density'], color='red',
530              legend_label='discrete')+
531     plot(W(x).subs(model1),(x,-15,5),
532              legend_label='continuous'))
533\end{sageblock}
534\begin{figure}[H] \centering
535\sageplot[width=.9\textwidth]{_}
536\caption{Initial Spatial Density}
537\end{figure}
538
539\subsubsection{Distribution of Velocity}
540
541The particles have an identical initial velocity in each world of $\sage{v.subs(model1)}$.
542\begin{sageblock}
543v0=array([v.subs(model1) for i in range(N)],dtype=rdf)
544\end{sageblock}
545
546\subsection{Energy of the Ensemble}
547\noindent
548Total classical potential
549\begin{sageblock}
550def PU(m,x): return sum(V(i) for i in x)
551\end{sageblock}
552\noindent
553Total kinetic energy
554\begin{sageblock}
555def KU(m,v): return RDF(0.5)*m*(v.dot(v))
556\end{sageblock}
557\noindent
558Total quantum potential
559\begin{sageblock}
560def QU(m,hbar,x): return qq(m,hbar,x).sum()
561\end{sageblock}
562
563\subsection{Acceleration}
564
565Acceleration is the combined effect of both classical and quantum forces.
566\begin{sageblock}
567def A(hbar,m,x): return (F(x) + ff(hbar,m,x))/m
568\end{sageblock}
569
570\subsection{Step Size Controller}
571
572Continuous integrating step size controller \cite{hairer2005}
573\begin{sagesilent}[fancy]
574def G(a,v):
575    alpha = RDF(-100.0) # sensitivity
576    eps = RDF(10.0^-10)
577    return -alpha*(a.dot(v))/max(v.dot(v),eps)
578\end{sagesilent}
579\includesagecollection{fancy}
580
581\subsection{Integration}
582Integrate using the Störmer–Verlet algorithm with adaptive step size.
583
584\begin{sagesilent}[fancy]
585t_end = 11600; t_samples = 100; x_samples = 100
586t_sample = t_end/t_samples
587x_sample = max(1,int(N)/int(x_samples))
588x_start = max(0,int(N-x_samples*x_sample)/int(2))
589XS = range(x_start,N-x_start-1,x_sample)
590ds0 = RDF(0.007) # initial step size
591dsn = 1 # initial step size divider
592xb = boundary(p0)
593while True:
594    try: # step size
595        ds = ds0/dsn
596        rho = RDF(1.0) # step density
597        dt = ds/rho
598        m0 = RDF(2000.0)
599        hbar0 = RDF(1.0)
600        # initial values
601        x = boundary(p0); v = v0.copy(); a = A(hbar0,m0,x)
602        rho = rho - RDF(0.5)*G(a,v)*dt
603        t = 0;  T = [t]
604        XX = {t:x[XS].copy()}; XV = {t:v[XS].copy()}
605        KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)
606        E1 = KK+PP+QQ  # total energy
607        Rho = {t:rho}; TK = {t:KK}; TP = {t:PP}; TQ = {t:QQ}; E = {t:E1}
608        try:
609            while t<t_end:
610                print "t =", t, "rho =", rho
611                t1 = t + t_sample
612                while t<t1:
613                    rho = rho + G(a,v)*dt
614                    dt = ds/rho; t += dt
615                    x += v*dt + 0.5*a*dt^2
616                    v += 0.5*a*dt; a = A(hbar0,m0,x); v += 0.5*a*dt
617                    if isnan(rho) or rho>100 or rho<0.01:
618                        raise ValueError(
619                            "Step control failed at %s. rho=%s"%(t,rho))
620                # Check no-crossing
621                if (x[1:]-x[:N-1]).min()<0:
622                    raise ValueError("crossing")
623                # Check total energy conservation
624                KK = KU(m0,v); PP = PU(m0,x); QQ = QU(hbar0,m0,x)
625                E2 = KK+PP+QQ
626                if abs(E1-E2)>0.01:
627                    raise ValueError(
628                        "Energy conservation bound failed at %s. Delta E:"\
629                        "|%s| > 0.01."%(t,E1-E2))
630                E1 = E2
631                T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()
632                Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1
633        except KeyboardInterrupt:
634            print "Interrupted at %s ..."%(t)
635            T += [t]; XX[t] = x[XS].copy(); XV[t] = v[XS].copy()
636            Rho[t] = rho; TK[t] = KK; TP[t] = PP; TQ[t] = QQ; E[t] = E1
637        tmax = t
638        break
639    except ValueError as msg:
640        print msg
641        dsn = dsn + 1
642        print "Trying a shorter initial step size: %s."%(ds0/dsn)
643        continue
644print "t =", tmax, "rho =",rho
645\end{sagesilent}
646\includesagecollection{fancy}
647
648The adaptation of the integration stepsize during the simulation is given by \textbf{\rm rho}. A value greater than 1 indicates a smaller stepsize.
649\begin{sageblock}
650_ = line2d([[t,Rho[t]] for t in T],axes_labels=['time','rho'])
651\end{sageblock}
652\begin{figure}[H] \centering
653\sageplot[width=.9\textwidth]{_}
654\caption{Step Density}
655\end{figure}
656
657Energy is exchanged between classical, quantum mechanical and kinetic forms due to the interaction of the particle with the potential barrier (in each world) as well as the interaction between worlds.
658\begin{sageblock}
659_ = (line2d([[t,TK[t]] for t in T], axes_labels=['time','energy'],
660        legend_label='kinetic', color='green') +
661    line2d([[t,TQ[t]] for t in T], axes_labels=['time','energy'],
662        legend_label='quantum potential',color='red') +
663    line2d([[t,TP[t]] for t in T], axes_labels=['time','energy'],
664        legend_label='classical potential',color='blue'))
665\end{sageblock}
666\begin{figure}[H] \centering
667\sageplot[width=.9\textwidth]{_}
668\caption{Quantum Potential, Classical Potential and Kinetic Energy}
669\end{figure}
670
671Conservation of energy is an important indicator of the quality of the numerical integration.
672Note in particular how the changes in total energy over time remain well-bounded.
673\begin{sageblock}
674_ = line2d([[t,E[t]] for t in T],axes_labels=['time','energy'],
675        legend_label='total energy',color='black')
676\end{sageblock}
677\begin{figure}[H] \centering
678\sageplot[width=.9\textwidth]{_}
679\caption{Energy Conservation}
680\end{figure}
681
682The trajectory of the particle in each world is shown below.  Notice how the trajectories do not cross, yet in some worlds the particle is scattered from the barrier while in others the particle tunnels across the barrier pushed'' by its counterpart in other worlds.
683\begin{sageblock}
684_ = (line([[t,XX[t].mean()] for t in T], color='black',
685          thickness=2, legend_label="average",
686          axes_labels=['time','location'], axes_labels_size=1)+
687     sum(line([[t,XX[t][i]] for t in T], color=hue((0.5*i)/len(XS)), alpha=0.5,
688              axes_labels=['time','location'], axes_labels_size=1)
689                  for i in range(len(XS))))
690\end{sageblock}
691\begin{figure}[H] \centering
692\sageplot[width=.9\textwidth]{_}
693\caption{Trajectories}
694\end{figure}
695
696The distribution at the end of the simulation shows the reflected and transmitted "wave packets" separating from each other.
697\begin{sageblock}
698_ = list_plot([[x[i-1],2/N/(x[i+1]-x[i-1])] for i in range(1,N-1) ],
699              axes_labels=['location','density'])
700\end{sageblock}
701\begin{figure}[H] \centering
702\sageplot[width=.9\textwidth]{_}
703\caption{Final Spatial Density}
704\end{figure}
705
706Transmission probability
707\begin{sageblock}
708_ = line([[t,(XV[t]>0).mean()] for t in T ],
709#_ = line([[t,(XX[t]>0).mean()] for t in T ],
710         axes_labels=['time','transmission'])
711\end{sageblock}
712\begin{figure}[H] \centering
713\sageplot[width=.9\textwidth]{_}
714\caption{Evolution of Spatial Density}
715\end{figure}
716Final $\sage{(v>0).mean()}$.
717
718\newpage
719\begin{appendices}
720\section{Sage Preliminaries}\label{SageAppendix}
721
722It is convenient to define a compact notation for functions:
723\includesagecollection{prelim}
724\includecollection{prelim1}
725\includesagecollection{prelim2}
726This routine shifts a numpy array view left or right by $n$ cells.
727\includesagecollection{shift}
728
729\end{appendices}
730\newpage
731\listoffigures
732\bibliographystyle{../bibstyles/utphys}
733\noindent
734\bibliography{../bibliography}
735
736\end{document}
737
738
739%sagemathcloud={"latex_command":"pdflatex -synctex=1 -interact=nonstopmode 'poirier-problem1.tex'"}