Trac #30318: Dot and cross products along a differentiable map

This ticket makes the methods `dot_product()`, `cross_product()` and
`norm()` of class `VectorField` work for vector fields defined along a
differentiable map, the codomain of which is a Riemannian manifold.
Previously, these methods worked only for vector fields ''on'' a
Riemannian manifold, i.e. along the identity map.
An important subcase is of course that of a curve in a Riemannian
manifold.

For instance, considering a helix parametrized by its arc length:
{{{
sage: E.<x,y,z> = EuclideanSpace()
sage: R.<s> = RealLine()
sage: C = E.curve((2*cos(s/3), 2*sin(s/3), sqrt(5)*s/3), (s, 0, 6*pi),
....:             name='C')
}}}
we have now
{{{
sage: T = C.tangent_vector_field()
sage: T.display()
C' = -2/3*sin(1/3*s) e_x + 2/3*cos(1/3*s) e_y + 1/3*sqrt(5) e_z
sage: norm(T)
Scalar field |C'| on the Real interval (0, 6*pi)
sage: norm(T).expr()
1
}}}
Introducing the unit normal vector N via the derivative of T:
{{{
sage: I = C.domain()
sage: Tp = I.vector_field([diff(T[i], s) for i in E.irange()],
dest_map=C,
....:                     name="T'")
sage: N = Tp / norm(Tp)
}}}
we get the binormal vector as the cross product of T and N:
{{{
sage: B = T.cross_product(N)
sage: B
Vector field along the Real interval (0, 6*pi) with values on the
 Euclidean space E^3
sage: B.display()
1/3*sqrt(5)*sin(1/3*s) e_x - 1/3*sqrt(5)*cos(1/3*s) e_y + 2/3 e_z
}}}
We can then form the Frenet-Serret frame:
{{{
sage: FS = I.vector_frame(('T', 'N', 'B'), (T, N, B),
....:                     symbol_dual=('t', 'n', 'b'))
sage: FS
Vector frame ((0, 6*pi), (T,N,B)) with values on the Euclidean space E^3
}}}
and check that it is orthonormal:
{{{
sage: [[u.dot(v).expr() for v in FS] for u in FS]
[[1, 0, 0], [0, 1, 0], [0, 0, 1]]
}}}
The Frenet-Serret formulas are obtained as:
{{{
sage: Np = I.vector_field([diff(N[i], s) for i in E.irange()],
....:                     dest_map=C, name="N'")
sage: Bp = I.vector_field([diff(B[i], s) for i in E.irange()],
....:                     dest_map=C, name="B'")
sage: for v in (Tp, Np, Bp):
....:     v.display(FS)
....:
T' = 2/9 N
N' = -2/9 T + 1/9*sqrt(5) B
B' = -1/9*sqrt(5) N
}}}

URL: https://trac.sagemath.org/30318
Reported by: egourgoulhon
Ticket author(s): Eric Gourgoulhon
Reviewer(s): Travis Scrimshaw

src/sage/manifolds/differentiable/curve.py

@@ -245,6 +245,108 @@ class DifferentiableCurve(DiffMap):
         R --> R
            t |--> cos(t)^2
 
+    .. RUBRIC:: Curvature and torsion of a curve in a Riemannian manifold
+
+    Let us consider a helix `C` in the Euclidean space `\mathbb{E}^3`
+    parametrized by its arc length `s`::
+
+        sage: E.<x,y,z> = EuclideanSpace()
+        sage: R.<s> = RealLine()
+        sage: C = E.curve((2*cos(s/3), 2*sin(s/3), sqrt(5)*s/3), (s, 0, 6*pi),
+        ....:             name='C')
+
+    Its tangent vector field is::
+
+        sage: T = C.tangent_vector_field()
+        sage: T.display()
+        C' = -2/3*sin(1/3*s) e_x + 2/3*cos(1/3*s) e_y + 1/3*sqrt(5) e_z
+
+    Since `C` is parametrized by its arc length `s`, `T` is a unit vector (with
+    respect to the Euclidean metric of `\mathbb{E}^3`)::
+
+        sage: norm(T)
+        Scalar field |C'| on the Real interval (0, 6*pi)
+        sage: norm(T).display()
+        |C'|: (0, 6*pi) --> R
+           s |--> 1
+
+    Vector fields along `C` are defined by the method
+    :meth:`~sage.manifolds.differentiable.manifold.DifferentiableManifold.vector_field`
+    of the domain of `C` with the keyword argument ``dest_map`` set to `C`. For
+    instance the derivative vector `T'=\mathrm{d}T/\mathrm{d}s` is::
+
+        sage: I = C.domain(); I
+        Real interval (0, 6*pi)
+        sage: Tp = I.vector_field([diff(T[i], s) for i in E.irange()], dest_map=C,
+        ....:                     name="T'")
+        sage: Tp.display()
+        T' = -2/9*cos(1/3*s) e_x - 2/9*sin(1/3*s) e_y
+
+    The norm of `T'` is the curvature of the helix::
+
+        sage: kappa = norm(Tp)
+        sage: kappa
+        Scalar field |T'| on the Real interval (0, 6*pi)
+        sage: kappa.expr()
+        2/9
+
+    The unit normal vector along `C` is::
+
+        sage: N = Tp / kappa
+        sage: N.display()
+        -cos(1/3*s) e_x - sin(1/3*s) e_y
+
+    while the binormal vector along `C` is `B = T \times N`::
+
+        sage: B = T.cross_product(N)
+        sage: B
+        Vector field along the Real interval (0, 6*pi) with values on the
+         Euclidean space E^3
+        sage: B.display()
+        1/3*sqrt(5)*sin(1/3*s) e_x - 1/3*sqrt(5)*cos(1/3*s) e_y + 2/3 e_z
+
+    The three vector fields `(T, N, B)` form the **Frenet-Serret frame** along
+    `C`::
+
+        sage: FS = I.vector_frame(('T', 'N', 'B'), (T, N, B),
+        ....:                     symbol_dual=('t', 'n', 'b'))
+        sage: FS
+        Vector frame ((0, 6*pi), (T,N,B)) with values on the Euclidean space E^3
+
+    The Frenet-Serret frame is orthonormal::
+
+        sage: matrix([[u.dot(v).expr() for v in FS] for u in FS])
+        [1 0 0]
+        [0 1 0]
+        [0 0 1]
+
+    The derivative vectors `N'` and `B'`::
+
+        sage: Np = I.vector_field([diff(N[i], s) for i in E.irange()],
+        ....:                     dest_map=C, name="N'")
+        sage: Np.display()
+        N' = 1/3*sin(1/3*s) e_x - 1/3*cos(1/3*s) e_y
+        sage: Bp = I.vector_field([diff(B[i], s) for i in E.irange()],
+        ....:                     dest_map=C, name="B'")
+        sage: Bp.display()
+        B' = 1/9*sqrt(5)*cos(1/3*s) e_x + 1/9*sqrt(5)*sin(1/3*s) e_y
+
+    The Frenet-Serret formulas::
+
+        sage: for v in (Tp, Np, Bp):
+        ....:     v.display(FS)
+        ....:
+        T' = 2/9 N
+        N' = -2/9 T + 1/9*sqrt(5) B
+        B' = -1/9*sqrt(5) N
+
+    The torsion of `C` is obtained as the third component of `N'` in the
+    Frenet-Serret frame::
+
+        sage: tau = Np[FS, 3]
+        sage: tau
+        1/9*sqrt(5)
+
     """
     def __init__(self, parent, coord_expression=None, name=None,
                  latex_name=None, is_isomorphism=False, is_identity=False):

src/sage/manifolds/differentiable/pseudo_riemannian_submanifold.py

@@ -724,8 +724,7 @@ def calc_normal(chart):
             n_comp = (n_form.contract(*args) / factorial(self._dim)).contract(
                 self.ambient_metric().inverse().along(self._immersion))
             if self._ambient._dim - self._dim == 1:
-                n_comp = n_comp / n_comp.norm(
-                    self.ambient_metric().along(self._immersion))
+                n_comp = n_comp / n_comp.norm(self.ambient_metric())
 
             norm_rst = self._normal.restrict(chart.domain())
             norm_rst.add_comp(max_frame.restrict(chart.domain()))[:] = n_comp[:]

src/sage/manifolds/differentiable/vectorfield.py

@@ -1128,10 +1128,44 @@ def dot_product(self, other, metric=None):
             h(u,v): E^2 --> R
                (x, y) |--> -(x^3*y - x*y^3)/((x^2 + 1)*y^2 + x^2 + 1)
 
+        Scalar product of two vector fields along a curve (a lemniscate of
+        Gerono)::
+
+            sage: R.<t> = RealLine()
+            sage: C = M.curve([sin(t), sin(2*t)/2], (t, 0, 2*pi), name='C')
+            sage: u = C.tangent_vector_field(name='u')
+            sage: u.display()
+            u = cos(t) e_x + (2*cos(t)^2 - 1) e_y
+            sage: I = C.domain(); I
+            Real interval (0, 2*pi)
+            sage: v = I.vector_field(cos(t), -1, dest_map=C, name='v')
+            sage: v.display()
+            v = cos(t) e_x - e_y
+            sage: s = u.dot_product(v); s
+            Scalar field u.v on the Real interval (0, 2*pi)
+            sage: s.display()
+            u.v: (0, 2*pi) --> R
+               t |--> sin(t)^2
+
+        Scalar product between a vector field along the curve and a vector
+        field on the ambient Euclidean plane::
+
+            sage: e_x = M.cartesian_frame()[1]
+            sage: s = u.dot_product(e_x); s
+            Scalar field u.e_x on the Real interval (0, 2*pi)
+            sage: s.display()
+            u.e_x: (0, 2*pi) --> R
+               t |--> cos(t)
+
         """
         default_metric = metric is None
         if default_metric:
-            metric = self._domain.metric()
+            metric = self._ambient_domain.metric()
+        dest_map = self.parent().destination_map()
+        if dest_map != metric.parent().base_module().destination_map():
+            metric = metric.along(dest_map)
+        if dest_map != other.parent().destination_map():
+            other = other.along(dest_map)
         resu = metric(self, other)
         # From the above operation the name of resu is "g(u,v')" where
         # g = metric._name, u = self._name, v = other._name
@@ -1208,10 +1242,26 @@ def norm(self, metric=None):
             |v|_h: E^2 --> R
                (x, y) |--> sqrt((2*x^2 + 1)*y^2 + x^2)/(sqrt(x^2 + 1)*sqrt(y^2 + 1))
 
+        Norm of the tangent vector field to a curve (a lemniscate of Gerono)::
+
+            sage: R.<t> = RealLine()
+            sage: C = M.curve([sin(t), sin(2*t)/2], (t, 0, 2*pi), name='C')
+            sage: v = C.tangent_vector_field()
+            sage: v.display()
+            C' = cos(t) e_x + (2*cos(t)^2 - 1) e_y
+            sage: s = v.norm(); s
+            Scalar field |C'| on the Real interval (0, 2*pi)
+            sage: s.display()
+            |C'|: (0, 2*pi) --> R
+               t |--> sqrt(4*cos(t)^4 - 3*cos(t)^2 + 1)
+
         """
         default_metric = metric is None
         if default_metric:
-            metric = self._domain.metric()
+            metric = self._ambient_domain.metric()
+        dest_map = self.parent().destination_map()
+        if dest_map != metric.parent().base_module().destination_map():
+            metric = metric.along(dest_map)
         resu = metric(self, self).sqrt()
         if self._name is not None:
             if default_metric:
@@ -1299,14 +1349,50 @@ def cross_product(self, other, metric=None):
             sage: w.display()
             -(x^2 + y^2)*sqrt(x^2 + 1)/(sqrt(y^2 + 1)*sqrt(z^2 + 1)) e_z
 
+        Cross product of two vector fields along a curve (arc of a helix)::
+
+            sage: R.<t> = RealLine()
+            sage: C = M.curve((cos(t), sin(t), t), (t, 0, 2*pi), name='C')
+            sage: u = C.tangent_vector_field()
+            sage: u.display()
+            C' = -sin(t) e_x + cos(t) e_y + e_z
+            sage: I = C.domain(); I
+            Real interval (0, 2*pi)
+            sage: v = I.vector_field(-cos(t), sin(t), 0, dest_map=C)
+            sage: v.display()
+            -cos(t) e_x + sin(t) e_y
+            sage: w = u.cross_product(v); w
+            Vector field along the Real interval (0, 2*pi) with values on the
+             Euclidean space E^3
+            sage: w.parent().destination_map()
+            Curve C in the Euclidean space E^3
+            sage: w.display()
+            -sin(t) e_x - cos(t) e_y + (2*cos(t)^2 - 1) e_z
+
+        Cross product between a vector field along the curve and a vector field
+        on the ambient Euclidean space::
+
+            sage: e_x = M.cartesian_frame()[1]
+            sage: w = u.cross_product(e_x); w
+            Vector field C' x e_x along the Real interval (0, 2*pi) with values
+             on the Euclidean space E^3
+            sage: w.display()
+            C' x e_x = e_y - cos(t) e_z
+
         """
-        if self._domain.dim() != 3:
+        if self._ambient_domain.dim() != 3:
             raise ValueError("the cross product is not defined in dimension " +
                              "different from 3")
         default_metric = metric is None
         if default_metric:
-            metric = self._domain.metric()
-        eps = metric.volume_form(1)
+            metric = self._ambient_domain.metric()
+        dest_map = self.parent().destination_map()
+        if dest_map == metric.parent().base_module().destination_map():
+            eps = metric.volume_form(1)
+        else:
+            eps = metric.volume_form(1).along(dest_map)
+        if dest_map != other.parent().destination_map():
+            other = other.along(dest_map)
         resu = eps.contract(1, 2, self.wedge(other), 0, 1) / 2
         # The result is named "u x v" only for a default metric:
         if (default_metric and self._name is not None and