Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

Laplace-Beltrami misbehaving

Two misbehaviors are apparent in the following: 1) objects not displaying in declared default frame, and 2) unexpected zero after applying Laplace-Beltrami operator to a scalar function in one of the two charts/frames. Q: Am I setting up the charts and frames incorrectly? Q: Missing something else? (I have heavily elided the sagemath output.) Thank you.

Preliminaries:

sage: version()
'SageMath version 8.2.beta7, Release Date: 2018-03-03'
sage: Parallelism().set(nproc=2)
sage: Manifold.options.omit_function_arguments=True
sage: from sage.manifolds.utilities import ExpressionNice
sage: var('rho12,rho13,rho23', domain='real')
sage: assume(rho12>0, rho13>0, rho23>0)
sage: var('r12,r13,r23', domain='real')
sage: assume(r12>0, r13>0, r23>0, r23<r12+r13, r13<r12+r23, r12<r13+r23)
sage: var('m1 m2 m3', domain='real')
sage: var('mu12,mu13,mu23', domain='real')
sage: assume(m1>0, m2>0, m3>0)
sage: m1=1; m2=1; m3=1
sage: mu12 = (m1*m2)/(m1+m2); mu23 = (m2*m3)/(m2+m3); mu13 = (m1*m3)/(m1+m3)
sage: rho12 = r12^2; rho13 = r13^2; rho23 = r23^2

sage: Ginv = matrix([[2*rho12/mu12, 1/m1*(rho12+rho13-rho23), 1/m2*(rho12+rho23-rho13)],\
...           [1/m1*(rho12+rho13-rho23), 2*rho13/mu13, 1/m3*(rho13+rho23-rho12)],\
...           [1/m2*(rho12+rho23-rho13), 1/m3*(rho13+rho23-rho12), 2*rho23/mu23]])
...
sage: Ginv
[               4*r12^2  r12^2 + r13^2 - r23^2  r12^2 - r13^2 + r23^2]
[ r12^2 + r13^2 - r23^2                4*r13^2 -r12^2 + r13^2 + r23^2]
[ r12^2 - r13^2 + r23^2 -r12^2 + r13^2 + r23^2                4*r23^2]

sage: G = Ginv.inverse(); #G.simplify_full()

Establish manifold, Riemann metric, charts and frames...

sage: M = Manifold(3,'R^3',field='real',start_index=1)
sage: U = M.open_subset('U')
sage: Rho.<r12,r13,r23> = U.chart("r12:(0,+oo) r13:(0,+oo) r23:(0,+oo)")
sage: Rho.add_restrictions([r23<r12+r13, r13<r12+r23, r12<r13+r23])
sage: g = M.riemannian_metric('g');
sage: g[:] = G[:].simplify_full()
sage: ginv = g.inverse(); ginv.display()
inv_g = 4*r12^2 d/dr12*d/dr12 + (r12^2 + r13^2 - r23^2) d/dr12*d/dr13 + (r12^2 - r13^2 + r23^2) d/dr12*d/dr23 + (r12^2 + r13^2 - r23^2) d/dr13*d/dr12 + 4*r13^2 d/dr13*d/dr13 + (-r12^2 + r13^2 + r23^2) d/dr13*d/dr23 + (r12^2 - r13^2 + r23^2) d/dr23*d/dr12 + (-r12^2 + r13^2 + r23^2) d/dr23*d/dr13 + 4*r23^2 d/dr23*d/dr23

sage: Rho.domain(); Rho.parent()
Open subset U of the 3-dimensional differentiable manifold R^3
<class 'sage.manifolds.differentiable.chart.RealDiffChart'>

sage: # Pekeris coordinates
sage: Tau.<u1,u2,u3> = U.chart("u1:(0,+oo) u2:(0,+oo) u3:(0,+oo)")
sage: Rho_Tau = Rho.transition_map(Tau, ((r12+r13-r23)/2, (r12-r13+r23)/2, (-r12+r13+r23)/2))
sage: Tau_Rho = Rho_Tau.inverse()
sage: U.set_default_chart(Tau)
sage: U.set_default_frame(Tau.frame())
sage: Rho_Tau.display(); Tau_Rho.display();
u1 = 1/2*r12 + 1/2*r13 - 1/2*r23
u2 = 1/2*r12 - 1/2*r13 + 1/2*r23
u3 = -1/2*r12 + 1/2*r13 + 1/2*r23
r12 = u1 + u2
r13 = u1 + u3
r23 = u2 + u3

sage: U.atlas(); Rho.frame(); Tau.frame()
[Chart (U, (r12, r13, r23)), Chart (U, (u1, u2, u3))]
Coordinate frame (U, (d/dr12,d/dr13,d/dr23))
Coordinate frame (U, (d/du1,d/du2,d/du3))

sage: g.display_comp(frame=Tau.frame())
g_u1,u1 = 1/48*(3*u2^4 + 12*u2*u3^3 + 3*u3^4 + 2*(8*u1*u2 + 9*u2^2)*u3^2 + 4*(4*u1^2*u2 + 4*u1*u2^2 + 3*u2^3)*u3)/(u1*u2*u3^4 + 2*(u1^2*u2 + u1*u2^2)*u3^3 + (2*u1^3*u2 + 3*u1^2*u2^2 + 2*u1*u2^3)*u3^2 + (u1^4*u2 + 2*u1^3*u2^2 + 2*u1^2*u2^3 + u1*u2^4)*u3) 
... [elided]
g_u3,u3 = 1/48*(3*u1^4 + 12*u1^3*u2 + 18*u1^2*u2^2 + 12*u1*u2^3 + 3*u2^4 + 16*u1*u2*u3^2 + 16*(u1^2*u2 + u1*u2^2)*u3)/(u1*u2*u3^4 + 2*(u1^2*u2 + u1*u2^2)*u3^3 + (2*u1^3*u2 + 3*u1^2*u2^2 + 2*u1*u2^3)*u3^2 + (u1^4*u2 + 2*u1^3*u2^2 + 2*u1^2*u2^3 + u1*u2^4)*u3)

First problem:

sage: #Why does this cell not compute and display in the Tau frame (declared default above)?
sage: ginv = g.inverse(); ginv.display_comp()
inv_g^r12,r12 = 4*u1^2 + 8*u1*u2 + 4*u2^2 
inv_g^r12,r13 = 2*u1^2 + 2*u1*u2 + 2*(u1 - u2)*u3 
inv_g^r12,r23 = 2*u1*u2 + 2*u2^2 - 2*(u1 - u2)*u3 
inv_g^r13,r12 = 2*u1^2 + 2*u1*u2 + 2*(u1 - u2)*u3 
inv_g^r13,r13 = 4*u1^2 + 8*u1*u3 + 4*u3^2 
inv_g^r13,r23 = -2*u1*u2 + 2*(u1 + u2)*u3 + 2*u3^2 
inv_g^r23,r12 = 2*u1*u2 + 2*u2^2 - 2*(u1 - u2)*u3 
inv_g^r23,r13 = -2*u1*u2 + 2*(u1 + u2)*u3 + 2*u3^2 
inv_g^r23,r23 = 4*u2^2 + 8*u2*u3 + 4*u3^2 

sage: #Display in Tau frame can be forced.
sage: ginv.display_comp(frame=Tau.frame())
inv_g^u1,u1 = 3*u1^2 + 3*u1*u2 + u2^2 + (3*u1 - u2)*u3 + u3^2 
inv_g^u1,u2 = u1*u2 - (u1 + u2)*u3 - u3^2 
inv_g^u1,u3 = -u1*u2 - u2^2 + (u1 - u2)*u3 
inv_g^u2,u1 = u1*u2 - (u1 + u2)*u3 - u3^2 
inv_g^u2,u2 = u1^2 + 3*u1*u2 + 3*u2^2 - (u1 - 3*u2)*u3 + u3^2 
inv_g^u2,u3 = -u1^2 - u1*u2 - (u1 - u2)*u3 
inv_g^u3,u1 = -u1*u2 - u2^2 + (u1 - u2)*u3 
inv_g^u3,u2 = -u1^2 - u1*u2 - (u1 - u2)*u3 
inv_g^u3,u3 = u1^2 - u1*u2 + u2^2 + 3*(u1 + u2)*u3 + 3*u3^2 

sage: nabla = g.connection()
sage: # nabla also not displayed in Tau frame!
sage: nabla.display()
Gam^r12_r12,r12 = -1/6*(8*r12^11 - r12^10*r13 - 8*r12^9*r13^2 + 3*r12^8*r13^3 - 24*r12^7*r13^4 - 2*r12^6*r13^5 + 40*r12^5*r13^6 - 2*r12^4*r13^7 - 16*r12^3*r13^8 + 3*r12^2*r13^9 - r13^11 + 13*r13*r23^10 - r23^11 + (3*r12^2 + 29*r13^2)*r23^9 - (16*r12^3 + 29*r12^2*r13 + 96*r12*r13^2 + 17*r13^3)*r23^8 - 2*(r12^4 - 4*r12^2*r13^2 + 23*r13^4)*r23^7 + 2*(20*r12^5 + 5*r12^4*r13 + 32*r12^3*r13^2 + 12*r12^2*r13^3 + 48*r12*r13^4 + 11*r13^5)*r23^6 - 2*(r12^6 + 7*r12^4*r13^2 + 3*r12^2*r13^4 - 11*r13^6)*r23^5 - 2*(12*r12^7 - 7*r12^6*r13 - 28*r12^5*r13^2 - 3*r12^4*r13^3 - 144*r12^3*r13^4 + 3*r12^2*r13^5 - 48*r12*r13^6 + 23*r13^7)*r23^4 + (3*r12^8 - 16*r12^6*r13^2 + 6*r12^4*r13^4 + 24*r12^2*r13^6 - 17*r13^8)*r23^3 - (8*r12^9 + 7*r12^8*r13 + 16*r12^7*r13^2 + 16*r12^6*r13^3 - 56*r12^5*r13^4 + 14*r12^4*r13^5 - 64*r12^3*r13^6 - 8*r12^2*r13^7 + 96*r12*r13^8 - 29*r13^9)*r23^2 - (r12^10 + 7*r12^8*r13^2 - 14*r12^6*r13^4 - 10*r12^4*r13^6 + 29*r12^2*r13^8 - 13*r13^10)*r23)/(r12^12 - 2*r12^10*r13^2 - r12^8*r13^4 + 4*r12^6*r13^6 - r12^4*r13^8 - 2*r12^2*r13^10 + r13^12 + r23^12 - 2*(r12^2 + r13^2)*r23^10 - (r12^4 + 10*r12^2*r13^2 + r13^4)*r23^8 + 4*(r12^6 + 3*r12^4*r13^2 + 3*r12^2*r13^4 + r13^6)*r23^6 - (r12^8 - 12*r12^6*r13^2 - 42*r12^4*r13^4 - 12*r12^2*r13^6 + r13^8)*r23^4 - 2*(r12^10 + 5*r12^8*r13^2 - 6*r12^6*r13^4 - 6*r12^4*r13^6 + 5*r12^2*r13^8 + r13^10)*r23^2) 
... [elided]

Second problem: [Appears to me to be incomplete conversion to the Tau frame.]

sage: psi = U.scalar_field({Rho: function('Psi')(r12,r13,r23)}, name='psi',latex_name='\Psi'); psi
Scalar field psi on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: LBpsi = nabla(nabla(psi).up(g)).trace(); LBpsi
Scalar field on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: DelPsi = LBpsi.expr().factor(); ExpressionNice(DelPsi)
1/2*(8*u1^6*u2*u3*d^2(Psi)/d(u1 + u2)^2 + 32*u1^5*u2^2*u3*d^2(Psi)/d(u1 + u2)^2 + 56*u1^4*u2^3*u3*d^2(Psi)/d(u1 + u2)^2 + 56*u1^3*u2^4*u3*d^2(Psi)/d(u1 + u2)^2 + 32*u1^2*u2^5*u3*d^2(Psi)/d(u1 + u2)^2 + 8*u1*u2^6*u3*d^2(Psi)/d(u1 + u2)^2 + 16*u1^5*u2*u3^2*d^2(Psi)/d(u1 + u2)^2 + 56*u1^4*u2^2*u3^2*d^2(Psi)/d(u1 + u2)^2 + 
... [elided]
u3) + 4*u2^2*u3^5*d(Psi)/d(u2 + u3) - u1*u3^6*d(Psi)/d(u2 + u3) + u2*u3^6*d(Psi)/d(u2 + u3))/((u1^2 + u1*u2 + u2^2 + u1*u3 + u2*u3 + u3^2)*(u1 + u2 + u3)*u1*u2*u3)

Leading to the unexpected zero:

sage: psi = U.scalar_field({Rho: function('Psi')(u1,u2,u3)}, name='psi',latex_name='\Psi'); psi
Scalar field psi on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: LBpsi = nabla(nabla(psi).up(g)).trace(); LBpsi
Scalar field zero on the Open subset U of the 3-dimensional differentiable manifold R^3

Laplace-Beltrami misbehaving

Two misbehaviors are apparent in the following: 1) objects not displaying in declared default frame, and 2) unexpected zero after applying Laplace-Beltrami operator to a scalar function in one of the two charts/frames. Q: Am I setting up the charts and frames incorrectly? Q: Missing something else? (I have heavily elided the sagemath output.) Thank you.

Preliminaries:

sage: version()
'SageMath version 8.2.beta7, Release Date: 2018-03-03'
sage: Parallelism().set(nproc=2)
sage: Manifold.options.omit_function_arguments=True
sage: from sage.manifolds.utilities import ExpressionNice
sage: var('rho12,rho13,rho23', domain='real')
sage: assume(rho12>0, rho13>0, rho23>0)
sage: var('r12,r13,r23', domain='real')
sage: assume(r12>0, r13>0, r23>0, r23<r12+r13, r13<r12+r23, r12<r13+r23)
sage: var('m1 m2 m3', domain='real')
sage: var('mu12,mu13,mu23', domain='real')
sage: assume(m1>0, m2>0, m3>0)
sage: m1=1; m2=1; m3=1
sage: mu12 = (m1*m2)/(m1+m2); mu23 = (m2*m3)/(m2+m3); mu13 = (m1*m3)/(m1+m3)
sage: rho12 = r12^2; rho13 = r13^2; rho23 = r23^2

sage: Ginv = matrix([[2*rho12/mu12, 1/m1*(rho12+rho13-rho23), 1/m2*(rho12+rho23-rho13)],\
...           [1/m1*(rho12+rho13-rho23), 2*rho13/mu13, 1/m3*(rho13+rho23-rho12)],\
...           [1/m2*(rho12+rho23-rho13), 1/m3*(rho13+rho23-rho12), 2*rho23/mu23]])
...
sage: Ginv
[               4*r12^2  r12^2 + r13^2 - r23^2  r12^2 - r13^2 + r23^2]
[ r12^2 + r13^2 - r23^2                4*r13^2 -r12^2 + r13^2 + r23^2]
[ r12^2 - r13^2 + r23^2 -r12^2 + r13^2 + r23^2                4*r23^2]

sage: G = Ginv.inverse(); #G.simplify_full()

Establish manifold, Riemann metric, charts and frames...

sage: M = Manifold(3,'R^3',field='real',start_index=1)
sage: U = M.open_subset('U')
sage: Rho.<r12,r13,r23> = U.chart("r12:(0,+oo) r13:(0,+oo) r23:(0,+oo)")
sage: Rho.add_restrictions([r23<r12+r13, r13<r12+r23, r12<r13+r23])
sage: g = M.riemannian_metric('g');
sage: g[:] = G[:].simplify_full()
sage: ginv = g.inverse(); ginv.display()
inv_g = 4*r12^2 d/dr12*d/dr12 + (r12^2 + r13^2 - r23^2) d/dr12*d/dr13 + (r12^2 - r13^2 + r23^2) d/dr12*d/dr23 + (r12^2 + r13^2 - r23^2) d/dr13*d/dr12 + 4*r13^2 d/dr13*d/dr13 + (-r12^2 + r13^2 + r23^2) d/dr13*d/dr23 + (r12^2 - r13^2 + r23^2) d/dr23*d/dr12 + (-r12^2 + r13^2 + r23^2) d/dr23*d/dr13 + 4*r23^2 d/dr23*d/dr23

sage: Rho.domain(); Rho.parent()
Open subset U of the 3-dimensional differentiable manifold R^3
<class 'sage.manifolds.differentiable.chart.RealDiffChart'>

sage: # Pekeris coordinates
sage: Tau.<u1,u2,u3> = U.chart("u1:(0,+oo) u2:(0,+oo) u3:(0,+oo)")
sage: Rho_Tau = Rho.transition_map(Tau, ((r12+r13-r23)/2, (r12-r13+r23)/2, (-r12+r13+r23)/2))
sage: Tau_Rho = Rho_Tau.inverse()
sage: U.set_default_chart(Tau)
sage: U.set_default_frame(Tau.frame())
sage: Rho_Tau.display(); Tau_Rho.display();
u1 = 1/2*r12 + 1/2*r13 - 1/2*r23
u2 = 1/2*r12 - 1/2*r13 + 1/2*r23
u3 = -1/2*r12 + 1/2*r13 + 1/2*r23
r12 = u1 + u2
r13 = u1 + u3
r23 = u2 + u3

sage: U.atlas(); Rho.frame(); Tau.frame()
[Chart (U, (r12, r13, r23)), Chart (U, (u1, u2, u3))]
Coordinate frame (U, (d/dr12,d/dr13,d/dr23))
Coordinate frame (U, (d/du1,d/du2,d/du3))

sage: g.display_comp(frame=Tau.frame())
g_u1,u1 = 1/48*(3*u2^4 + 12*u2*u3^3 + 3*u3^4 + 2*(8*u1*u2 + 9*u2^2)*u3^2 + 4*(4*u1^2*u2 + 4*u1*u2^2 + 3*u2^3)*u3)/(u1*u2*u3^4 + 2*(u1^2*u2 + u1*u2^2)*u3^3 + (2*u1^3*u2 + 3*u1^2*u2^2 + 2*u1*u2^3)*u3^2 + (u1^4*u2 + 2*u1^3*u2^2 + 2*u1^2*u2^3 + u1*u2^4)*u3) 
... [elided]
g_u3,u3 = 1/48*(3*u1^4 + 12*u1^3*u2 + 18*u1^2*u2^2 + 12*u1*u2^3 + 3*u2^4 + 16*u1*u2*u3^2 + 16*(u1^2*u2 + u1*u2^2)*u3)/(u1*u2*u3^4 + 2*(u1^2*u2 + u1*u2^2)*u3^3 + (2*u1^3*u2 + 3*u1^2*u2^2 + 2*u1*u2^3)*u3^2 + (u1^4*u2 + 2*u1^3*u2^2 + 2*u1^2*u2^3 + u1*u2^4)*u3)

First problem:

sage: #Why does this cell not compute and display in the Tau frame (declared default above)?
sage: ginv = g.inverse(); ginv.display_comp()
inv_g^r12,r12 = 4*u1^2 + 8*u1*u2 + 4*u2^2 
inv_g^r12,r13 = 2*u1^2 + 2*u1*u2 + 2*(u1 - u2)*u3 
inv_g^r12,r23 = 2*u1*u2 + 2*u2^2 - 2*(u1 - u2)*u3 
inv_g^r13,r12 = 2*u1^2 + 2*u1*u2 + 2*(u1 - u2)*u3 
inv_g^r13,r13 = 4*u1^2 + 8*u1*u3 + 4*u3^2 
inv_g^r13,r23 = -2*u1*u2 + 2*(u1 + u2)*u3 + 2*u3^2 
inv_g^r23,r12 = 2*u1*u2 + 2*u2^2 - 2*(u1 - u2)*u3 
inv_g^r23,r13 = -2*u1*u2 + 2*(u1 + u2)*u3 + 2*u3^2 
inv_g^r23,r23 = 4*u2^2 + 8*u2*u3 + 4*u3^2 

sage: #Display in Tau frame can be forced.
sage: ginv.display_comp(frame=Tau.frame())
inv_g^u1,u1 = 3*u1^2 + 3*u1*u2 + u2^2 + (3*u1 - u2)*u3 + u3^2 
inv_g^u1,u2 = u1*u2 - (u1 + u2)*u3 - u3^2 
inv_g^u1,u3 = -u1*u2 - u2^2 + (u1 - u2)*u3 
inv_g^u2,u1 = u1*u2 - (u1 + u2)*u3 - u3^2 
inv_g^u2,u2 = u1^2 + 3*u1*u2 + 3*u2^2 - (u1 - 3*u2)*u3 + u3^2 
inv_g^u2,u3 = -u1^2 - u1*u2 - (u1 - u2)*u3 
inv_g^u3,u1 = -u1*u2 - u2^2 + (u1 - u2)*u3 
inv_g^u3,u2 = -u1^2 - u1*u2 - (u1 - u2)*u3 
inv_g^u3,u3 = u1^2 - u1*u2 + u2^2 + 3*(u1 + u2)*u3 + 3*u3^2 

sage: nabla = g.connection()
sage: # nabla also not displayed in Tau frame!
sage: nabla.display()
Gam^r12_r12,r12 = -1/6*(8*r12^11 - r12^10*r13 - 8*r12^9*r13^2 + 3*r12^8*r13^3 - 24*r12^7*r13^4 - 2*r12^6*r13^5 + 40*r12^5*r13^6 - 2*r12^4*r13^7 - 16*r12^3*r13^8 + 3*r12^2*r13^9 - r13^11 + 13*r13*r23^10 - r23^11 + (3*r12^2 + 29*r13^2)*r23^9 - (16*r12^3 + 29*r12^2*r13 + 96*r12*r13^2 + 17*r13^3)*r23^8 - 2*(r12^4 - 4*r12^2*r13^2 + 23*r13^4)*r23^7 + 2*(20*r12^5 + 5*r12^4*r13 + 32*r12^3*r13^2 + 12*r12^2*r13^3 + 48*r12*r13^4 + 11*r13^5)*r23^6 - 2*(r12^6 + 7*r12^4*r13^2 + 3*r12^2*r13^4 - 11*r13^6)*r23^5 - 2*(12*r12^7 - 7*r12^6*r13 - 28*r12^5*r13^2 - 3*r12^4*r13^3 - 144*r12^3*r13^4 + 3*r12^2*r13^5 - 48*r12*r13^6 + 23*r13^7)*r23^4 + (3*r12^8 - 16*r12^6*r13^2 + 6*r12^4*r13^4 + 24*r12^2*r13^6 - 17*r13^8)*r23^3 - (8*r12^9 + 7*r12^8*r13 + 16*r12^7*r13^2 + 16*r12^6*r13^3 - 56*r12^5*r13^4 + 14*r12^4*r13^5 - 64*r12^3*r13^6 - 8*r12^2*r13^7 + 96*r12*r13^8 - 29*r13^9)*r23^2 - (r12^10 + 7*r12^8*r13^2 - 14*r12^6*r13^4 - 10*r12^4*r13^6 + 29*r12^2*r13^8 - 13*r13^10)*r23)/(r12^12 - 2*r12^10*r13^2 - r12^8*r13^4 + 4*r12^6*r13^6 - r12^4*r13^8 - 2*r12^2*r13^10 + r13^12 + r23^12 - 2*(r12^2 + r13^2)*r23^10 - (r12^4 + 10*r12^2*r13^2 + r13^4)*r23^8 + 4*(r12^6 + 3*r12^4*r13^2 + 3*r12^2*r13^4 + r13^6)*r23^6 - (r12^8 - 12*r12^6*r13^2 - 42*r12^4*r13^4 - 12*r12^2*r13^6 + r13^8)*r23^4 - 2*(r12^10 + 5*r12^8*r13^2 - 6*r12^6*r13^4 - 6*r12^4*r13^6 + 5*r12^2*r13^8 + r13^10)*r23^2) 
... [elided]

Second problem: [Appears to me to be incomplete conversion to the Tau frame.]

sage: psi = U.scalar_field({Rho: function('Psi')(r12,r13,r23)}, name='psi',latex_name='\Psi'); psi
Scalar field psi on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: LBpsi = nabla(nabla(psi).up(g)).trace(); LBpsi
Scalar field on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: DelPsi = LBpsi.expr().factor(); ExpressionNice(DelPsi)
1/2*(8*u1^6*u2*u3*d^2(Psi)/d(u1 + u2)^2 + 32*u1^5*u2^2*u3*d^2(Psi)/d(u1 + u2)^2 + 56*u1^4*u2^3*u3*d^2(Psi)/d(u1 + u2)^2 + 56*u1^3*u2^4*u3*d^2(Psi)/d(u1 + u2)^2 + 32*u1^2*u2^5*u3*d^2(Psi)/d(u1 + u2)^2 + 8*u1*u2^6*u3*d^2(Psi)/d(u1 + u2)^2 + 16*u1^5*u2*u3^2*d^2(Psi)/d(u1 + u2)^2 + 56*u1^4*u2^2*u3^2*d^2(Psi)/d(u1 + u2)^2 + 
... [elided]
u3) + 4*u2^2*u3^5*d(Psi)/d(u2 + u3) - u1*u3^6*d(Psi)/d(u2 + u3) + u2*u3^6*d(Psi)/d(u2 + u3))/((u1^2 + u1*u2 + u2^2 + u1*u3 + u2*u3 + u3^2)*(u1 + u2 + u3)*u1*u2*u3)

Leading to the unexpected zero:

sage: psi = U.scalar_field({Rho: function('Psi')(u1,u2,u3)}, name='psi',latex_name='\Psi'); psi
Scalar field psi on the Open subset U of the 3-dimensional differentiable manifold R^3

sage: LBpsi = nabla(nabla(psi).up(g)).trace(); LBpsi
Scalar field zero on the Open subset U of the 3-dimensional differentiable manifold R^3