Trac #28462: Bug in computing the second fundamental form of a Riemannian submanifold

As reported in this [https://ask.sagemath.org/question/47786
ask.sagemath question], we have currently (Sage 8.8 and 8.9.beta9):
{{{
sage: P = Manifold(3, 'P', structure='Riemannian')
sage: Q = Manifold(2, 'Q', ambient=P, structure='Riemannian')
sage: CP.<x,y,z> = P.chart()
sage: CQ.<u,v> = Q.chart()
sage: g = P.metric()
sage: c = 2/(1 + y^2 + z^2)
sage: g[0,0], g[1,1], g[2,2] = 1, c^2, c^2
sage: phi = Q.diff_map(P, (u+v, u, v))
sage: phi_inv = P.diff_map(Q, (y, z))
sage: Q.set_embedding(phi, inverse=phi_inv)
sage: Q.second_fundamental_form()
TypeError: unable to convert 1/2*sqrt(2)*(u^4 + v^4 + 2*(u^2 + 1)*v^2 +
2*u^2 + 1)*y/
(sqrt(u^4 + v^4 + 2*(u^2 + 1)*v^2 + 2*u^2 + 3)*(y^2 + z^2 + 1))
+ 1/2*sqrt(2)*(u^4 + v^4 + 2*(u^2 + 1)*v^2 + 2*u^2 + 1)*z/
(sqrt(u^4 + v^4 + 2*(u^2 + 1)*v^2 + 2*u^2 + 3)*(y^2 + z^2 + 1))
to an integer
}}}
This results from the mention of the ring `SR` missing in the
declaration of a matrix involved in the computation.

The issue is fixed by this ticket.

URL: https://trac.sagemath.org/28462
Reported by: egourgoulhon
Ticket author(s): Eric Gourgoulhon
Reviewer(s): Florentin Jaffredo

Author: Release Manager

src/sage/manifolds/differentiable/pseudo_riemannian_submanifold.py

@@ -1058,6 +1058,23 @@ def second_fundamental_form(self):
             sage: N.extrinsic_curvature().display()  # long time
             K = -4/(x^4 + 8*x^2 + 16) dx*dx
 
+        An example with a non-Euclidean ambient metric::
+
+            sage: M = Manifold(2, 'M', structure='Riemannian')
+            sage: N = Manifold(1, 'N', ambient=M, structure='Riemannian',
+            ....:              start_index=1)
+            sage: CM.<x,y> = M.chart()
+            sage: CN.<u> = N.chart()
+            sage: g = M.metric()
+            sage: g[0, 0], g[1, 1] = 1, 1/(1 + y^2)^2
+            sage: phi = N.diff_map(M, (u, u))
+            sage: N.set_embedding(phi)
+            sage: N.second_fundamental_form()
+            Field of symmetric bilinear forms K on the 1-dimensional Riemannian
+             submanifold N embedded in the 2-dimensional Riemannian manifold M
+            sage: N.second_fundamental_form().display()
+            K = 2*sqrt(u^4 + 2*u^2 + 2)*u/(u^6 + 3*u^4 + 4*u^2 + 2) du*du
+
         """
         if self._ambient._dim - self._dim != 1:
             raise ValueError("second_fundamental_form is defined only for"
@@ -1080,15 +1097,19 @@ def second_fundamental_form(self):
                 n = self.normal()
 
                 for chart in self.atlas():
-                    gamma_n = matrix(self._dim + 1, self._dim + 1)
+                    gamma_n = matrix(SR, self._dim + 1, self._dim + 1)
+                    subs = dict(zip(self._ambient.default_chart()[:],
+                                    self._immersion.expression(chart)))
                     for i in range(self._dim + 1):
                         for j in range(self._dim + 1):
-                            gamma_n[i, j] = sum(
-                                nab[self._ambient.frames()[0], :][i][j][k].expr() *
+                            Gam_ij = [nab[self._ambient.frames()[0],
+                                          :][i][j][k].expr().subs(subs)
+                                      for k in range(self._dim + 1)]
+                            gamma_n[i, j] = chart.simplify(sum(
+                                Gam_ij[k] *
                                 n.restrict(chart.domain()).comp(
-                                    n.restrict(chart.domain())._fmodule.bases()[0])
-                                [:][k].expr() for k in
-                                range(self._dim + 1))
+                                  n.restrict(chart.domain())._fmodule.bases()[0])
+                                [:][k].expr() for k in range(self._dim + 1)))
                     dXdu = self._immersion.differential_functions(chart)
                     dNdu = matrix(SR, self._dim + 1, self._dim)
                     for i in range(self._dim + 1):
@@ -1100,7 +1121,11 @@ def second_fundamental_form(self):
                         self._immersion.restrict(chart.domain())).restrict(
                         chart.domain())[:, chart]
                     K = dXdu.transpose() * g * (dNdu + gamma_n * dXdu)
-                    resu[chart.frame(), :] = K
+                    si = self._sindex
+                    for i in self.irange():
+                        for j in self.irange(i):  # since K is symmetric
+                            resu[chart.frame(), i, j, chart] = chart.simplify(
+                                                      K[i - si, j - si].expr())
 
             self._second_fundamental_form = resu
         return self._second_fundamental_form