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 with the Jupyter Notebook server, via the command sage -n jupyter
%display latex
Let us restore the 2-sphere manifold constructed in the preceeding worksheet:
M = Manifold(2, 'M')
U = M.open_subset('U')
XU.<x,y> = U.chart()
V = M.open_subset('V')
XV.<xp,yp> = V.chart("xp:x' yp:y'")
M.declare_union(U,V)
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)
XV_to_XU = XU_to_XV.inverse()
M.atlas()
We also reconstruct the point $p$:
p = U((1,2), chart=XU, name='p')
print(p)
Point p on the 2-dimensional differentiable manifold M
and the embedding $\mathbb{S}^2 \to \mathbb{R}^3$:
R3 = Manifold(3, 'R^3', r'\mathbb{R}^3')
XR3.<X,Y,Z> = R3.chart()
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()
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)
Finally we reconstruct the scalar field $f$:
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()
and assign the Python variable CM
to the algebra of scalar fields:
CM = M.scalar_field_algebra()
CM
The tangent space at the point $p$ introduced above is generated by
Tp = M.tangent_space(p)
Tp
It is a vector space over $\mathbb{R}$, which is represented by Sage's Symbolic Ring:
print(Tp.category())
Category of finite dimensional vector spaces over Symbolic Ring
The dimension of $T_p M$ is the same as that of $M$:
dim(Tp)
Tangent spaces are implemented as a class inherited from TangentSpace
via the category framework:
type(Tp)
The class TangentSpace
actually inherits from the generic class
FiniteRankFreeModule
, which, in SageMath, is devoted to free modules of finite rank
without any distinguished basis:
isinstance(Tp, FiniteRankFreeModule)
Two bases of $T_p M$ are already available: those generated by the derivations
at $p$ along the coordinates of charts XU
and XV
respectively:
Tp.bases()
None of these bases is distinguished, but one if the default one, which simply means that it is the basis to be considered if the basis argument is skipped in some methods:
Tp.default_basis()
A tangent vector is created as an element of the tangent space by the
standard SageMath procedure
new\_element = parent(...)
, where ...
stands for some material sufficient to construct the element:
vp = Tp((-3, 2), name='v')
print(vp)
Tangent vector v at Point p on the 2-dimensional differentiable manifold M
Since the basis is not specified, the pair $(-3,2)$ refers to components with respect to the default basis:
vp.display()
We have of course
vp.parent()
vp in Tp
As other manifold objects, tangent vectors have some plotting capabilities:
graph += vp.plot(chart=XR3, mapping=Phi, scale=0.5, color='gold')
show(graph, viewer='threejs', online=True)
The main attribute of the object vp
representing the vector $v$ is the dictionary _components
, which stores the components of $v$ in various bases of $T_p M$:
vp._components
As we can see above, the keys of the dictionary _components
are the bases of $T_p M$, while the values belongs to the class Components devoted to store ring elements indexed by integers or tuples of integers:
vpc = vp._components[Tp.default_basis()]
vpc
type(vpc)
The components themselves are stored in the dictionary _comp
of the Components
object, with the indices as keys:
vpc._comp
The $C^\infty(M)$-module of vector fields on $M$, $\mathfrak{X}(M)$, is obtained as
YM = M.vector_field_module()
YM
YM.category()
YM.base_ring() is CM
$\mathfrak{X}(M)$ is not a free module (at least its SageMath implementation does not
belong to the class FiniteRankFreeModule
):
isinstance(YM, FiniteRankFreeModule)
This is because $M=\mathbb{S}^2$ is not a parallelizable manifold:
M.is_manifestly_parallelizable()
Via the category mechanism,
the module $\mathfrak{X}(M)$ is implemented by a dynamically-generated subclass
of the class VectorFieldModule
, which is devoted to modules of vector fields
on non-parallelizable manifolds:
type(YM)
On the contrary, the set $\mathfrak{X}(U)$ of vector fields on $U$ is a free module of finite rank over the algebra $C^\infty(U)$:
YU = U.vector_field_module()
isinstance(YU, FiniteRankFreeModule)
YU.base_ring()
This is because the open subset $U$ is a parallelizable manifold:
U.is_manifestly_parallelizable()
being a coordinate chart domain:
U.is_manifestly_coordinate_domain()
We can check that in the atlas of $U$, at least one chart has $U$ for domain:
U.atlas()
The rank of $\mathfrak{X}(U)$ is the manifold's dimension:
rank(YU)
Via the category mechanism,
the free module $\mathfrak{X}(U)$ is implemented by a dynamically-generated subclass
of the class VectorFieldFreeModule
, which is devoted to modules of vector fields
on parallelizable manifolds:
type(YU)
The class VectorFieldFreeModule
is itself a subclass
of the generic class FiniteRankFreeModule
:
class_graph(
sage.manifolds.differentiable.vectorfield_module.VectorFieldFreeModule
).plot()
Since $U$ is a chart domain, the free module $\mathfrak{X}(U)$ is automatically endowed with a basis: the coordinate frame associated to the chart:
YU.bases()
Let us denote by eU
this frame. We can set eU = YU.bases()[0]
or
alternatively
eU = YU.default_basis()
eU
Another equivalent instruction would have been eU = U.default_frame()
.
Similarly, $\mathfrak{X}(V)$ is a free module, endowed with the coordinate frame
associated to stereographic coordinates from the South pole, which we
denote by eV
:
YV = V.vector_field_module()
YV.bases()
eV = YV.default_basis()
eV
If we consider the intersection $W=U\cap V$, we notice its module of vector fields is endowed with two bases, reflecting the fact that $W$ is covered by two charts: $(W,(x,y))$ and $(W,(x',y'))$:
W = U.intersection(V)
YW = W.vector_field_module()
YW.bases()
Let us denote by eUW
and eUV
these two bases, which are
actually the restrictions of the vector frames eU
and eV
to $W$:
eUW = eU.restrict(W)
eVW = eV.restrict(W)
YW.bases() == [eUW, eVW]
The free module $\mathfrak{X}(W)$ is also automatically endowed with automorphisms connecting the two bases, i.e. change-of-frame operators:
W.changes_of_frame()
The first of them is
P = W.change_of_frame(eUW, eVW)
P
It belongs to the general linear group of the free module $\mathfrak{X}(W)$:
P.parent()
and its matrix is deduced from the Jacobian matrix of the transition map XV
$\to$ XU
:
P[:]
We introduce a vector field $v$ on $M$ by
v = M.vector_field(name='v')
v[eU, 0] = f.restrict(U)
v[eU, 1] = -2
v.display(eU)
Notice that at this stage, we have defined $v$ only on $U$, by setting
its components in the vector frame eU
, either explicitely as scalar
fields, like the component $v^0$ set to the restriction of $f$ to $U$ or
to some symbolic expression, like for the component $v^1$: the $-2$
will be coerced to the constant scalar field of value $-2$.
We can ask for the scalar-field value of a components via the double-bracket
operator, since eU
is the default frame on $M$, we don't have to specify
it:
v[[0]]
v[[0]].display()
Note that the single bracket operator returns a chart function of the component:
v[0]
The restriction of $v$ to $W$ is of course
v.restrict(W).display(eUW)
Since we have a second vector frame on $W$, namely eVW
, and the
change-of-frame automorphisms are known, we can ask for the components
of $v$ with respect to that frame:
v.restrict(W).display(eVW)
Notice that the components are expressed in terms of the coordinates $(x,y)$
since they form the default chart on $W$. To have them expressed in
terms of the coordinates $(x',y')$, we have to add the restriction of
the chart
$(V,(x',y'))$ to $W$ as the second argument of the method display()
:
v.restrict(W).display(eVW, XV.restrict(W))
We extend the expression of $v$ to the full vector frame XV
by continuation of this expression:
v.add_comp_by_continuation(eV, W, chart=XV)
We have then
v.display(eV)
At this stage, the vector field $v$ is defined in all $M$. According to the hairy ball theorem\index{hairy ball theorem}, it has to vanish somewhere. Let us show that this occurs at the North pole, by first introducing the latter, as the point of stereographic coordinates $(x',y')=(0,0)$:
N = M((0,0), chart=XV, name='N')
print(N)
Point N on the 2-dimensional differentiable manifold M
As a check, we verify that the image of $N$ by the canonical embedding $\Phi: \mathbb{S}^2 \to \mathbb{R}^3$ is the point of Cartesian coordinates $(0,0,1)$:
XR3(Phi(N))
The vanishing of $\left. v\right| _N$:
v.at(N).display()
On the other hand, $v$ does not vanish at the point $p$ introduced above:
v.at(p).display()
We may plot the vector field $v$ in terms of the stereographic coordinates from the North pole:
v.plot(chart=XU, chart_domain=XU, max_range=2,
number_values=5, scale=0.4, aspect_ratio=1)
or in term of those from the South pole:
v.plot(chart=XV, chart_domain=XV, max_range=2,
number_values=9, scale=0.05, aspect_ratio=1)
Thanks to the embedding $\Phi$, we may also have a 3D plot of $v$ atop of the 3D plot already obtained:
graph_v = v.plot(chart=XR3, mapping=Phi, chart_domain=XU,
number_values=7, scale=0.2) + \
v.plot(chart=XR3, mapping=Phi, chart_domain=XV,
number_values=7, scale=0.2)
show(graph + graph_v, viewer='threejs', online=True)
Note that the sampling, performed on the two charts XU
and XV
is not uniform on the sphere. A better sampling would be acheived by introducing
spherical coordinates.
Vector fields on $M$ are implemented via the class VectorField
:
isinstance(v, sage.manifolds.differentiable.vectorfield.VectorField)
Since $M$ is not parallelizable, the representation of the vector field
$v$ via its restrictions to parallelizable open subsets; they are stored in the dictionary _restrictions
, whose keys are the open subsets:
v._restrictions
Let us consider one of these restrictions, for instance the restriction to $U$:
vU = v._restrictions[U]
vU is v.restrict(U)
Since $U$ is a parallelizable open subset, vU
belongs to the class
VectorFieldParal
:
isinstance(vU, sage.manifolds.differentiable.vectorfield.VectorFieldParal)
The main attribute of vU
is the dictionary _components
, which stores the components of vU
with respect to (possibly various) vector frames on $U$, the keys of that dictionary being the vector frames:
vU._components
v._restrictions[W]._components
The values of the dictionary _components
belong to the same class Components as that used for the storage of components of tangent vectors (cf. the example of vp
above):
vUc = vU._components[eU]
vUc
type(vUc)
As already mentioned above, the components themselves are stored in a dictionary, _comp
, whose keys are the indices:
vUc._comp
The difference with the tangent vector case is that the values are now scalar fields, i.e. elements of $C^\infty(U)$ in the present case. This is of course in agreement with the treatment of $\mathfrak{X}(U)$ as a free module over $C^\infty(U)$.
Let us perform some algebraic operation on $v$:
w = v + f*v
w
The code for the addition is accessible via
v.__add__??
As for the addition of scalar field (cf. this notebook), we see that __add__()
is implemented at the level
of the generic class Element
from which VectorField
inherits.
When both operands of the addition have the same parent, as here, __add__()
invokes the
method _add_()
(note the single underscore on each side of add
). This operator is
implemented at the level of TensorField
, as it can be checked from the source code:
v._add_??
The first step in the addition of two vector fields is to search in the
restrictions of both vector fields for common domains: this is performed via the method _common_subdomains
. Then the addition is performed at the level of the restrictions. The rest of the code is simply the set up of the vector field object containing the result.
Recursively, the addition performed will reach a level at which
the domains are parallelizable. Then a different method _add_()
, will
be involved, as we can check on vU
:
vU._add_??
We see that this method _add_()
is implemented
at the level of tensors on free modules, i.e. in the class
FreeModuleTensor, from which VectorFieldParal
inherits. Here the free module is clearly $\mathfrak{X}(U)$. The addition amounts to adding the components in a
a basis of the free module in which both operands have known components. Such
a basis is returned by the method common_basis
invoked in the first line.
If necessary, this method can use change-of-basis formulas to compute the
components of self
or other
in a common basis.
The addition of the components in the returned basis involves the method __add__()
of
class Components
; we can examine the corresponding code via the Components
object vUc
:
vUc.__add__??
We note that the computation can be parallelized on the
components if the user has turned on parallelization (this is done with the command
Parallelism().set(nproc=8)
(for 8 threads)). Focusing on the sequential code, we see that the addition is
performed component by component. Each component being an element of
$C^\infty(U)$ (the base ring of $\mathfrak{X}(U)$), this addition is that
of scalar fields, as presented in this notebook.
The action of $v$ on $f$ is defined pointwise by considering $v$ at each point $p\in M$ as a derivation (the very definition of a tangent vector); the result is then a scalar field $v(f)$ on $M$:
vf = v(f)
vf
vf.display()
df = f.differential()
df
print(df)
1-form df on the 2-dimensional differentiable manifold M
A 1-form is a tensor field of type $(0,1)$:
df.tensor_type()
while a vector field is a tensor field of type $(1,0)$:
v.tensor_type()
Specific 1-forms are those forming the dual basis (coframe) of a given vector frame: for instance for the vector frame eU
= $(U, (\partial_{x}, \partial_{y}))$, we have
eU.dual_basis()
print(eU.dual_basis()[0])
1-form dx on the Open subset U of the 2-dimensional differentiable manifold M
Since eU
is the default frame on $M$, the default display of $\mathrm{d}f$ is performed in terms of eU
's coframe:
df.display()
df[0] == diff(f.expr(), x)
df[1] == diff(f.expr(), y)
df.display(eV)
df[eV,0,XV]
df[eV,0,XV] == diff(f.expr(XV), xp)
df[eV,1,XV] == diff(f.expr(XV), yp)
The parent of $\mathrm{d}f$ is the set $\Omega^1(M)$ of all 1-forms on $M$, considered as a $C^\infty(M)$-module:
print(df.parent())
df.parent()
Module Omega^1(M) of 1-forms on the 2-dimensional differentiable manifold M
df.parent().base_ring()
This module is actually the dual of YM
= $\mathfrak{X}(M)$:
YM.dual()
A 1-form acts on vector fields, yielding a scalar field:
print(df(v))
Scalar field df(v) on the 2-dimensional differentiable manifold M
This scalar field is nothing but the result of the action of $v$ on $f$, considering $v$ at each point $p\in M$ as a derivation (the very definition of a tangent vector):
df(v) == v(f)
We construct a tensor of type $(1,1)$ by taking the tensor product $v\otimes \mathrm{d}f$:
t = v * df
t
t.display()
t.display(eV)
t.display_comp()
The parent of t
is the set $\mathcal{T}^{(1,1)}(M)$ of all type-$(1,1)$ tensor fields on $M$,
considered as a $C^\infty(M)$-module:
print(t.parent())
t.parent()
Module T^(1,1)(M) of type-(1,1) tensors fields on the 2-dimensional differentiable manifold M
t.parent().base_ring()
t._restrictions
As for vector fields, since $M$ is not parallelizable, the $C^\infty(M)$-module $\mathcal{T}^{(1,1)}(M)$ is not free and the tensor fields are described by their restrictions to parallelizable subdomains:
print(t._restrictions[U].parent())
Free module T^(1,1)(U) of type-(1,1) tensors fields on the Open subset U of the 2-dimensional differentiable manifold M
t._restrictions[U].parent().base_ring()
The standard metric on $M=\mathbb{S}^2$ is that induced by the Euclidean metric of $\mathbb{R}^3$. Let us start by defining the latter:
h = R3.metric('h')
h[0,0], h[1,1], h[2, 2] = 1, 1, 1
h.display()
The metric $g$ on $M$ is the pullback of $h$ associated with the embedding $\Phi$:
g = M.metric('g')
g.set( Phi.pullback(h) )
print(g)
Riemannian metric g on the 2-dimensional differentiable manifold M
Note that we could have defined $g$ intrinsically, i.e. by providing its components in the two frames eU
and eV
, as we did for the metric $h$ on $\mathbb{R}^3$. Instead, we have chosen to get it as the pullback of $h$, as an example of pullback associated with some differential map.
The metric is a symmetric tensor field of type (0,2):
g.tensor_type()
The expression of the metric in terms of the default frame on $M$ (eU
):
g.display()
We may factorize the metric components:
g[0,0].factor() ; g[1,1].factor()
g.display()
A matrix view of the components of $g$ in the manifold's default frame:
g[:]
Display in terms of the vector frame $(V, (\partial_{x'}, \partial_{y'}))$:
g.display(eV)
The metric acts on vector field pairs, resulting in a scalar field:
print(g(v,v))
Scalar field g(v,v) on the 2-dimensional differentiable manifold M
g(v,v).parent()
g(v,v).display()
The Levi-Civita connection associated with the metric $g$:
nab = g.connection()
print(nab)
nab
Levi-Civita connection nabla_g associated with the Riemannian metric g on the 2-dimensional differentiable manifold M
The nonzero Christoffel symbols of $g$ (skipping those that can be deduced by symmetry on the last two indices) w.r.t. the chart XU
:
g.christoffel_symbols_display(chart=XU)
$\nabla_g$ acting on the vector field $v$:
Dv = nab(v)
print(Dv)
Tensor field nabla_g(v) of type (1,1) on the 2-dimensional differentiable manifold M
Dv.display()
The Riemann tensor associated with the metric $g$:
Riem = g.riemann()
print(Riem)
Riem.display()
Tensor field Riem(g) of type (1,3) on the 2-dimensional differentiable manifold M
The components of the Riemann tensor in the default frame on $M$:
Riem.display_comp()
print(Riem.parent())
Module T^(1,3)(M) of type-(1,3) tensors fields on the 2-dimensional differentiable manifold M
Riem.symmetries()
no symmetry; antisymmetry: (2, 3)
The Riemann tensor associated with the Euclidean metric $h$ on $\mathbb{R}^3$ is identically zero:
h.riemann().display()
The Ricci tensor and the Ricci scalar:
Ric = g.ricci()
Ric.display()
R = g.ricci_scalar()
R.display()
Hence we recover the fact that $(\mathbb{S}^2,g)$ is a Riemannian manifold of constant positive curvature.
In dimension 2, the Riemann curvature tensor is entirely determined by the Ricci scalar $R$ according to $$ R^i_{\ \, jlk} = \frac{R}{2} \left( \delta^i_{\ \, k} g_{jl} - \delta^i_{\ \, l} g_{jk} \right)$$ Let us check this formula here, under the form $R^i_{\ \, jlk} = -R g_{j[k} \delta^i_{\ \, l]}$:
delta = M.tangent_identity_field()
Riem == - R*(g*delta).antisymmetrize(2,3)
Similarly the relation $\mathrm{Ric} = (R/2)\; g$ must hold:
Ric == (R/2)*g
The Levi-Civita tensor associated with $g$:
eps = g.volume_form()
print(eps)
eps.display()
2-form eps_g on the 2-dimensional differentiable manifold M
The exterior derivative of the 2-form $\epsilon_g$:
print(eps.exterior_derivative())
3-form deps_g on the 2-dimensional differentiable manifold M
Of course, since $M$ has dimension 2, all 3-forms vanish identically:
eps.exterior_derivative().display()