CoCalc Public Filesdiscretization.tex
Author: Bill Page
Views : 162
Compute Environment: Ubuntu 18.04 (Deprecated)
1\documentclass[notitlepage]{article}
2\usepackage{titling}
3\usepackage[utf8]{inputenc}
4\usepackage{titling}
5\usepackage{fullpage}
6\usepackage{mathtools}
7\usepackage{sagetex}
8\usepackage{url}
9%\usepackage{amsmath}
10\usepackage{cite}
11\usepackage{hyperref}
12\usepackage{xcolor}
13\hypersetup{
16    citecolor={blue!50!black},
17    urlcolor={blue!80!black}
18}
19\usepackage{breqn}
20\usepackage{collect}
21\setlength{\parskip}{1ex}
22
23\title{Discretization of Poirier's Quantum Dynamical Equations of Motion}
24\author{Bill Page \small{(\href{mailto:bill.page@newsynthesis.org}{bill.page@newsynthesis.org})}}
25\date{\today}
26
27\begin{document}
28
29\maketitle
30\thispagestyle{empty}
31
32\begin{abstract}
33\textbf{From Poirier’s Bohmian Mechanics without Wavefunctions, to Hall’s Many Interacting Worlds in More Than One Dimension}
34
35In this note we show that Hall's "toy" model of quantum mechanics as many interacting interacting classical worlds 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.  Using the same approach we can construct similar models for the two- and three-dimensional case. We observe that the discretization of Poirier's expression for the quantum force is identical to the gradient of the discretization of the quantum potential in both the 1-D case and in the 2 and 3-D cases when using a symmetric difference operator.
36\end{abstract}
37
38\setcounter{section}{-1}
39\section{Preliminary}
40\begin{sagesilent}
41import numpy
42import mpmath
43\end{sagesilent}
44In this document we have used the computer algebra system
45Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of
46SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}
47and
48mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}
49for both symbolic and numeric calculations. The source code for this document can be found
50online\footnote{\url{https://cloud.sagemath.com/projects/b04b5777-e269-4c8f-a4b8-b21dbe1c93c6/files/discretization.tex}}. The inset bolded text below and many of the formulas in this paper are input and output produced directly by the Sage system.
51\begin{sageblock}
52from numpy import array,concatenate,isnan
53from mpmath import erfinv
54\end{sageblock}
55It is convenient to define a compact notation for functions and some special symbols:
56
57\definecollection{prelim}
58\begin{sagesilent}
59begincol=chr(92)+"begin { collect* } { prelim }"
60endcol=chr(92)+"end { collect* }"
61\end{sagesilent}
62
63ok { \sagestr{begincol} }
64\begin{sageblock}
65def fun(*y): # fun(x,y)(f)(x,y)=f
66    cache={}
67    def upd(f,x,y): # memoized
68        try: return cache[x]
69        except KeyError:
70            try: cache[x] = f.subs(dict([[y[i],x[i]] for i in range(len(y))]))
71            except AttributeError: # ask forgiveness
72                def sub(f): return f.subs(dict([[y[i],x[i]] for i in range(len(y))]))
73                cache[x] = map(sub,f)
74        return cache[x]
75    return lambda f: lambda *x: upd(f,x,y)
76def argscript(self, *args): return "%s_{%s}"%(self.name(),','.join(map(repr, args)))
77def functions(*names):
78    return map(lambda nam:function(nam, print_latex_func=argscript),names)
79def detJscript(self, *args):
80    return "\\begin{vmatrix}\\scriptstyle{J}\\end{vmatrix}_{%s}"% (
81        ','.join(map(repr, args)))
82detJ = function('detJ', print_latex_func=detJscript)
83def squared(x): return x.dot(x)  # vector (numpy array) x^2
84hbar = var('hbar',latex_name='\\hbar')
85mu = var('mu',latex_name='\\mu')
86\end{sageblock}
87\sagestr{endcol}
88
89
90The 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.
91\begin{sageblock}
92def fix(it):
93    return it.subs(dict([[function('x')(n+i,m+j),X[0](n+i,m+j)]
94                            for i in range(-1,2) for j in range(-1,2)] +
95                        [[function('y')(n+i,m+j),X[1](n+i,m+j)]
96                            for i in range(-1,2) for j in range(-1,2)]))
97\end{sageblock}
98
99\section{2-Dimensional System}
100In Schiff and Poirier \cite{schiff2012communication} eq. (18), the quantum dynamical equation of motion
101
102\begin{equation}
103\mu \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
104\end{equation}
105when specialized to 2-D and discretized by replacing derivatives with difference operators is the 2-D analog of Hall’s \cite{hall2014quantum} “toy” expression eq. (24).
106
107By 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.
108
109Let $E =[x_{n,m},y_{n,m}]$ represent the configuration space of one particle in 2-D space or two particles in 1-D space in the world with co-ordinates $n,m$. In Sage this can be represented as:
110\begin{sageblock}
111X = functions('x','y')
112d = len(X)
113C = var('n,m') # C=(n,m)
114E = tuple([x(*C) for x in X])
115\end{sageblock}
116\begin{equation} \sage{E} \end{equation}
117
118Forward and backward difference approximations to the derivatives are
119\begin{sageblock}
120def DM(x,n):return x-x.subs(n==n-1)
121def DP(x,n):return x.subs(n==n+1)-x
122\end{sageblock}
123Similarly the central difference approximation to this derivative is
124\begin{sageblock}
125def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2
126\end{sageblock}
127
128To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes:
129\begin{sageblock}
130J=fun(n,m)(matrix([[DM(E[j],C[i]) for i in range(d)] for j in range(d)]))
131\end{sageblock}
132\begin{equation} \sage{J(n,m)} \end{equation}
133It 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
134\begin{sageblock}
136\end{sageblock}
137\begin{equation} \sage{A(n,m)} \end{equation}
138The inverse is
139\begin{sageblock}
140K=fun(n,m)(A(n,m)/detJ(n,m))
141\end{sageblock}
142\begin{equation} \sage{K(n,m)} \end{equation}
143The determinant of the Jacobian is
144\begin{sageblock}
145Jdet=fun(n,m)(fix(J(n,m).determinant().expand()))
146dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-2,3) for j in range(-2,3)}
147\end{sageblock}
148% \sage{Jdet(n,m)} \end{equation}
149
150It is $\sage{fix(factor(J(n,m)^(-1)))==K(n,m).subs(detJ(n,m)==Jdet(n,m))}$ that $J(n,m)^{-1}=K(n,m) \text{ where }{\sage{detJ(n,m)}=\operatorname{Jdet}(n,m)}$.
151
152\subsection{Quantum Force}
153
154The Schiff and Poirier expression eq. (18) in 2-D has many more terms than in the 1-D case!
155
156\begin{sageblock}
157R=fun(n,m) (vector([ (hbar^2/4/mu)*sum( sum( sum( sum(
158 DP( K(n,m)[k,i] * K(n,m)[p,j] * DP( DM( K(n,m)[l,j],C[l] ),C[k] ),C[p] )
159     for p in range(d)) for k in range(d)) for j in range(d)) for l in range(d))
160         for i in range(d)])
161\end{sageblock}
162
163\subsection{Quantum Potential}
164
165\begin{sageblock}
166U=fun(n,m)(
167  sum( sum( sum(
168    -(hbar^2/4/mu) *
169    (K(n,m)[k,j] * DM( DP( K(n,m)[l,j],C[l] ),C[k] ) +
170      1/2 * DP(K(n,m)[l,j],D[l]) * DP(K(n,m)[k,j],C[k]))
171        for k in range(d)) for j in range(d)) for l in range(d)) )
172\end{sageblock}
173
174\begin{dmath}
175\sage{ eq20.subs({hbar:1,mu:1}) }
176\end{dmath}
177
178\begin{sageblock}
179U2=fun(n,m)(
180  sum( sum( sum( -(hbar^2/4/mu)*(1/2 * DP(K(n,m)[l,j],C[l]) * DP(K(n,m)[k,j],C[k]))
181    for k in range(d)) for j in range(d)) for l in range(d)) )
182\end{sageblock}
183
184\begin{dmath}
185\sage{ eq20n2 }
186\end{dmath}
187
188\begin{sageblock}
189U3=fun(n,m)(
190  sum( sum( sum( -(hbar^2/4/mu)*(K(n,m)[k,j] * DM( DP( K(n,m)[l,j],C[l] ),C[k] ))
191    for k in range(d)) for j in range(d)) for l in range(d)) )
192\end{sageblock}
193
194\begin{dmath}
195\sage{ eq20n3 }
196\end{dmath}
197
198It is $\sage{bool(eq20 == eq20n2+eq20n3)}$ that $eq20 = eq20n2+eq20n3$.
199
200In 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.
201
202\begin{sageblock}
203xyk = [X[i](n+j,m+k)
204          for j in range(-2,3) for k in range(-2,3) for i in range(d)]
205xyv = [var(X[i].name()+str(2+j)+str(2+k))
206          for j in range(-2,3) for k in range(-2,3) for i in range(d)]
207\end{sageblock}
208
209\begin{sageblock}
210xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
211ynm = X[1](n,m).subs(dict(zip(xyk,xyv)))
212\end{sageblock}
213\begin{dmath}
214\sage{[xnm,ynm]}
215\end{dmath}
216\bibliographystyle{bibstyles/utphys}
217\bibliography{bibliography}
218\end{document}
219
220We can compute the discretized quantum force as the gradient of the discretized quantum potential symbolically. We would like to compare this with the discretization of Poirier's expression for the force.
221
222Let $U(n,m)$ be the potential at $n,m$ with $U2(n,m)$ and $U3(n,m)$ being just the 2nd and 3rd order terms respectively.
223\begin{sageblock}
224def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
225def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
226def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
227\end{sageblock}
229\begin{sageblock}
230dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-2,3) for j in range(-1,3)}
231rrx = -sum(U(n+i,m+j) for i in range(-2,3) for j in range(-2,3)
232          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
233rry = -sum([U(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
234          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
235\end{sageblock}
236Or computing the gradients of the 2nd and 3rd order terms separately:
237\begin{sageblock}
238rr2x = -sum([U2(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
239           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
240rr2y = -sum([U2(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
241           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
242rr3x = -sum([U3(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
243           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
244rr3y = -sum([U3(n+i,m+j) for i in range(-2,3) for j in range(-2,3)]
245           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
246\end{sageblock}
247
248It is $\sage{bool(rr2x.expand()+rr3x.expand() == rrx.expand())}$ that $rrx = rr2x+rr3x$.
249
250It is $\sage{bool(rr2y.expand()+rr3y.expand() == rry.expand())}$ that $rrx = rr2y+rr3y$.
251
252Surprisingly, just like in the 1-D case, 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 is differs from the 2nd term by only a sign.
253
254It is $\sage{bool(rr3x.expand() == -2*rr2x.expand())}$ that $rr3x = -2 rr2x$.
255
256It is $\sage{bool(rr3y.expand() == -2*rr2y.expand())}$ that $rr3y = -2 rr2y$.
257
258But unlike the 1-D case, the result for the full expression is not the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression (eq18n) for the quantum force!
259
260\begin{sageblock}
261test=dict(zip(xyv,[random() for i in xyv]))
262\end{sageblock}
263
264It is $\sage{bool(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rrx.subs(test))}$ that $eq18[0] = rrx$.
265
266It is $\sage{bool(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rry.subs(test))}$ that $eq18[1] = rry$.
267
268Repeating all of the above using a symmetric difference approximation to the derivatives we find however that the two expressions for the quantum force turn out to be the same.
269
270\begin{sageblock}
271@CachedFunction
272def J(n,m):
273    return matrix([[D(E[j],i).subs({var('n'):n,var('m'):m})
274             for i in range(d)]
275               for j in range(d)])
276@CachedFunction
277def A(n,m): return fix(J(n,m).adjoint()).subs({var('n'):n,var('m'):m})
278@CachedFunction
279def K(n,m): return (A(n,m)/detJ(n,m)).subs({var('n'):n,var('m'):m})
280@CachedFunction
281def Jdet(n,m): return J(n,m).determinant().expand()
282\end{sageblock}
283\begin{equation} \sage{J(n,m)} \end{equation}
284\begin{equation} \sage{A(n,m)} \end{equation}
285\begin{equation} \sage{K(n,m)} \end{equation}
286\begin{sageblock}
287eq18 = [ (hbar^2/4/mu)*sum( sum( sum( sum(
288 D( K(n,m)[k,i] * K(n,m)[p,j] * D( D( K(n,m)[l,j],l ),k ),p )
289   for p in range(d))
290     for k in range(d))
291       for j in range(d))
292         for l in range(d)) for i in range(d)]
293eq20 = sum( sum( sum(
294 -(hbar^2/4/mu) *
295   (K(n,m)[k,j] * D( D( K(n,m)[l,j],l ),k ) +
296     1/2 * D(K(n,m)[l,j],l) * D(K(n,m)[k,j],k))
297       for k in range(d)) for j in range(d)) for l in range(d))
298eq20n2 = sum( sum( sum( -(hbar^2/4/mu)*(1/2 * DP(K(n,m)[l,j],l) * DP(K(n,m)[k,j],k))
299   for k in range(d)) for j in range(d)) for l in range(d))
300eq20n3 = sum( sum( sum( -(hbar^2/4/mu)*(K(n,m)[k,j] * DM( DP( K(n,m)[l,j],l ),k ))
301   for k in range(d)) for j in range(d)) for l in range(d))
302\end{sageblock}
303
304\begin{dmath} \sage{eq20} \end{dmath}
305
306\begin{sageblock}
307xyk = [X[i](n+j,m+k)
308          for j in range(-4,5) for k in range(-4,5) for i in range(d)]
309xyv = [var(X[i].name()+str(5+j)+str(5+k))
310          for j in range(-4,5) for k in range(-4,5) for i in range(d)]
311xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
312ynm = X[1](n,m).subs(dict(zip(xyk,xyv)))
313\end{sageblock}
314
315\begin{dmath}
316\sage{[xnm,ynm]}
317\end{dmath}
318
319\begin{sageblock}
320def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
321def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
322def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
323dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-6,7) for j in range(-6,7)}
324rrx = -sum(U(n+i,m+j) for i in range(-3,4) for j in range(-3,4)
325          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
326rry = -sum([U(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
327          ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
328rr2x = -sum([U2(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
329           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
330rr2y = -sum([U2(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
331           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
332rr3x = -sum([U3(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
333           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(xnm)
334rr3y = -sum([U3(n+i,m+j) for i in range(-3,4) for j in range(-3,4)]
335           ).subs(dets).subs(dict(zip(xyk,xyv))).diff(ynm)
336
337\end{sageblock}
338
339\begin{sageblock}
340vals=[RDF(random()) for i in xyv]
341test=dict(zip(xyv,vals))
342constants = {hbar:0.2,mu:1.0}
343\end{sageblock}
344
345
346It is $\sage{squared(map(lambda x:x.subs(test).subs(constants),array([rr2x,rr2y]) - array([rrx,rry])))<1e-15}$ that $[rrx,rry] = [rr2x,rr2y]+[rr3x,rr3y]$.
347
348The force term from the 3rd order term is still exactly -2 times the force term obtained from the 2nd order term.
349
350It is $\sage{squared(map(lambda x:x.subs(test).subs(constants),array([rr3x,rr3y]) - -2*array([rr2x,rr2y])))<1e-15}$ that $[rr3x,rr3y] = -2 [rr2x,rr2y]$.
351
352But now the full expression is the same as the expression we obtained directly by substituting difference operators for derivatives in Poirier's expression (eq18n) for the quantum force. It is very time consuming to show this symbolically but numerical serve as a practical demonstration.
353
354\begin{sageblock}
355ff0=fast_callable(rrx.subs(constants), vars=xyv, domain=RDF)
356ff1=fast_callable(rry.subs(constants), vars=xyv, domain=RDF)
357def ff(p): return array([ff0(*p),ff1(*p)])
358gg0=fast_callable(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
359gg1=fast_callable(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
360def gg(p): return array([gg0(*p),gg1(*p)])
361\end{sageblock}
362
363It is $\sage{ abs(ff0(*vals)-gg0(*vals)) < 1.0e-12 and abs(ff1(*vals)-gg1(*vals)) < 1.0e-12 }$ that $eq18 = [rrx,rry]$.
364
365\section{Simulation}
366
367Using the expressions derived above we can generate efficient numerical functions for use in the solution of the equations of motion.
368
369Classical Potential
370\begin{sageblock}
371    P0(x)=-9/(2*cosh(2*x)^2)
372\end{sageblock}
373\sageplot[width=.5\textwidth]{
374	plot(P0,(-3,3),axes_labels=['distance','energy'],
375	axes_labels_size=1,legend_label='potential')
376}
377
378Classical Force
379\begin{sageblock}
380	F0(x)= -(18*sinh(2*x))/(cosh(2*x)^3)
381\end{sageblock}
382\sageplot[width=.5\textwidth]{
383    plot(F0,(-3,3),axes_labels=['distance','force'],axes_labels_size=1)
384}
385
386Quantum Potential
387
388$\sage{U(n,m).subs(dets).subs(dict(zip(xyk,xyv))).subs(constants).variables()}$
389
390%\begin{sageblock}
391%ru0=fast_callable(U(n,m).subs(dets).subs(dict(zip(xyk,xyv))).subs(constants),
392%vars=xyv, domain=RDF)
393%def ru(x): return ru0(*x)
394%\end{sageblock}
395
396\bibliographystyle{bibstyles/utphys}
397\bibliography{bibliography}
398\end{document}
399%sagemathcloud={"latex_command":"pdflatex -synctex=1 -interact=nonstopmode 'discretization.tex'"}
400