CoCalc Shared Filesdiscretization.texOpen in CoCalc with one click!
Author: Bill Page
Views : 26
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{
14
colorlinks,
15
linkcolor={red!50!black},
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
35
In 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}
41
import numpy
42
import mpmath
43
\end{sagesilent}
44
In this document we have used the computer algebra system
45
Sage\footnote{\sagestr{version()}\cite{SageMath2016} on SageMathCloud\cite{SageMathCloud2016}} with components of
46
SciPy\footnote{numpy Version \sagestr{numpy.version.version} \cite{SciPy2016}}
47
and
48
mpmath\footnote{mpmath Version \sagestr{mpmath.__version__} \cite{mpmath2016}}
49
for both symbolic and numeric calculations. The source code for this document can be found
50
online\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}
52
from numpy import array,concatenate,isnan
53
from mpmath import erfinv
54
\end{sageblock}
55
It is convenient to define a compact notation for functions and some special symbols:
56
57
\definecollection{prelim}
58
\begin{sagesilent}
59
begincol=chr(92)+"begin { collect* } { prelim }"
60
endcol=chr(92)+"end { collect* }"
61
\end{sagesilent}
62
63
ok { \sagestr{begincol} }
64
\begin{sageblock}
65
def 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)
76
def argscript(self, *args): return "%s_{%s}"%(self.name(),','.join(map(repr, args)))
77
def functions(*names):
78
return map(lambda nam:function(nam, print_latex_func=argscript),names)
79
def detJscript(self, *args):
80
return "\\begin{vmatrix}\\scriptstyle{J}\\end{vmatrix}_{%s}"% (
81
','.join(map(repr, args)))
82
detJ = function('detJ', print_latex_func=detJscript)
83
def squared(x): return x.dot(x) # vector (numpy array) x^2
84
hbar = var('hbar',latex_name='\\hbar')
85
mu = var('mu',latex_name='\\mu')
86
\end{sageblock}
87
\sagestr{endcol}
88
89
90
The 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}
92
def 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}
100
In 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}
105
when 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
107
By 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
109
Let $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}
111
X = functions('x','y')
112
d = len(X)
113
C = var('n,m') # C=(n,m)
114
E = tuple([x(*C) for x in X])
115
\end{sageblock}
116
\begin{equation} \sage{E} \end{equation}
117
118
Forward and backward difference approximations to the derivatives are
119
\begin{sageblock}
120
def DM(x,n):return x-x.subs(n==n-1)
121
def DP(x,n):return x.subs(n==n+1)-x
122
\end{sageblock}
123
Similarly the central difference approximation to this derivative is
124
\begin{sageblock}
125
def D(x,n):return x.subs(n==n+1)/2-x.subs(n==n-1)/2
126
\end{sageblock}
127
128
To derive Hall's equations for a finite number of interacting worlds we replace derivatives with difference operators so the Jacobian becomes:
129
\begin{sageblock}
130
J=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}
133
It 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}
135
A=fun(n,m)(fix(J(n,m).adjoint()))
136
\end{sageblock}
137
\begin{equation} \sage{A(n,m)} \end{equation}
138
The inverse is
139
\begin{sageblock}
140
K=fun(n,m)(A(n,m)/detJ(n,m))
141
\end{sageblock}
142
\begin{equation} \sage{K(n,m)} \end{equation}
143
The determinant of the Jacobian is
144
\begin{sageblock}
145
Jdet=fun(n,m)(fix(J(n,m).determinant().expand()))
146
dets = {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
%\begin{equation} \sage{Jdet(n,m)} \end{equation}
149
150
It 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
154
The Schiff and Poirier expression eq. (18) in 2-D has many more terms than in the 1-D case!
155
156
\begin{sageblock}
157
R=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}
166
U=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}
179
U2=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}
189
U3=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
198
It is $\sage{bool(eq20 == eq20n2+eq20n3)}$ that $eq20 = eq20n2+eq20n3$.
199
200
In 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}
203
xyk = [X[i](n+j,m+k)
204
for j in range(-2,3) for k in range(-2,3) for i in range(d)]
205
xyv = [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}
210
xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
211
ynm = 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
220
We 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
222
Let $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}
224
def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
225
def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
226
def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
227
\end{sageblock}
228
Gradients:
229
\begin{sageblock}
230
dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-2,3) for j in range(-1,3)}
231
rrx = -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)
233
rry = -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}
236
Or computing the gradients of the 2nd and 3rd order terms separately:
237
\begin{sageblock}
238
rr2x = -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)
240
rr2y = -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)
242
rr3x = -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)
244
rr3y = -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
248
It is $\sage{bool(rr2x.expand()+rr3x.expand() == rrx.expand())}$ that $rrx = rr2x+rr3x$.
249
250
It is $\sage{bool(rr2y.expand()+rr3y.expand() == rry.expand())}$ that $rrx = rr2y+rr3y$.
251
252
Surprisingly, 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
254
It is $\sage{bool(rr3x.expand() == -2*rr2x.expand())}$ that $rr3x = -2 rr2x$.
255
256
It is $ \sage{bool(rr3y.expand() == -2*rr2y.expand())}$ that $rr3y = -2 rr2y$.
257
258
But 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}
261
test=dict(zip(xyv,[random() for i in xyv]))
262
\end{sageblock}
263
264
It is $\sage{bool(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rrx.subs(test))}$ that $eq18[0] = rrx$.
265
266
It is $\sage{bool(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(test)==rry.subs(test))}$ that $eq18[1] = rry$.
267
268
Repeating 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
272
def 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
277
def A(n,m): return fix(J(n,m).adjoint()).subs({var('n'):n,var('m'):m})
278
@CachedFunction
279
def K(n,m): return (A(n,m)/detJ(n,m)).subs({var('n'):n,var('m'):m})
280
@CachedFunction
281
def 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}
287
eq18 = [ (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)]
293
eq20 = 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))
298
eq20n2 = 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))
300
eq20n3 = 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}
307
xyk = [X[i](n+j,m+k)
308
for j in range(-4,5) for k in range(-4,5) for i in range(d)]
309
xyv = [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)]
311
xnm = X[0](n,m).subs(dict(zip(xyk,xyv)))
312
ynm = 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}
320
def U(n,m): return eq20.subs({var('n'):n,var('m'):m})
321
def U2(n,m): return eq20n2.subs({var('n'):n,var('m'):m})
322
def U3(n,m): return eq20n3.subs({var('n'):n,var('m'):m})
323
dets = {detJ(n+i,m+j):Jdet(n+i,m+j) for i in range(-6,7) for j in range(-6,7)}
324
rrx = -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)
326
rry = -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)
328
rr2x = -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)
330
rr2y = -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)
332
rr3x = -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)
334
rr3y = -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}
340
vals=[RDF(random()) for i in xyv]
341
test=dict(zip(xyv,vals))
342
constants = {hbar:0.2,mu:1.0}
343
\end{sageblock}
344
345
346
It 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
348
The force term from the 3rd order term is still exactly -2 times the force term obtained from the 2nd order term.
349
350
It 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
352
But 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}
355
ff0=fast_callable(rrx.subs(constants), vars=xyv, domain=RDF)
356
ff1=fast_callable(rry.subs(constants), vars=xyv, domain=RDF)
357
def ff(p): return array([ff0(*p),ff1(*p)])
358
gg0=fast_callable(eq18[0].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
359
gg1=fast_callable(eq18[1].subs(dets).subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
360
def gg(p): return array([gg0(*p),gg1(*p)])
361
\end{sageblock}
362
363
It 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
367
Using the expressions derived above we can generate efficient numerical functions for use in the solution of the equations of motion.
368
369
Classical 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
378
Classical 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
386
Quantum 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