Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

Each moving point in the following diagram represents the position of two particles in a single world (four worlds are shown). As the point approaches the x=y line (shown in red), the particles approach each other and the classical potential is at a maximum. In three worlds the particles overcome the potential barrier and continue. But one world (the last to approach) is reflected back. Although initially the two particles in each world approach each other with the same velocity, the particles in the last world to approach the barrier have lost energy to their counterparts in other worlds due to the quantum mechanical (interworld) interaction. When viewed in the projection we see two particles rebound instead of passing each other.

Views: 110

From Poirier's Bohmian Mechanics without Wavefunctions, to Hall's Many Interacting Worlds -- in More Than One Dimension

Ref:

  1. Quantum Mechanics Without Wavefunctions Jeremy Schiff and Bill Poirier J. Chem. Phys. 136, 031102 (2012)

  2. Quantum Phenomena Modeled by Interactions between Many Classical Worlds Michael J. W. Hall, Dirk-André Deckert and Howard M. Wiseman, PHYSICAL REVIEW X 4, 041013 (23 October 2014)

  3. Verlet integration (Wikipedia)

  4. Explicit, Time Reversible, Adaptive Step Size Control Ernst Hairer and Gustaf Söderlind SIAM Journal on Scientific Computing. 2005, vol. 26, no. 6, p. 1838-1851

#%typeset_mode True
from numpy import array,concatenate,isnan
from mpmath import erfinv
hbar=var('hbar',latex_name='\\hbar')
mu=var('mu',latex_name='\mu')

2-D

Poirier’s many-D generalization, eq. (18), when specialized to 2-D and discretized by replacing derivatives with difference operators is the 2-D analog of Hall’s “toy” expression eq. (24). Let x(n,m)x(n,m) and y(n,m)y(n,m) represent the position of a particle (one particle in 2-D space) or the position of two particles (in 1-D configuration space) at an unspecified time in the world indexed by nn and mm.

vars = ['x','y']; d = len(vars)
X = map(function,vars); x,y = X
n,m = var('n,m'); ind=[n,m]
# position of particle in world (n,m)
Enm=map(lambda x:x(*ind),X);Enm
[x(n, m), y(n, m)]

For Hall's finite number of interacting worlds we will use the following difference operators instead of derivatives.

def Dminus(x,i):return(x-x.subs(ind[i]==ind[i]-1))
def Dplus(x,i):return(x.subs(ind[i]==ind[i]+1)-x)

The Jacobian replacing derivatives with differences.

Jnm = matrix(map(lambda e:[Dminus(e,i) for i in range(d)],Enm));Jnm
[-x(n - 1, m) + x(n, m) -x(n, m - 1) + x(n, m)] [-y(n - 1, m) + y(n, m) -y(n, m - 1) + y(n, m)] 
# inverse Jacobian
Knm = matrix(map(lambda x: map(lambda y:y.factor(),x),Jnm^(-1)))

# determinant
Dnm=Jnm.determinant(); Dnm
-(x(n, m - 1) - x(n, m))*(y(n - 1, m) - y(n, m)) + (x(n - 1, m) - x(n, m))*(y(n, m - 1) - y(n, m)) 
# minors
Mnm=matrix(map(lambda x: map(lambda y:(y*Dnm).full_simplify(),x),Knm));Mnm
[-y(n, m - 1) + y(n, m) x(n, m - 1) - x(n, m)] [ y(n - 1, m) - y(n, m) -x(n - 1, m) + x(n, m)] 
# inverse is minors divided by determinant
def Jdet(n,m): return Dnm.subs({var('n'):n,var('m'):m})
map(lambda x:map(lambda y:y.full_simplify(),x),Mnm/Dnm-Knm)==[[0,0],[0,0]]
True

Quantum Force

The Schiff and Poirier expression eq. (18) in 2-D has many more terms than in the 1-D case!

# introduce symbolic determinant
det=function('det')
eq18n = [ (hbar^2/4/mu)*sum([ sum([ sum([ sum([
 Dplus( Mnm[k,i]/det(n,m)*Mnm[p,j]/det(n,m)*Dplus( Dminus( Mnm[l,j]/det(n,m),l ),k ),p )
   for p in range(d)]) for k in range(d)]) for j in range(d)]) for l in range(d)]) for i in range(d)]
# Too large to conveniently display
# eq18n

# the denominators of both components are the same
show(eq18n[0].denominator())
bool(eq18n[0].denominator()==eq18n[1].denominator())
4 μdet(n+2,m−1)det(n+2,m)det(n+1,m+1)det(n+1,m−1)det(n+1,m)3det(n−1,m+2)det(n−1,m+1)det(n−1,m)det(n,m+2)det(n,m+1)3det(n,m−1)det(n,m)3\displaystyle 4 \, {{\mu}} {\rm det}\left(n + 2, m - 1\right) {\rm det}\left(n + 2, m\right) {\rm det}\left(n + 1, m + 1\right) {\rm det}\left(n + 1, m - 1\right) {\rm det}\left(n + 1, m\right)^{3} {\rm det}\left(n - 1, m + 2\right) {\rm det}\left(n - 1, m + 1\right) {\rm det}\left(n - 1, m\right) {\rm det}\left(n, m + 2\right) {\rm det}\left(n, m + 1\right)^{3} {\rm det}\left(n, m - 1\right) {\rm det}\left(n, m\right)^{3}
True 
# expression components are rational functions
eq18numer0=eq18n[0].numerator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)})
eq18denom0=eq18n[0].denominator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)})
eq18numer1=eq18n[1].numerator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)})
eq18denom1=eq18n[1].denominator().subs({det(n+i,m+j):Jdet(n+i,m+j) for i in range(-1,3) for j in range(-1,3)})

# denominator explicitly
show(eq18denom0)
4 ((x(n+2,m−1)−x(n+1,m−1))(y(n+2,m−1)−y(n+2,m−2))−(x(n+2,m−1)−x(n+2,m−2))(y(n+2,m−1)−y(n+1,m−1)))((x(n+2,m)−x(n+1,m))(y(n+2,m−1)−y(n+2,m))−(x(n+2,m−1)−x(n+2,m))(y(n+2,m)−y(n+1,m)))((x(n+1,m+1)−x(n,m+1))(y(n+1,m+1)−y(n+1,m))−(x(n+1,m+1)−x(n+1,m))(y(n+1,m+1)−y(n,m+1)))((x(n+1,m−1)−x(n,m−1))(y(n+1,m−1)−y(n+1,m−2))−(x(n+1,m−1)−x(n+1,m−2))(y(n+1,m−1)−y(n,m−1)))((x(n+1,m)−x(n,m))(y(n+1,m−1)−y(n+1,m))−(x(n+1,m−1)−x(n+1,m))(y(n+1,m)−y(n,m)))3((x(n−1,m+2)−x(n−2,m+2))(y(n−1,m+2)−y(n−1,m+1))−(x(n−1,m+2)−x(n−1,m+1))(y(n−1,m+2)−y(n−2,m+2)))((x(n,m+2)−x(n,m+1))(y(n−1,m+2)−y(n,m+2))−(x(n−1,m+2)−x(n,m+2))(y(n,m+2)−y(n,m+1)))((x(n−1,m+1)−x(n−2,m+1))(y(n−1,m+1)−y(n−1,m))−(x(n−1,m+1)−x(n−1,m))(y(n−1,m+1)−y(n−2,m+1)))((x(n,m+1)−x(n,m))(y(n−1,m+1)−y(n,m+1))−(x(n−1,m+1)−x(n,m+1))(y(n,m+1)−y(n,m)))3((x(n−1,m)−x(n−2,m))(y(n−1,m−1)−y(n−1,m))−(x(n−1,m−1)−x(n−1,m))(y(n−1,m)−y(n−2,m)))((x(n,m−1)−x(n,m−2))(y(n−1,m−1)−y(n,m−1))−(x(n−1,m−1)−x(n,m−1))(y(n,m−1)−y(n,m−2)))((x(n,m−1)−x(n,m))(y(n−1,m)−y(n,m))−(x(n−1,m)−x(n,m))(y(n,m−1)−y(n,m)))3μ\displaystyle 4 \, {\left({\left(x\left(n + 2, m - 1\right) - x\left(n + 1, m - 1\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 2, m - 2\right)\right)} - {\left(x\left(n + 2, m - 1\right) - x\left(n + 2, m - 2\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 1, m - 1\right)\right)}\right)} {\left({\left(x\left(n + 2, m\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 2, m - 1\right) - y\left(n + 2, m\right)\right)} - {\left(x\left(n + 2, m - 1\right) - x\left(n + 2, m\right)\right)} {\left(y\left(n + 2, m\right) - y\left(n + 1, m\right)\right)}\right)} {\left({\left(x\left(n + 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n + 1, m + 1\right) - y\left(n + 1, m\right)\right)} - {\left(x\left(n + 1, m + 1\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 1, m + 1\right) - y\left(n, m + 1\right)\right)}\right)} {\left({\left(x\left(n + 1, m - 1\right) - x\left(n, m - 1\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m - 2\right)\right)} - {\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m - 2\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n, m - 1\right)\right)}\right)} {\left({\left(x\left(n + 1, m\right) - x\left(n, m\right)\right)} {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} - {\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(y\left(n + 1, m\right) - y\left(n, m\right)\right)}\right)}^{3} {\left({\left(x\left(n - 1, m + 2\right) - x\left(n - 2, m + 2\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n - 1, m + 1\right)\right)} - {\left(x\left(n - 1, m + 2\right) - x\left(n - 1, m + 1\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n - 2, m + 2\right)\right)}\right)} {\left({\left(x\left(n, m + 2\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n - 1, m + 2\right) - y\left(n, m + 2\right)\right)} - {\left(x\left(n - 1, m + 2\right) - x\left(n, m + 2\right)\right)} {\left(y\left(n, m + 2\right) - y\left(n, m + 1\right)\right)}\right)} {\left({\left(x\left(n - 1, m + 1\right) - x\left(n - 2, m + 1\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n - 1, m\right)\right)} - {\left(x\left(n - 1, m + 1\right) - x\left(n - 1, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n - 2, m + 1\right)\right)}\right)} {\left({\left(x\left(n, m + 1\right) - x\left(n, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)} - {\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(y\left(n, m + 1\right) - y\left(n, m\right)\right)}\right)}^{3} {\left({\left(x\left(n - 1, m\right) - x\left(n - 2, m\right)\right)} {\left(y\left(n - 1, m - 1\right) - y\left(n - 1, m\right)\right)} - {\left(x\left(n - 1, m - 1\right) - x\left(n - 1, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n - 2, m\right)\right)}\right)} {\left({\left(x\left(n, m - 1\right) - x\left(n, m - 2\right)\right)} {\left(y\left(n - 1, m - 1\right) - y\left(n, m - 1\right)\right)} - {\left(x\left(n - 1, m - 1\right) - x\left(n, m - 1\right)\right)} {\left(y\left(n, m - 1\right) - y\left(n, m - 2\right)\right)}\right)} {\left({\left(x\left(n, m - 1\right) - x\left(n, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m\right)\right)} - {\left(x\left(n - 1, m\right) - x\left(n, m\right)\right)} {\left(y\left(n, m - 1\right) - y\left(n, m\right)\right)}\right)}^{3} {{\mu}}

Create a numeric function that depends on just the nearest neighbors. This requires a function of 2x5x5 = 50 scalar quantities.

xyk=[X[i](n+j,m+k) for j in [-2,-1,0,1,2] for k in [-2,-1,0,1,2] for i in range(d)]
xyv=[var(vars[i]+str(2+j)+str(2+k)) for j in [-2,-1,0,1,2] for k in [-2,-1,0,1,2] for i in range(d)]
show(zip(xyk,xyv))
[(x(n−2,m−2)\displaystyle x\left(n - 2, m - 2\right), x00\displaystyle x_{00}), (y(n−2,m−2)\displaystyle y\left(n - 2, m - 2\right), y00\displaystyle y_{00}), (x(n−2,m−1)\displaystyle x\left(n - 2, m - 1\right), x01\displaystyle x_{01}), (y(n−2,m−1)\displaystyle y\left(n - 2, m - 1\right), y01\displaystyle y_{01}), (x(n−2,m)\displaystyle x\left(n - 2, m\right), x02\displaystyle x_{02}), (y(n−2,m)\displaystyle y\left(n - 2, m\right), y02\displaystyle y_{02}), (x(n−2,m+1)\displaystyle x\left(n - 2, m + 1\right), x03\displaystyle x_{03}), (y(n−2,m+1)\displaystyle y\left(n - 2, m + 1\right), y03\displaystyle y_{03}), (x(n−2,m+2)\displaystyle x\left(n - 2, m + 2\right), x04\displaystyle x_{04}), (y(n−2,m+2)\displaystyle y\left(n - 2, m + 2\right), y04\displaystyle y_{04}), (x(n−1,m−2)\displaystyle x\left(n - 1, m - 2\right), x10\displaystyle x_{10}), (y(n−1,m−2)\displaystyle y\left(n - 1, m - 2\right), y10\displaystyle y_{10}), (x(n−1,m−1)\displaystyle x\left(n - 1, m - 1\right), x11\displaystyle x_{11}), (y(n−1,m−1)\displaystyle y\left(n - 1, m - 1\right), y11\displaystyle y_{11}), (x(n−1,m)\displaystyle x\left(n - 1, m\right), x12\displaystyle x_{12}), (y(n−1,m)\displaystyle y\left(n - 1, m\right), y12\displaystyle y_{12}), (x(n−1,m+1)\displaystyle x\left(n - 1, m + 1\right), x13\displaystyle x_{13}), (y(n−1,m+1)\displaystyle y\left(n - 1, m + 1\right), y13\displaystyle y_{13}), (x(n−1,m+2)\displaystyle x\left(n - 1, m + 2\right), x14\displaystyle x_{14}), (y(n−1,m+2)\displaystyle y\left(n - 1, m + 2\right), y14\displaystyle y_{14}), (x(n,m−2)\displaystyle x\left(n, m - 2\right), x20\displaystyle x_{20}), (y(n,m−2)\displaystyle y\left(n, m - 2\right), y20\displaystyle y_{20}), (x(n,m−1)\displaystyle x\left(n, m - 1\right), x21\displaystyle x_{21}), (y(n,m−1)\displaystyle y\left(n, m - 1\right), y21\displaystyle y_{21}), (x(n,m)\displaystyle x\left(n, m\right), x22\displaystyle x_{22}), (y(n,m)\displaystyle y\left(n, m\right), y22\displaystyle y_{22}), (x(n,m+1)\displaystyle x\left(n, m + 1\right), x23\displaystyle x_{23}), (y(n,m+1)\displaystyle y\left(n, m + 1\right), y23\displaystyle y_{23}), (x(n,m+2)\displaystyle x\left(n, m + 2\right), x24\displaystyle x_{24}), (y(n,m+2)\displaystyle y\left(n, m + 2\right), y24\displaystyle y_{24}), (x(n+1,m−2)\displaystyle x\left(n + 1, m - 2\right), x30\displaystyle x_{30}), (y(n+1,m−2)\displaystyle y\left(n + 1, m - 2\right), y30\displaystyle y_{30}), (x(n+1,m−1)\displaystyle x\left(n + 1, m - 1\right), x31\displaystyle x_{31}), (y(n+1,m−1)\displaystyle y\left(n + 1, m - 1\right), y31\displaystyle y_{31}), (x(n+1,m)\displaystyle x\left(n + 1, m\right), x32\displaystyle x_{32}), (y(n+1,m)\displaystyle y\left(n + 1, m\right), y32\displaystyle y_{32}), (x(n+1,m+1)\displaystyle x\left(n + 1, m + 1\right), x33\displaystyle x_{33}), (y(n+1,m+1)\displaystyle y\left(n + 1, m + 1\right), y33\displaystyle y_{33}), (x(n+1,m+2)\displaystyle x\left(n + 1, m + 2\right), x34\displaystyle x_{34}), (y(n+1,m+2)\displaystyle y\left(n + 1, m + 2\right), y34\displaystyle y_{34}), (x(n+2,m−2)\displaystyle x\left(n + 2, m - 2\right), x40\displaystyle x_{40}), (y(n+2,m−2)\displaystyle y\left(n + 2, m - 2\right), y40\displaystyle y_{40}), (x(n+2,m−1)\displaystyle x\left(n + 2, m - 1\right), x41\displaystyle x_{41}), (y(n+2,m−1)\displaystyle y\left(n + 2, m - 1\right), y41\displaystyle y_{41}), (x(n+2,m)\displaystyle x\left(n + 2, m\right), x42\displaystyle x_{42}), (y(n+2,m)\displaystyle y\left(n + 2, m\right), y42\displaystyle y_{42}), (x(n+2,m+1)\displaystyle x\left(n + 2, m + 1\right), x43\displaystyle x_{43}), (y(n+2,m+1)\displaystyle y\left(n + 2, m + 1\right), y43\displaystyle y_{43}), (x(n+2,m+2)\displaystyle x\left(n + 2, m + 2\right), x44\displaystyle x_{44}), (y(n+2,m+2)\displaystyle y\left(n + 2, m + 2\right), y44\displaystyle y_{44})]
# suitable values
constants = {hbar:0.2,mu:1.0}
# compile numeric functions
ff0=fast_callable(eq18numer0.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
ff1=fast_callable(eq18numer1.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
ffd=fast_callable(eq18denom0.subs(dict(zip(xyk,xyv))).subs(constants), vars=xyv, domain=RDF)
def ff(p): den=ffd(*p); return array([ff0(*p)/den,ff1(*p)/den])

Quantum Potential###

Why does Poirier eq. (20) include a second term? If we include it then total energy T+V+QT+V+Q is not conserved!

eq20nm = sum([ sum([ sum([
 (hbar^2/4/mu)*(Knm[k,j]*Dplus( Dplus( Knm[l,j],l ),k )+1/2*Dplus(Knm[l,j],l)*Dplus(Knm[k,j],k))
   for k in range(d)]) for j in range(d)]) for l in range(d)])
#too complex to display conveniently
#eq20nm

Expand terms using symbolic determinants.

eq20nm1 = sum([ sum([ sum([
 (hbar^2/4/mu)*(Mnm[k,j]/det(n,m)*Dplus( Dplus( Mnm[l,j]/det(n,m),l ),k ))
   for k in range(d)]) for j in range(d)]) for l in range(d)])
show(eq20nm1.subs({hbar:1,mu:1}).expand().combine().collect(det(n,m)).collect(det(n+2,m)).collect(det(n+1,m+1)).collect(det(n+1,m)).collect(det(n,m+1)) )
x(n−1,m)2−2 x(n−1,m)x(n,m−1)+x(n,m−1)2+y(n−1,m)2−2 y(n−1,m)y(n,m−1)+y(n,m−1)24 det(n,m)2+x(n+2,m−1)x(n,m−1)−x(n+2,m)x(n,m−1)−x(n+2,m−1)x(n,m)+x(n+2,m)x(n,m)+y(n+2,m−1)y(n,m−1)−y(n+2,m)y(n,m−1)−y(n+2,m−1)y(n,m)+y(n+2,m)y(n,m)4 det(n+2,m)det(n,m)+x(n+1,m+1)x(n−1,m)−x(n+1,m)x(n−1,m)+x(n+1,m+1)x(n,m−1)−x(n,m+1)x(n,m−1)−2 x(n+1,m+1)x(n,m)+x(n+1,m)x(n,m)+x(n,m+1)x(n,m)+y(n+1,m+1)y(n−1,m)−y(n+1,m)y(n−1,m)+y(n+1,m+1)y(n,m−1)−y(n,m+1)y(n,m−1)−2 y(n+1,m+1)y(n,m)+y(n+1,m)y(n,m)+y(n,m+1)y(n,m)4 det(n+1,m+1)det(n,m)+x(n+1,m−1)x(n−1,m)−x(n+1,m)x(n−1,m)−2 x(n+1,m−1)x(n,m−1)+x(n+1,m)x(n,m−1)+x(n+1,m−1)x(n,m)+x(n,m−1)x(n,m)−x(n,m)2+y(n+1,m−1)y(n−1,m)−y(n+1,m)y(n−1,m)−2 y(n+1,m−1)y(n,m−1)+y(n+1,m)y(n,m−1)+y(n+1,m−1)y(n,m)+y(n,m−1)y(n,m)−y(n,m)24 det(n+1,m)det(n,m)+x(n−1,m+2)x(n−1,m)−x(n−1,m)x(n,m+2)−x(n−1,m+2)x(n,m)+x(n,m+2)x(n,m)+y(n−1,m+2)y(n−1,m)−y(n−1,m)y(n,m+2)−y(n−1,m+2)y(n,m)+y(n,m+2)y(n,m)4 det(n,m+2)det(n,m)−2 x(n−1,m+1)x(n−1,m)−x(n−1,m)x(n,m+1)−x(n−1,m+1)x(n,m−1)+x(n,m+1)x(n,m−1)−x(n−1,m+1)x(n,m)−x(n−1,m)x(n,m)+x(n,m)2+2 y(n−1,m+1)y(n−1,m)−y(n−1,m)y(n,m+1)−y(n−1,m+1)y(n,m−1)+y(n,m+1)y(n,m−1)−y(n−1,m+1)y(n,m)−y(n−1,m)y(n,m)+y(n,m)24 det(n,m+1)det(n,m)\displaystyle \frac{x\left(n - 1, m\right)^{2} - 2 \, x\left(n - 1, m\right) x\left(n, m - 1\right) + x\left(n, m - 1\right)^{2} + y\left(n - 1, m\right)^{2} - 2 \, y\left(n - 1, m\right) y\left(n, m - 1\right) + y\left(n, m - 1\right)^{2}}{4 \, {\rm det}\left(n, m\right)^{2}} + \frac{x\left(n + 2, m - 1\right) x\left(n, m - 1\right) - x\left(n + 2, m\right) x\left(n, m - 1\right) - x\left(n + 2, m - 1\right) x\left(n, m\right) + x\left(n + 2, m\right) x\left(n, m\right) + y\left(n + 2, m - 1\right) y\left(n, m - 1\right) - y\left(n + 2, m\right) y\left(n, m - 1\right) - y\left(n + 2, m - 1\right) y\left(n, m\right) + y\left(n + 2, m\right) y\left(n, m\right)}{4 \, {\rm det}\left(n + 2, m\right) {\rm det}\left(n, m\right)} + \frac{x\left(n + 1, m + 1\right) x\left(n - 1, m\right) - x\left(n + 1, m\right) x\left(n - 1, m\right) + x\left(n + 1, m + 1\right) x\left(n, m - 1\right) - x\left(n, m + 1\right) x\left(n, m - 1\right) - 2 \, x\left(n + 1, m + 1\right) x\left(n, m\right) + x\left(n + 1, m\right) x\left(n, m\right) + x\left(n, m + 1\right) x\left(n, m\right) + y\left(n + 1, m + 1\right) y\left(n - 1, m\right) - y\left(n + 1, m\right) y\left(n - 1, m\right) + y\left(n + 1, m + 1\right) y\left(n, m - 1\right) - y\left(n, m + 1\right) y\left(n, m - 1\right) - 2 \, y\left(n + 1, m + 1\right) y\left(n, m\right) + y\left(n + 1, m\right) y\left(n, m\right) + y\left(n, m + 1\right) y\left(n, m\right)}{4 \, {\rm det}\left(n + 1, m + 1\right) {\rm det}\left(n, m\right)} + \frac{x\left(n + 1, m - 1\right) x\left(n - 1, m\right) - x\left(n + 1, m\right) x\left(n - 1, m\right) - 2 \, x\left(n + 1, m - 1\right) x\left(n, m - 1\right) + x\left(n + 1, m\right) x\left(n, m - 1\right) + x\left(n + 1, m - 1\right) x\left(n, m\right) + x\left(n, m - 1\right) x\left(n, m\right) - x\left(n, m\right)^{2} + y\left(n + 1, m - 1\right) y\left(n - 1, m\right) - y\left(n + 1, m\right) y\left(n - 1, m\right) - 2 \, y\left(n + 1, m - 1\right) y\left(n, m - 1\right) + y\left(n + 1, m\right) y\left(n, m - 1\right) + y\left(n + 1, m - 1\right) y\left(n, m\right) + y\left(n, m - 1\right) y\left(n, m\right) - y\left(n, m\right)^{2}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m\right)} + \frac{x\left(n - 1, m + 2\right) x\left(n - 1, m\right) - x\left(n - 1, m\right) x\left(n, m + 2\right) - x\left(n - 1, m + 2\right) x\left(n, m\right) + x\left(n, m + 2\right) x\left(n, m\right) + y\left(n - 1, m + 2\right) y\left(n - 1, m\right) - y\left(n - 1, m\right) y\left(n, m + 2\right) - y\left(n - 1, m + 2\right) y\left(n, m\right) + y\left(n, m + 2\right) y\left(n, m\right)}{4 \, {\rm det}\left(n, m + 2\right) {\rm det}\left(n, m\right)} - \frac{2 \, x\left(n - 1, m + 1\right) x\left(n - 1, m\right) - x\left(n - 1, m\right) x\left(n, m + 1\right) - x\left(n - 1, m + 1\right) x\left(n, m - 1\right) + x\left(n, m + 1\right) x\left(n, m - 1\right) - x\left(n - 1, m + 1\right) x\left(n, m\right) - x\left(n - 1, m\right) x\left(n, m\right) + x\left(n, m\right)^{2} + 2 \, y\left(n - 1, m + 1\right) y\left(n - 1, m\right) - y\left(n - 1, m\right) y\left(n, m + 1\right) - y\left(n - 1, m + 1\right) y\left(n, m - 1\right) + y\left(n, m + 1\right) y\left(n, m - 1\right) - y\left(n - 1, m + 1\right) y\left(n, m\right) - y\left(n - 1, m\right) y\left(n, m\right) + y\left(n, m\right)^{2}}{4 \, {\rm det}\left(n, m + 1\right) {\rm det}\left(n, m\right)}
eq20nm2 = sum([ sum([ sum([
 (hbar^2/4/mu)*(1/2*Dplus(Mnm[l,j]/det(n,m),l)*Dplus(Mnm[k,j]/det(n,m),k))
   for k in range(d)]) for j in range(d)]) for l in range(d)])

show(eq20nm2.subs({hbar:1,mu:1}))
18 (x(n+1,m−1)−x(n+1,m)det(n+1,m)−x(n,m−1)−x(n,m)det(n,m))2+18 (y(n+1,m−1)−y(n+1,m)det(n+1,m)−y(n,m−1)−y(n,m)det(n,m))2−14 (x(n+1,m−1)−x(n+1,m)det(n+1,m)−x(n,m−1)−x(n,m)det(n,m))(x(n−1,m+1)−x(n,m+1)det(n,m+1)−x(n−1,m)−x(n,m)det(n,m))+18 (x(n−1,m+1)−x(n,m+1)det(n,m+1)−x(n−1,m)−x(n,m)det(n,m))2−14 (y(n+1,m−1)−y(n+1,m)det(n+1,m)−y(n,m−1)−y(n,m)det(n,m))(y(n−1,m+1)−y(n,m+1)det(n,m+1)−y(n−1,m)−y(n,m)det(n,m))+18 (y(n−1,m+1)−y(n,m+1)det(n,m+1)−y(n−1,m)−y(n,m)det(n,m))2\displaystyle \frac{1}{8} \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} + \frac{1}{8} \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} - \frac{1}{4} \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} + \frac{1}{8} \, {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2} - \frac{1}{4} \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{{\rm det}\left(n + 1, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)} + \frac{1}{8} \, {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{{\rm det}\left(n, m + 1\right)} - \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{{\rm det}\left(n, m\right)}\right)}^{2}

Manually reorganizing terms, we can put it into another form.

eq20nm4=((x(n+1,m-1)-x(n+1,m))^2+(y(n+1,m-1)-y(n+1,m))^2)/8/det(n+1,m)^2 + ((x(n-1,m+1)-x(n,m+1))^2+(y(n-1,m+1)-y(n,m+1))^2)/8/det(n,m+1)^2 + ((x(n-1,m)-x(n,m-1))^2+(y(n-1,m)-y(n,m-1))^2)/8/det(n,m)^2 - ((x(n-1,m+1)-x(n,m+1))*(x(n-1,m)-x(n,m-1)) + (y(n-1,m+1)-y(n,m+1))*(y(n-1,m)-y(n,m-1)))/4/det(n,m+1)/det(n,m) + ((x(n+1,m-1)-x(n+1,m))*(x(n-1,m)-x(n,m-1)) + (y(n+1,m-1)-y(n+1,m))*(y(n-1,m)-y(n,m-1)))/4/det(n+1,m)/det(n,m) - ((x(n+1,m-1)-x(n+1,m))*(x(n-1,m+1)-x(n,m+1)) + (y(n+1,m-1)-y(n+1,m))*(y(n-1,m+1)-y(n,m+1)))/4/det(n+1,m)/det(n,m+1)
show(eq20nm4)
(x(n+1,m−1)−x(n+1,m))2+(y(n+1,m−1)−y(n+1,m))28 det(n+1,m)2+(x(n−1,m+1)−x(n,m+1))2+(y(n−1,m+1)−y(n,m+1))28 det(n,m+1)2−(x(n+1,m−1)−x(n+1,m))(x(n−1,m+1)−x(n,m+1))+(y(n+1,m−1)−y(n+1,m))(y(n−1,m+1)−y(n,m+1))4 det(n+1,m)det(n,m+1)+(x(n−1,m)−x(n,m−1))2+(y(n−1,m)−y(n,m−1))28 det(n,m)2+(x(n+1,m−1)−x(n+1,m))(x(n−1,m)−x(n,m−1))+(y(n+1,m−1)−y(n+1,m))(y(n−1,m)−y(n,m−1))4 det(n+1,m)det(n,m)−(x(n−1,m+1)−x(n,m+1))(x(n−1,m)−x(n,m−1))+(y(n−1,m+1)−y(n,m+1))(y(n−1,m)−y(n,m−1))4 det(n,m+1)det(n,m)\displaystyle \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)}^{2} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)}^{2}}{8 \, {\rm det}\left(n + 1, m\right)^{2}} + \frac{{\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)}^{2} + {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)}^{2}}{8 \, {\rm det}\left(n, m + 1\right)^{2}} - \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m + 1\right)} + \frac{{\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)}^{2} + {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}^{2}}{8 \, {\rm det}\left(n, m\right)^{2}} + \frac{{\left(x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)\right)} {\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)} + {\left(y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}}{4 \, {\rm det}\left(n + 1, m\right) {\rm det}\left(n, m\right)} - \frac{{\left(x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)\right)} {\left(x\left(n - 1, m\right) - x\left(n, m - 1\right)\right)} + {\left(y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)\right)} {\left(y\left(n - 1, m\right) - y\left(n, m - 1\right)\right)}}{4 \, {\rm det}\left(n, m + 1\right) {\rm det}\left(n, m\right)}
bool((eq20nm2.subs({hbar:1,mu:1})-eq20nm4).full_simplify()==0)
True

But for the computation we don't really need the determinants. (Only use the second term.)

eq20nm0 = sum([ sum([ sum([
 (hbar^2/4/mu)*(1/2*Dplus(Knm[l,j],l)*Dplus(Knm[k,j],k))
   for k in range(d)]) for j in range(d)]) for l in range(d)])
show(eq20nm0.collect(hbar).collect(mu))
((x(n+1,m−1)−x(n+1,m)x(n+1,m)y(n+1,m−1)−x(n,m)y(n+1,m−1)−x(n+1,m−1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m−1)y(n,m)−x(n+1,m)y(n,m)−x(n,m−1)−x(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))2+2 (x(n+1,m−1)−x(n+1,m)x(n+1,m)y(n+1,m−1)−x(n,m)y(n+1,m−1)−x(n+1,m−1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m−1)y(n,m)−x(n+1,m)y(n,m)−x(n,m−1)−x(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))(x(n−1,m+1)−x(n,m+1)x(n,m+1)y(n−1,m+1)−x(n,m)y(n−1,m+1)−x(n−1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n−1,m+1)y(n,m)−x(n,m+1)y(n,m)+x(n−1,m)−x(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))+(x(n−1,m+1)−x(n,m+1)x(n,m+1)y(n−1,m+1)−x(n,m)y(n−1,m+1)−x(n−1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n−1,m+1)y(n,m)−x(n,m+1)y(n,m)+x(n−1,m)−x(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))2+(y(n+1,m−1)−y(n+1,m)x(n+1,m)y(n+1,m−1)−x(n,m)y(n+1,m−1)−x(n+1,m−1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m−1)y(n,m)−x(n+1,m)y(n,m)−y(n,m−1)−y(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))2+2 (y(n+1,m−1)−y(n+1,m)x(n+1,m)y(n+1,m−1)−x(n,m)y(n+1,m−1)−x(n+1,m−1)y(n+1,m)+x(n,m)y(n+1,m)+x(n+1,m−1)y(n,m)−x(n+1,m)y(n,m)−y(n,m−1)−y(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))(y(n−1,m+1)−y(n,m+1)x(n,m+1)y(n−1,m+1)−x(n,m)y(n−1,m+1)−x(n−1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n−1,m+1)y(n,m)−x(n,m+1)y(n,m)+y(n−1,m)−y(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))+(y(n−1,m+1)−y(n,m+1)x(n,m+1)y(n−1,m+1)−x(n,m)y(n−1,m+1)−x(n−1,m+1)y(n,m+1)+x(n,m)y(n,m+1)+x(n−1,m+1)y(n,m)−x(n,m+1)y(n,m)+y(n−1,m)−y(n,m)x(n,m−1)y(n−1,m)−x(n,m)y(n−1,m)−x(n−1,m)y(n,m−1)+x(n,m)y(n,m−1)+x(n−1,m)y(n,m)−x(n,m−1)y(n,m))2)ℏ28 μ\displaystyle \frac{{\left({\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + 2 \, {\left(\frac{x\left(n + 1, m - 1\right) - x\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{x\left(n, m - 1\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} + {\left(\frac{x\left(n - 1, m + 1\right) - x\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{x\left(n - 1, m\right) - x\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2} + 2 \, {\left(\frac{y\left(n + 1, m - 1\right) - y\left(n + 1, m\right)}{x\left(n + 1, m\right) y\left(n + 1, m - 1\right) - x\left(n, m\right) y\left(n + 1, m - 1\right) - x\left(n + 1, m - 1\right) y\left(n + 1, m\right) + x\left(n, m\right) y\left(n + 1, m\right) + x\left(n + 1, m - 1\right) y\left(n, m\right) - x\left(n + 1, m\right) y\left(n, m\right)} - \frac{y\left(n, m - 1\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)} + {\left(\frac{y\left(n - 1, m + 1\right) - y\left(n, m + 1\right)}{x\left(n, m + 1\right) y\left(n - 1, m + 1\right) - x\left(n, m\right) y\left(n - 1, m + 1\right) - x\left(n - 1, m + 1\right) y\left(n, m + 1\right) + x\left(n, m\right) y\left(n, m + 1\right) + x\left(n - 1, m + 1\right) y\left(n, m\right) - x\left(n, m + 1\right) y\left(n, m\right)} + \frac{y\left(n - 1, m\right) - y\left(n, m\right)}{x\left(n, m - 1\right) y\left(n - 1, m\right) - x\left(n, m\right) y\left(n - 1, m\right) - x\left(n - 1, m\right) y\left(n, m - 1\right) + x\left(n, m\right) y\left(n, m - 1\right) + x\left(n - 1, m\right) y\left(n, m\right) - x\left(n, m - 1\right) y\left(n, m\right)}\right)}^{2}\right)} {{\hbar}}^{2}}{8 \, {{\mu}}}

Numeric form.

eq20nm3=eq20nm0.subs(dict(zip(xyk,xyv)))
ru=fast_callable(eq20nm3.subs({hbar:0.2,mu:1.0}), vars=xyv, domain=RDF)

Classical force

F0(x)= 9*sech(5*x)^2*tanh(5*x)

plot(F0,(-3,3),axes_labels=['distance','force'],axes_labels_size=1)

Classical potential

P0(x)=9/(10*cosh(5*x)^2)

plot(P0,(-3,3),axes_labels=['distance','energy'],axes_labels_size=1,legend_label='potential')

Numeric form.

F1=fast_callable(F0(x11), vars=[x11], domain=RDF)
def FF(xy): FF1=F1(xy[0]-xy[1]); return array([FF1,-FF1])
def PO(xy): return P0(xy[0]-xy[1])

Initial Conditions

Consider two particles in a one dimensional space over 4 worlds. The initial configuration space is shown below.

N=2; M=2
r0=array([[[(i-0.5)/N,(j-0.5)/M] for j in range(-M/2+1,M/2+1)] for i in range(-N/2+1,N/2+1)])
p0=array(map(lambda y:map(lambda x:array([1.0*RDF(erfinv(2.0*x[0])-2.0+x[1]/2),RDF(erfinv(2.0*x[1])+2.0+x[0]/2)]),y),r0))
list_plot([r0[i,j] for i in range(N) for j in range(M)])+list_plot([p0[i,j] for i in range(N) for j in range(M)],color='red',figsize=[5,5])

The particles approach each other with the same initial velocity in each world.

v0=array([[[RDF(1.0),RDF(-1.0)] for j in range(M)] for i in range(N)])

In order to compute the quantum force on the particles in world n,mn,m we introduce ficticious worlds at "approximate infinity" to represent open boundry conditions.

inf = RDF(10.0^12);
inf2 = RDF(10.0^6);
def rf1(XY,n,m):
    global N,M
    p=[[XY[n+i,m+j][0],XY[n+i,m+j][1]] if n+i>=0 and m+j>=0 and n+i<N and m+j<M else [
      -inf if n+i<-1 else -inf2 if n+i<=-1 else inf if n+i>N else inf2 if n+i>=N else XY[n+i,m][0],
      -inf if m+j<-1 else -inf2 if m+j<=-1 else inf if m+j>M else inf2 if m+j>=M else XY[n,m+j][1]]
        for i in [-2,-1,0,1,2] for j in [-2,-1,0,1,2]]
    p=reduce(lambda x,y:x+y,p)
    r=ff(p)
    if isnan(r[0]): print "r[0] isnan"
    if isnan(r[1]): print "r[1] isnan"
    return r

Total classical potential

def PU(x): return sum([PO(x[i,j]) for i in range(N) for j in range(M)])

Total quantum Potential

def ru1(XY,n,m):
    global N,M
    nm=n+m+1
    p=[[XY[n+i,m+j][0],XY[n+i,m+j][1]] if n+i>=0 and m+j>=0 and n+i<N and m+j<M else [
      -inf if n+i<-1 else -inf2 if n+i<=-1 else inf if n+i>N else inf2 if n+i>=N else XY[n+i,m][0],
      -inf if m+j<-1 else -inf2 if m+j<=-1 else inf if m+j>M else inf2 if m+j>=M else XY[n,m+j][1]]
        for i in [-2,-1,0,1,2] for j in [-2,-1,0,1,2]]
    p=reduce(lambda x,y:x+y,p)
    return ru(*p)
def QU(x): return sum([ru1(x,i,j) for i in range(N) for j in range(M)])

Total kinetic energy

def K(x): return RDF(0.5)*sum([sum(x[i,j]^2) for i in range(N) for j in range(M)])

Simulation

Run the integration using Störmer–Verlet discretization with a continuous integrating stepsize controller.

%time
t_end = 6.0; step = .04
# initial step size
ds = RDF(0.001)
# controller
alpha = RDF(30.0); beta = RDF(0.0)
# step density
rho = RDF(1.0); dt=ds/rho
def G(p,q):
    a = sum([sum(p[i,j]*q[i,j]) for j in range(M) for i in range(N)])
    b = sum([sum(q[i,j]*q[i,j]) for j in range(M) for i in range(N)])
    eps = RDF(10.0^-10)
    return -alpha*a/max(b,eps)+beta*(1.0-ds/dt)
#acceleration = (quantum force + classical force)/mass
mm = RDF(1.0)
def A(x): return RDF(0.5)*array([[(rf1(x,i,j)+FF(XY[i,j]))/mm for j in range(M)] for i in range(N)])
# initial values
XY = p0.copy(); VV = v0.copy(); AX = A(XY)
rho = rho - RDF(0.5)*G(AX,VV)*dt
t = 0;  tt = [t]
XV = [ [ [ [XY[i,j][0],XY[i,j][1]] for i in range(N)] for j in range(M)] ]
KK=K(VV); PP=PU(XY); QQ=QU(XY)
T1=KK+PP+QQ
rhos = [rho]; TK={t:KK}; TP={t:PP};TQ={t:QQ}; T=[T1]
while t<t_end:
    t1 = t+step
    while t<t1:
        rho = rho + G(AX,VV)*dt
        if isnan(rho) or rho>100 or rho<0.01:
            print "Step control failed at %s. rho=%s"%(t,rho)
            t=t_end; break
        dt = ds/rho; t += dt
        # note 0.5 in A above
        XY += VV*dt + AX*dt^2
        VV += AX*dt; AX = A(XY); VV += AX*dt
    # Check total energy conservation
    KK=K(VV); PP=PU(XY); QQ=QU(XY)
    T2=KK+PP+QQ
    if abs(T1-T2)>0.1:
        print "Integration failed at %s. Try a shorter time step: %s."%(t,ds/2)
        #break
    T1=T2
    # record
    XV = XV + [ [ [ [XY[i,j][0],XY[i,j][1]] for i in range(N)] for j in range(M)] ]; tt += [t]
    rhos.append(rho); TK[t]=KK; TP[t]=PP; TQ[t]=QQ; T.append(T1)
CPU time: 33.48 s, Wall time: 41.55 s 
line2d(zip(tt,rhos),axes_labels=['time','rho'],legend_label='step density')
(line2d([[t,TK[t]] for t in tt],axes_labels=['time','energy'],legend_label='kinetic',color='green')+
 line2d([[t,TQ[t]] for t in tt],axes_labels=['time','energy'],legend_label='quantum potential',color='red')+
 line2d([[t,TP[t]] for t in tt],axes_labels=['time','energy'],legend_label='classical potential',color='blue'))
line2d(zip(tt,T),axes_labels=['time','energy'],legend_label='total',color='black')

line2d([[t,TQ[t]] for t in tt],axes_labels=['time','energy'],legend_label='quantum potential',color='red')

You might want something like this: https://addons.mozilla.org/En-US/firefox/addon/toggle-animated-gifs/ to manage the animation.

Start/stop GIF animations through a keyboard shortcut or by clicking them. You can also restart animations from the beginning, or disable animations by default. %md Each point in the following diagram represents the position of two particles. When the point approaches the x=y line

Each moving point in the following diagram represents the position of two particles in a single world (four worlds are shown). As the point approaches the x=y line (shown in red), the particles approach each other and the classical potential is at a maximum. In three worlds the particles overcome the potential barrier and continue. But one world (the last to approach) is reflected back. Although initially the two particles in each world approach each other with the same velocity, the particles in the last world to approach the barrier have lost energy to their counterparts in other worlds due to the quantum mechanical (interworld) interaction. When viewed in the projection we see two particles rebound instead of passing each other.

animate([text("time:"+("%4.1f"%x[0]).rjust(6), (6,3),color='black') +
text("kinetic:"+("%.4e"%TK[x[0]]).rjust(10), (6,2.5),color='black') +
text("classic:"+("%.4e"%TP[x[0]]).rjust(10), (6,2),color='black') +
text("quantum:"+("%.4e"%TQ[x[0]]).rjust(10), (6,1.5),color='black') +
line([[-6,-6],[6,6]],color='red',legend_label='potential barrier') +
  list_plot([x[1][i][j] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=-6,ymax=6) for x in zip(tt,XV)]).show(gif=true,delay=20)

show(sum([line([[tt[i],XV[i][j][k][0]] for i in range(0,len(tt))],color = 'red') for j in range(M) for k in range(N)])+
     sum([line([[tt[i],XV[i][j][k][1]] for i in range(0,len(tt))],color = 'blue') for j in range(M) for k in range(N)])
  ,figsize=[10,7],axes_labels=['time','location'],axes_labels_size=1)

animate([text("time:"+("%4.1f"%x[0]).rjust(5), (3,3),color='black') +
  list_plot([[x[1][i][j][0],PO(x[1][i][j])] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=0,ymax=3,color='red')+
  list_plot([[x[1][i][j][1],PO(x[1][i][j])] for j in range(N) for i in range(M)],xmin=-6,xmax=6,ymin=0,ymax=3,color='blue') for x in zip(tt,XV)]).show(gif=true,delay=20)