# Better implementation for Laplace-Beltrami operator

I implemented the Laplace-Beltrami operator on a scalar function psi as * d * d psi. This works correctly in the manifolds R2, R4 and R6, but run-time increases super-exponentially with dimension. For R8, 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 R6 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

edit retag close merge delete

Sort by » oldest newest most voted

For a scalar field, it is much faster to compute the Laplace-Beltrami operator as $\nabla_\mu \nabla^\mu\psi$, instead of $\star\mathrm{d}\star\mathrm{d}\psi$. For your example:

sage: nabla = g.connection()
sage: LBpsi = nabla(nabla(psi).up(g)).trace()


takes only 12 s on my computer (and of course yields the same result), see the online notebook.

The slowness you've experienced is due to the lack of optimization in the Hodge dual computation (especially for the Hodge dual of a 6-form, as indicated by the timing in cell 21 of the notebook), which is all the more severe that the dimension is high. This certainly will be improved in some future version. Thanks for pointing it out!

more