Ask Your Question

Better implementation for Laplace-Beltrami operator

asked 2017-05-04 23:12:35 +0200

Richard_L gravatar image

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())


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 flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted

answered 2017-05-05 14:20:47 +0200

eric_g gravatar image

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!

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2017-05-04 23:12:35 +0200

Seen: 384 times

Last updated: May 05 '17