Better implementation for Laplace-Beltrami operator

edit

I implemented the Laplace-Beltrami operator on a scalar function psi as * d * d psi. This works correctly in the manifolds R², R⁴ and R⁶, but run-time increases super-exponentially with dimension. For R⁸, I gave up after 26 hours (and noticed that memory use was > 5.5 GB in core plus 4.5 GB swapped out). Below is the case of R⁶ for a known result: center-of-mass separation of two particles.

My question: is there a more efficient way to compute the Laplace-Beltrami operator on a scalar function?

sage: version()
'SageMath version 7.6, Release Date: 2017-03-25'
sage: Parallelism().set(nproc=1)
sage: M = Manifold(2*3,'R^6',field='real',start_index=1)
sage: m1, m2 = var('m1 m2', domain='real'); # Not domain='positive', avoiding Maxima bug
sage: #m1=1; m2=1; # for testing. Does impact execution time.
sage: m_CM = m1+m2; mu1 = m1/m_CM; mu2 = m2/m_CM; mu = m1*m2/m_CM

Remembering that the Schrödinger equation has terms $$\hbar / (2*m) \nabla \psi,$$ we adjust the metric accordingly.

sage: c_Cart.<x1,y1,z1,x2,y2,z2> = M.chart();
sage: g = M.riemannian_metric('g');
sage: g[1,1],g[2,2],g[3,3], g[4,4],g[5,5],g[6,6] = m1,m1,m1, m2,m2,m2; 
sage: g.display()
g = m1 dx1*dx1 + m1 dy1*dy1 + m1 dz1*dz1 + m2 dx2*dx2 + m2 dy2*dy2 + m2 dz2*dz2

Introduce center of mass coordinates, using total mass (m_CM) and reduced mass (mu).

sage: %time
sage: c_CM.<X,Y,Z,x,y,z> = M.chart();
sage: # Restrictions in following line seem to make no difference except for time.
sage: ch_Cart_CM = c_Cart.transition_map(c_CM, [mu1*x1+mu2*x2,mu1*y1+mu2*y2,mu1*z1+mu2*z2, x1-x2,y1-y2,z1-z2],  restrictions2 = [x!=0, y!=0, z!=0])
sage: ch_Cart_CM; ch_Cart_CM.inverse()
Change of coordinates from Chart (R^6, (x1, y1, z1, x2, y2, z2)) to Chart (R^6, (X, Y, Z, x, y, z))
Change of coordinates from Chart (R^6, (X, Y, Z, x, y, z)) to Chart (R^6, (x1, y1, z1, x2, y2, z2))
CPU time: 36.66 s,  Wall time: 36.67 s
sage: M.set_default_chart(c_CM)    # this saves a little typing later
sage: M.set_default_frame(c_CM.frame())

Wavefunction

sage: psiX = M.scalar_field({c_CM: function('Psi_X')(X,Y,Z)}, name='psiX',latex_name='\Psi_X'); psiX
Scalar field psiX on the 6-dimensional differentiable manifold R^6
sage: psix = M.scalar_field({c_CM: function('psi_x')(x,y,z)}, name='psix',latex_name='\psi_x'); psix
Scalar field psix on the 6-dimensional differentiable manifold R^6
sage: psi = psiX*psix

Here we build up the Laplace-Beltrami operator acting on ψ. We display each step, for pedagogical purposes.

sage: %time
sage: # The differential and the exterior derivative are equivalent when operating on a scalar function.
sage: dpsi = psi.differential(); dpsi; dpsi.parent(); #dpsi.display()
1-form dpsiX*psix on the 6-dimensional differentiable manifold R^6
Free module /\^1(R^6) of 1-forms on the 6-dimensional differentiable manifold R^6
CPU time: 0.28 s,  Wall time: 0.28 s
sage: %time
sage: sPsi = dpsi.hodge_dual(g); #sPsi.display()
CPU time: 135.07 s,  Wall time: 135.25 s
sage: %time
sage: dsPsi = sPsi.exterior_derivative(); #dsPsi.display()
CPU time: 4.72 s,  Wall time: 4.72 s
sage: %time
sage: LBpsi = dsPsi.hodge_dual(g); #LBpsi.display()
CPU time: 461.94 s,  Wall time: 461.96 s

 More briefly:

sage: %time
sage: LBpsi = psi.exterior_derivative().hodge_dual(g).exterior_derivative().hodge_dual(g); 
sage: #LBpsi.display(c_CM)
CPU time: 466.84 s,  Wall time: 466.85 s
sage: LBpsi.coord_function()
(m1*m2*(d^2(Psi_X)/dX^2 + d^2(Psi_X)/dY^2 + d^2(Psi_X)/dZ^2)*psi_x(x, y, z) + (m1^2*Psi_X(X, Y, Z) + 2*m1*m2*Psi_X(X, Y, Z) + m2^2*Psi_X(X, Y, Z))*d^2(psi_x)/dx^2 + (m1^2*Psi_X(X, Y, Z) + 2*m1*m2*Psi_X(X, Y, Z) + m2^2*Psi_X(X, Y, Z))*d^2(psi_x)/dy^2 + (m1^2*Psi_X(X, Y, Z) + 2*m1*m2*Psi_X(X, Y, Z) + m2^2*Psi_X(X, Y, Z))*d^2(psi_x)/dz^2)/(m1^2*m2 + m1*m2^2)

Although Sage has difficulty simplifying, we see that the above expression has the expected form, fully separating in center of mass coordinates.

sage: factor((m1^2*m2+m1*m2^2))
(m1 + m2)*m1*m2