Richard_L's profile - activity

I understand how to compute the Levi-Civita connection on a pseudo-Riemannian manifold. I would like now to use it to compute the covariant derivative in various directions, as hinted in the code commentary to class Levi-Civita Connections. Is there a already existing method which I am overlooking or must I use explicit sums over coefficient indices?

Thanks. It is a little long-winded, but following eric_g's suggestion above I have nabla(psi).contract(Rho.frame()[1])

Covariant derivative I understand how to compute the Levi-Civita connection on a pseudo-Riemannian manifold. I would lik

Thank you. I missed that.

I create a subset of a Riemannian manifold. It is orientable, but it does not have the attribute 'volume_form'. Why? Note that the metric (g) has the attribute 'volume_form'. [Some variables, all real, elided below]

# Define manifold
M = Manifold(3,'R^3',field='real',start_index=1)
U = M.open_subset('U')

# Add chart and metric
Rho.<r12,r13,r23> = U.chart("r12:(0,+oo) r13:(0,+oo) r23:(0,+oo)")
Rho.add_restrictions([r23<r12+r13, r13<r12+r23, r12<r13+r23])
g = U.riemannian_metric('g');
g[:] = G[:].simplify_full()

# Establish affine connection
nabla = g.connection()

# Check that we have an orientation
U.orientation()

# Volume form 
#VolForm = U.volume_form()  # FAILS: 'DifferentiableManifold_with_category' object has no attribute 'volume_form'
VolForm = g.volume_form()  # This succeeds as expected
VolForm.display()

The on-line documentation shows examples of manifolds with attribute 'volume_form', thus my question.

Set up manifold with Riemann metric:

'SageMath version 7.6, Release Date: 2017-03-25'
sage: M = Manifold(2,'R^2',field='real',start_index=1)
sage: c_Cart.<x1,x2> = M.chart();
sage: g = M.riemannian_metric('g');
sage: g[1,1],g[2,2] = 1,1;
sage: g.display()
g = dx1dx1 + dx2dx2

And a change of coordinates:

sage: c_CM.<x,x> = M.chart();
sage: m1, m2 = var('m1 m2', domain='positive');
sage: m_CM = m1+m2; mu1 = m1/m_CM; mu2 = m2/m_CM; mu = m1m2/m_CM
sage: # Restrictions in following line seem to make no difference ...
sage: ch_Cart_CM = c_Cart.transition_map(c_CM, [mu1x1+mu2*x2, x1-x2], restrictions1 = x1!=x2, restrictions2 = x!=0)
sage: ch_Cart_CM;
Change of coordinates from Chart (R^2, (x1, x2)) to Chart (R^2, (X, x))

Now define a connection on g:

sage: g.comp(c_CM.frame())[:]
[ 2 -(m1 - m2)/(m1 + m2)]
[ -(m1 - m2)/(m1 + m2) (m1^2 + m2^2)/(m1^2 + 2m1m2 + m2^2)]

And a connection on G:

sage: G = M.riemannian_metric('G');
sage: G[:]= g.comp(c_CM.frame())[:]
sage: G.display()
G = 2 dx1dx1 - (m1 - m2)/(m1 + m2) dx1dx2 - (m1 - m2)/(m1 + m2) dx2dx1 + (m1^2 + m2^2)/(m1^2 + 2m1m2 + m2^2) dx2dx2

Everything appears to be fine. Let's try the Riemann connection:

sage: g.connection()
Levi-Civita connection nabla_g associated with the Riemannian metric g on the 2-dimensional differentiable manifold R^2

Good. BUT

sage: #G.connection() # blows up (in Maxima?) w/ "Detected access to protected memory ..."

Is there something wrong in the way G was defined, or are we limited by Maxima (whose ECL routines contain the message)?

I have a complicated metric tensor. In order to compute the Laplacian in this metric, I need the metric tensor inverse. But this exceeds my memory capacity after ~36 hours of computing.

I would like instead to set the inverse with what I believe to be the correct entries. Is that possible? Reviewing the documentation and the source code, I see no indication of such an operation.

I wonder if g.inverse().comp(XI.frame())[:] = Ginv[:] would work, where g denotes the metric tensor, and Ginv the matrix of expressions which should be correct.

BTW, I have the metric and its inverse in one frame. It is after mapping into a second frame [XI] that the problem arises.

I discovered on StackOverflow that some jupyter installations require specifying the ip address of localhost. So,

sage -n jupyter --ip=127.0.0.1

does the trick!

This seems to be a quirk of jupyter, since the classic sage notebook interface finds localhost just fine.

When I try to start sage with the Jupyter notebook, I get this error:

:~/work> sage -n jupyter
┌────────────────────────────────────────────────────────────────────┐
│ SageMath version 8.8, Release Date: 2019-06-26                     │
│ Using Python 2.7.15. Type "help()" for help.                       │
└────────────────────────────────────────────────────────────────────┘
Please wait while the Sage Jupyter Notebook server starts...
[I 13:23:55.592 NotebookApp] Using MathJax: nbextensions/mathjax/MathJax.js
[I 13:23:55.773 NotebookApp] Writing notebook server cookie secret to /run/user/1002/jupyter/notebook_cookie_secret
Traceback (most recent call last):
  File "/home/rllozes/sage-8.8/src/bin/sage-notebook", line 267, in <module>
    launcher(unknown)
  File "/home/rllozes/sage-8.8/src/bin/sage-notebook", line 99, in __init__
    main(argv)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/jupyter_core/application.py", line 266, in launch_instance
    return super(JupyterApp, cls).launch_instance(argv=argv, **kwargs)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/traitlets/config/application.py", line 657, in launch_instance
    app.initialize(argv)
  File "<decorator-gen-7>", line 2, in initialize
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/traitlets/config/application.py", line 87, in catch_config_error
    return method(app, *args, **kwargs)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/notebook/notebookapp.py", line 1635, in initialize
    self.init_webapp()
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/notebook/notebookapp.py", line 1407, in init_webapp
    self.http_server.listen(port, self.ip)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/tornado/tcpserver.py", line 142, in listen
    sockets = bind_sockets(port, address=address)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/site-packages/tornado/netutil.py", line 197, in bind_sockets
    sock.bind(sockaddr)
  File "/home/rllozes/sage-8.8/local/lib/python2.7/socket.py", line 230, in meth
    return getattr(self._sock,name)(*args)

socket.error: [Errno 99] Cannot assign requested address

Any idea why I get a socket error?

If I read the source code correctly, a matrix defined over SR is inverted by the Moore-Penrose pseudo-inversion. For a real, symmetric, symbolic, full-rank matrix, is this algorithm optimal? Or is there another algorithm that would be a) less memory intensive, and/or b) less time-consuming. [In my particular case, the attempted inversion of a 4x4 matrix runs out of memory (8 GB) after approximately 15 hours.]

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

Thank you, Eric.

In working with the Laplacian of a scalar field, I want to collect the complete coefficient expression of each partial derivative. In 3-space, this can be done "by hand", but in 6-space, automation is essential to avoid error. Consider:

sage: version()
'SageMath version 8.5, Release Date: 2018-12-22'
sage: #Patch to avoid maxima bug
sage: maxima_calculus.eval("domain:real;")
'real'
sage: M = Manifold(3,'R^3',field='real',start_index=1)
sage: c_C.<x,y,z> = M.chart();
sage: c_C.add_restrictions(z>0)
sage: g = M.riemannian_metric('g');
sage: g[1,1],g[2,2],g[3,3] = 1,1,1; 
sage: g.display()
g = dx*dx + dy*dy + dz*dz
sage: c_S.<r,theta,phi> = M.chart(r'r:(0,+oo) theta:(0,pi/2):\theta phi:(0,2*pi):\phi');
sage: ch_C_S = c_C.transition_map(c_S, (sqrt(x^2+y^2+z^2), arccos(z/(sqrt(x^2+y^2+z^2))), arctan2(y,x)))
sage: ch_C_S
Change of coordinates from Chart (R^3, (x, y, z)) to Chart (R^3, (r, theta, phi))
sage: ch_C_S.set_inverse(r*sin(theta)*cos(phi), r*sin(theta)*sin(phi), r*cos(theta))
sage: M.set_default_chart(c_S)    # this saves a little typing later
sage: M.set_default_frame(c_S.frame())
sage: chi = M.scalar_field(function('chi')(r,theta,phi), name='chi', latex_name=r'\chi')
sage: delChi = chi.laplacian(g)
sage: ddChi = delChi.expr()
sage: ddChi.collect_common_factors()
((r^2*diff(chi(r, theta, phi), r, r) + 2*r*diff(chi(r, theta, phi), r) + diff(chi(r, theta, phi), theta, theta))*sin(theta)^2 + cos(theta)*sin(theta)*diff(chi(r, theta, phi), theta) + diff(chi(r, theta, phi), phi, phi))/(r^2*sin(theta)^2)

# We can find the coefficients of scalar variables (sort of):
sage: ddChi.coefficient(r,2); ddChi.coefficient(r,1)
diff(chi(r, theta, phi), r, r)/r^2
2*diff(chi(r, theta, phi), r)/r^2
sage: ddChi.coefficient(sin(theta)^2)    # From which we see that .coefficient() is a shallow operator.
(r^2*diff(chi(r, theta, phi), r, r) + 2*r*diff(chi(r, theta, phi), r) + diff(chi(r, theta, phi), theta, theta))/(r^2*sin(theta)^2)

In the above, we work on the expression for the Laplacian (living in the symbolic ring). But .coefficient() will not work with a partial derivative as argument. Searching through the functions available for delChi (the Laplacian), I see nothing applicable there either. Any suggestions?

Thanks. Too bad inverse is not implemented for matrices over the (invertible) matrix ring.

Consider:

A = matrix(SR, 2, var('a1,b1,c1,d1'))
B = matrix(SR, 2, var('a2,b2,c2,d2'))
C = matrix(SR, 2, var('a3,b3,c3,d3'))
D = matrix(2,2, [A, B, B.T, C])
D; D.is_invertible()

This gives the matrix D as a matrix of (fully expanded) matrices, and confirms that D is invertible. However:

D.inverse()

results in

Traceback (click to the left of this block for traceback) 
...
AttributeError: 'MatrixSpace_with_category' object has no attribute
'is_field'

Seemingly, one cannot build after all a matrix of matrices, but only of basic ring elements.

Even better, I would prefer to construct a matrix of the (non-commutative, invertible) symbolic variables A, B, C, so as to eliminate the clutter of their constituent elements.

Nils Bruin provided the answer over on sage-devel. It is to insert

maxima_calculus.eval("domain:real;")

before the variable declarations.

Thanks,but that does not work. Sage will send to Maxima something like:

domain: complex;
declare( t, real );

Alas, certain expressions involving t will now fail. What I'm looking for is a way to cause Maxima to execute "domain:real;".