This notebook accompanies the lecture Symbolic tensor calculus on manifolds at JNCF 2018.
Click here to download the notebook file (ipynb format). To run it, you must start SageMath within the Jupyter Notebook server, via the command sage -n jupyter
NB: a version of SageMath at least equal to 7.5 (8.2 for the SymPy part) is required to run this notebook.
First we set up the notebook to use LaTeX display:
%display latex
Manifolds are constructed via the global function Manifold
:
M = Manifold(2, 'M')
print(M)
2-dimensional differentiable manifold M
By default, Manifold
returns a manifold over $\mathbb{R}$:
M.base_field()
Other base fields must be specified with the optional keyword field
, like
M = Manifold(2, 'M', field='complex')
We may check that $M$ is a topological space:
M in Sets().Topological()
Actually, $M$ belongs to the following categories:
M.categories()
As we can see, by default, Manifold
constructs a smooth manifold.
If one would like to stick to the topological level, one should write
M = Manifold(2, 'M', structure='topological')
Smooth manifolds are implemented by the Python class DifferentiableManifold
:
type(M)
The type of M
appears as DifferentiableManifold-with-category
because it is actually a subclass of DifferentiableManifold
, which is dynamically generated by SageMath's category mechanism:
isinstance(M,
sage.manifolds.differentiable.manifold.DifferentiableManifold)
The class DifferentiableManifold
inherits from TopologicalManifold
:
isinstance(M, sage.manifolds.manifold.TopologicalManifold)
We declare a chart, along with the symbols used to denote the coordinates (here $x=x^0$ and $y=x^1$) by
U = M.open_subset('U')
XU.<x,y> = U.chart()
XU
Open subsets are implemented by a (dynamically generated) subclass of
DifferentiableManifold
, since they are manifolds in their own:
isinstance(U,
sage.manifolds.differentiable.manifold.DifferentiableManifold)
Points on $M$ are created by passing their coordinates in a given chart:
p = U((1,2), chart=XU, name='p')
print(p)
Point p on the 2-dimensional differentiable manifold M
The syntax U(...)
used to create $p$ as an element of $U$ reflects the
parent/element pattern employed in SageMath; indeed $U$ is the parent of $p$:
p.parent()
Points are implemented by the class ManifoldPoint
.
The principal attribute of this class is a Python dictionary storing the point's coordinates
in various charts:
p._coordinates
Of course, we can recover the point's coordinates by letting the chart act on the point:
XU(p)
Let us introduce a second chart on $M$:
V = M.open_subset('V')
XV.<xp,yp> = V.chart("xp:x' yp:y'")
XV
and declare that $M$ is covered by only two charts, i.e. that $M=U\cup V$:
M.declare_union(U,V)
M.atlas()
We define the transition map XU
$\to$ XV
on $W = U\cap V$ as follows:
XU_to_XV = XU.transition_map(XV,
(x/(x^2+y^2), y/(x^2+y^2)),
intersection_name='W',
restrictions1= x^2+y^2!=0,
restrictions2= xp^2+yp^2!=0)
XU_to_XV.display()
The value of the argument restrictions1
means that
$W = U\setminus \{S\}$, where $S$ is the point of coordinates $(x,y)=(0,0)$,
while the value of restrictions2
means that
$W = V\setminus \{N\}$, where $N$ is the point of coordinates $(x',y')=(0,0)$.
Since $M=U\cup V$, we have then
$$
U = M \setminus \{N\},\qquad
V = M \setminus \{S\},\quad\mbox{and}\quad
W = M \setminus \{N, S\} .
$$
The transition map XV
$\to$ XU
is obtained by computing the inverse of the one defined above:
XU_to_XV.inverse().display()
At this stage, the smooth manifold $M$ is fully specified, being covered by
one atlas with all transition maps specified. One may have recognized that
$M$ is nothing but the 2-dimensional sphere:
$$
M = \mathbb{S}^2 ,
$$
with XU
(resp. XV
) being
the chart of \defin{stereographic coordinates}\index{stereographic!coordinates}
from the North pole $N$ (resp. the South pole $S$).
Since the transition maps have been defined, we can ask for the coordinates $(x',y')$ of the point $p$:
XV(p)
This operation has updated the dictionary _coordinates
:
p._coordinates
As a example of smooth map, let us consider the canonical embedding of $M=\mathbb{S}^2$ into $\mathbb{R}^3$. We need first to introduce $\mathbb{R}^3$ as a 3-dimensional smooth manifold, endowed with a single chart:
R3 = Manifold(3, 'R^3', r'\mathbb{R}^3')
XR3.<X,Y,Z> = R3.chart()
XR3
The embedding $\Phi: \mathbb{S}^2 \to \mathbb{R}^3$ is then defined in terms of its coordinate expression in the two charts covering $M=\mathbb{S}^2$:
Phi = M.diff_map(R3, {(XU, XR3):
[2*x/(1+x^2+y^2), 2*y/(1+x^2+y^2),
(x^2+y^2-1)/(1+x^2+y^2)],
(XV, XR3):
[2*xp/(1+xp^2+yp^2), 2*yp/(1+xp^2+yp^2),
(1-xp^2-yp^2)/(1+xp^2+yp^2)]},
name='Phi', latex_name=r'\Phi')
Phi.display()
We may use $\Phi$ for graphical purposes, for instance to display the grids of the stereographic charts XU
(in red) and XV
(in green), with the point $p$ atop:
graph = XU.plot(chart=XR3, mapping=Phi, number_values=25,
label_axes=False) + \
XV.plot(chart=XR3, mapping=Phi, number_values=25,
color='green', label_axes=False) + \
p.plot(chart=XR3, mapping=Phi, label_offset=0.05)
show(graph, viewer='threejs', online=True)
Let us define a scalar field, i.e. a smooth map $M\to \mathbb{R}$:
f = M.scalar_field({XU: 1/(1+x^2+y^2), XV: (xp^2+yp^2)/(1+xp^2+yp^2)},
name='f')
f.display()
Scalar fields are implemented by the class DiffScalarField
and the function chart representations are
stored in the attribute _express
of this class, which is a
Python dictionary whose keys are the various charts defined on $M$:
f._express
One may wonder about the compatibility of the two coordinate expressions
provided in the definition of $f$. Actually, to enforce the compatibility, it is possible to declare the scalar field in a single chart, XU
say,
and then to obtain its expression in chart XV
by analytic continuation
from the expression in $W=U\cap V$, where both expressions are known, thanks
to the transition map XV
$\to$ XU
:
f0 = M.scalar_field({XU: 1/(1+x^2+y^2)})
f0.add_expr_by_continuation(XV, U.intersection(V))
f == f0
The representation of the scalar field in a given chart, i.e. the public access to the directory _express
, is obtained via the method coord_function()
:
fU = f.coord_function(XU)
fU.display()
fV = f.coord_function(XV)
fV.display()
Both fU
and fV
are instances of the class ChartFunction
:
isinstance(fU, sage.manifolds.chart_func.ChartFunction)
Mathematically, chart functions are real-valued functions on the codomain of the considered chart. They map coordinates to elements of the base field (here $\mathbb{R}$):
fU(1,2)
fU(*XU(p))
while scalar fields maps manifold points to $\mathbb{R}$:
f(p)
Internally, each chart function stores coordinate expressions with respect to various computational engines:
The coordinate expressions are stored in the dictionary _express
,
whose keys are strings identifying the computational engines. By default
only SageMath symbolic expressions, i.e. expressions pertaining
to the so-called Symbolic Ring (SR
), are stored:
fU._express
The public access to the dictionary _express
is performed via the
method expr()
:
fU.expr()
type(fU.expr())
Actually, fU.expr()
is a shortcut for fU.expr('SR')
since
SR
is the default symbolic engine. Note that the class
Expression
is that devoted to SageMath symbolic expressions.
The method expr()
can also be invoked to get the expression in
another symbolic engine, for instance SymPy:
fU.expr('sympy')
type(fU.expr('sympy'))
This operation has updated the dictionary _express
:
fU._express
The default calculus engine for chart functions of chart XU
can
changed thanks to the method set_calculus_method()
:
XU.set_calculus_method('sympy')
fU.expr()
We can have better (i.e. LaTeX) rendering of SymPy objects with
from sympy import init_printing
init_printing()
We don't use it here to make clear the distinction between SR
objects and SymPy ones.
Reverting to SageMath's symbolic engine:
XU.set_calculus_method('SR')
fU.expr()
Symbolic expressions can be accessed directly from the scalar field,
f.expr(XU)
being a shortcut for f.coord_function(XU).expr()
:
f.expr(XU)
f.expr(XV)
The commutative algebra of scalar fields on $M$ is
CM = M.scalar_field_algebra()
CM
CM.category()
As for the manifold classes, the actual Python class implementing
$C^\infty(M)$ is inherited from DiffScalarFieldAlgebra
via the
category mechanism, hence it bares the name DiffScalarFieldAlgebra-with-category
:
type(CM)
The class DiffScalarFieldAlgebra-with-category
is dynamically generated as a subclass of DiffScalarFieldAlgebra
with extra functionalities, like
for instance the method is_commutative()
:
CM.is_commutative()
Let us have a look at the code of this method; it suffices to use the double question mark:
CM.is_commutative??
We see from the File
field that the code belongs to the category part of SageMath, not to the manifold part, where the class DiffScalarFieldAlgebra
is defined.
Regarding the scalar field $f$ introduced above, we have of course
f in CM
Actually, in Sage language, CM
=$C^\infty(M)$ is the parent of f
:
f.parent() is CM
The zero element of the algebra $C^\infty(M)$ is
CM.zero().display()
while its unit element is
CM.one().display()
Let us consider some operation in the algebra $C^\infty(M)$:
h = f + 2*CM.one()
h.display()
h(p)
We are going to investigate how the above addition is performed. For the Python interpreter
h = f + 2*CM.one()
is equivalent to h = f.__add__(2*CM.one())
,
i.e. the +
operator amounts to calling the method __add__()
on f
.
To have a look at the source code of this method, we use the double question mark:
f.__add__??
We see that the method __add__()
is implemented at the level
of the class Element
from which DiffScalarField
inherits, via
CommutativeAlgebraElement
.
In the present case, left
= f
and right
= 2*CM.one()
have the same parent, namely the algebra CM
, so that the actual
result is computed in the line
return (<Element>left)._add_(right)
This invokes the method _add_()
(note the single underscore on each side of add
). This operator is
implemented at the level of ScalarField
, as it can be checked from the source code:
f._add_??
This reflects a general strategy in SageMath: the arithmetic Python operators __add__
, __sub__
, etc. are implemented at the
top-level class Element
, while the specific element subclasses,
like ScalarField
here, implement single-underscore methods _add_
, _sub_
, etc., which perform the actual computation
when both operands have the same parent.
Looking at the code, we notice that the first step is to search
for the charts in which both operands of the addition operator have a coordinate
expression. This is performed by the method common_charts()
;
in the current case, we get the two stereographic charts defined on $M$:
f.common_charts(2*CM.one())
In general, common_charts()
returns the charts for which both operands
have already a known coordinate expression or for which a coordinate
expression can be computed by a known transition map, as we can see on the source code:
f.common_charts??
Once the list of charts in which both operands have a coordinate expression,
the addition is performed at the chart function level.
The code for the addition of chart functions defined on the same chart
is (recall that fU
is the chart function representing $f$ in chart XU
):
fU._add_??
We notice that the addition is performed on the symbolic expression
with respect to the symbolic engine currently at work (SageMath/Pynac, SymPy,...), as returned by
the method expr()
.
Let us recall that the user can change the symbolic engine at any time
by means of the method set_calculus_method()
, applied either to
a chart or to an open subset (possibly M
itself).
Besides, we notice in the code above that the result of the symbolic addition
is automatically simplified, by means of the method _simplify
.
It invokes a chain of simplification functions, which depends on the
symbolic engine (See here
for details).
Let us now discuss the second case in the __add__
method of
Element
, namely the case for which the parents of both operands are
different. This case is treated via SageMath coercion model, which allows one to deal with additions like
h1 = f + 2
h1.display()
A priori, f + 2
is not a well defined operation, since the integer $2$ does not
belong to the algebra $C^\infty(M)$. However SageMath manages to treat it
because $2$ can be coerced (i.e. automatically and unambigously converted) via CM(2)
into a element of $C^\infty(M)$, namely the constant scalar field whose value is $2$:
CM(2).display()
This happens because there exists a coercion map from the parent of $2$, namely the ring of integers $\mathbb{Z}$
(denoted ZZ
in SageMath), to $C^\infty(M)$:
2.parent()
CM.has_coerce_map_from(ZZ)