This notebook demonstrates a few differential geometry capabilities of SageMath on the example of the 2-dimensional sphere, using SymPy as symbolic engine, instead of SageMath's default one (Pynac+Maxima). The relevent tools have been implemented through the SageManifolds project.

Click here to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, via the command `sage -n jupyter`

*NB:* a version of SageMath at least equal to 8.3 is required to run this worksheet:

In [1]:

```
version()
```

Out[1]:

First we set up the notebook to display mathematical objects using LaTeX rendering:

In [2]:

```
%display latex
from sympy import init_printing
init_printing()
```

We also define a viewer for 3D plots (use `'threejs'`

or `'jmol'`

for interactive 3D graphics):

In [3]:

```
viewer3D = 'threejs' # must be 'threejs', jmol', 'tachyon' or None (default)
```

We start by declaring $\mathbb{S}^2$ as a differentiable manifold of dimension 2 over $\mathbb{R}$:

In [4]:

```
S2 = Manifold(2, 'S^2', latex_name=r'\mathbb{S}^2', start_index=1)
```

The first argument, 2, is the dimension of the manifold, while the second argument is the symbol used to label the manifold.

The argument start_index sets the index range to be used on the manifold for labelling components w.r.t. a basis or a frame: start_index=1 corresponds to $\{1,2\}$; the default value is start_index=0 and yields to $\{0,1\}$.

In [5]:

```
print(S2)
```

In [6]:

```
S2
```

Out[6]:

The manifold is a `Parent`

object:

In [7]:

```
isinstance(S2, Parent)
```

Out[7]:

in the category of smooth manifolds over $\mathbb{R}$:

In [8]:

```
S2.category()
```

Out[8]:

We ask for all symbolic computations to be performed with SymPy, instead of SageMath's default symbolic engine (Pynac+Maxima, implemented via SageMath's Symbolic Ring):

In [9]:

```
S2.set_calculus_method('sympy')
```

The sphere cannot be covered by a single chart. At least two charts are necessary, for instance the charts associated with the stereographic projections from the North pole and the South pole respectively. Let us introduce the open subsets covered by these two charts:
$$ U := \mathbb{S}^2\setminus\{N\}, $$
$$ V := \mathbb{S}^2\setminus\{S\}, $$
where $N$ is a point of $\mathbb{S}^2$, which we shall call the *North pole*, and $S$ is the point of $U$ of stereographic coordinates $(0,0)$, which we call the *South pole*:

In [10]:

```
U = S2.open_subset('U') ; print(U)
```

In [11]:

```
V = S2.open_subset('V') ; print(V)
```

We declare that $\mathbb{S}^2 = U \cup V$:

In [12]:

```
S2.declare_union(U, V)
```

Then we declare the stereographic chart on $U$, denoting by $(x,y)$ the coordinates resulting from the stereographic projection from the North pole:

In [13]:

```
stereoN.<x,y> = U.chart()
```

The expression `.<x,y>`

in the left-hand side means that the Python variables `x`

and `y`

are set to the two coordinates of the chart. This allows one to refer subsequently to the coordinates by their names. In the present case, the function `chart()`

has no argument, which implies that the coordinate symbols will be `x`

and `y`

(i.e. exactly the characters appearing in the `<...>`

operator) and that each coordinate range is $(-\infty,+\infty)$. As we will see below, for other cases, an argument must be passed to `chart()`

to specify each coordinate symbol and range, as well as some specific LaTeX symbol.

In [14]:

```
stereoN
```

Out[14]:

The coordinates can be accessed individually, either by means of their indices in the chart ( following the convention `start_index=1`

set in the manifold's definition) or by their names as Python variables:

In [15]:

```
stereoN[1]
```

Out[15]:

In [16]:

```
y is stereoN[2]
```

Out[16]:

Similarly, we introduce on $V$ the coordinates $(x',y')$ corresponding to the stereographic projection from the South pole:

In [17]:

```
stereoS.<xp,yp> = V.chart("xp:x' yp:y'")
```

In this case, the string argument passed to `chart`

stipulates that the text-only names of the coordinates are xp and yp (same as the Python variables names defined within the `<...>`

operator in the left-hand side), while their LaTeX names are $x'$ and $y'$.

In [18]:

```
stereoS
```

Out[18]:

At this stage, the user's atlas on the manifold has two charts:

In [19]:

```
S2.atlas()
```

Out[19]:

We have to specify the **transition map** between the charts 'stereoN' = $(U,(x,y))$ and 'stereoS' = $(V,(x',y'))$; it is given by the standard inversion formulas:

In [20]:

```
stereoN_to_S = stereoN.transition_map(stereoS,
(x/(x^2+y^2), y/(x^2+y^2)),
intersection_name='W',
restrictions1= x^2+y^2!=0,
restrictions2= xp^2+xp^2!=0)
stereoN_to_S.display()
```

Out[20]:

In the above declaration, 'W' is the name given to the chart-overlap subset: $W := U\cap V$, the condition $x^2+y^2 \not=0$ defines $W$ as a subset of $U$, and the condition $x'^2+y'^2\not=0$ defines $W$ as a subset of $V$.

The inverse coordinate transformation is computed by means of the method `inverse()`

:

In [21]:

```
stereoS_to_N = stereoN_to_S.inverse()
stereoS_to_N.display()
```

Out[21]:

In the present case, the situation is of course perfectly symmetric regarding the coordinates $(x,y)$ and $(x',y')$.

At this stage, the user's atlas has four charts:

In [22]:

```
S2.atlas()
```

Out[22]:

Let us store $W = U\cap V$ into a Python variable for future use:

In [23]:

```
W = U.intersection(V)
```

Similarly we store the charts $(W,(x,y))$ (the restriction of $(U,(x,y))$ to $W$) and $(W,(x',y'))$ (the restriction of $(V,(x',y'))$ to $W$) into Python variables:

In [24]:

```
stereoN_W = stereoN.restrict(W)
stereoN_W
```

Out[24]:

In [25]:

```
stereoS_W = stereoS.restrict(W)
stereoS_W
```

Out[25]:

We may plot the chart $(W, (x',y'))$ in terms of itself, as a grid:

In [26]:

```
stereoS_W.plot()
```

Out[26]:

More interestingly, let us plot the stereographic chart $(x',y')$ in terms of the stereographic chart $(x,y)$ on the domain $W$ where both systems overlap (we split the plot in four parts to avoid the singularity at $(x',y')=(0,0)$):

In [27]:

```
graphSN1 = stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[-6,-0.02]})
graphSN2 = stereoS_W.plot(stereoN, ranges={xp:[-6,-0.02], yp:[0.02,6]})
graphSN3 = stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[-6,-0.02]})
graphSN4 = stereoS_W.plot(stereoN, ranges={xp:[0.02,6], yp:[0.02,6]})
show(graphSN1+graphSN2+graphSN3+graphSN4,
xmin=-1.5, xmax=1.5, ymin=-1.5, ymax=1.5)
```

The standard spherical (or polar) coordinates $(\theta,\phi)$ are defined on the open domain $A\subset W \subset \mathbb{S}^2$ that is the complement of the "origin meridian"; since the latter is the half-circle defined by $y=0$ and $x\geq 0$, we declare:

In [28]:

```
A = W.open_subset('A', coord_def={stereoN_W: (y!=0, x<0),
stereoS_W: (yp!=0, xp<0)})
print(A)
```

The restriction of the stereographic chart from the North pole to $A$ is

In [29]:

```
stereoN_A = stereoN_W.restrict(A)
stereoN_A
```

Out[29]:

We then declare the chart $(A,(\theta,\phi))$ by specifying the intervals $(0,\pi)$ and $(0,2\pi)$ spanned by respectively $\theta$ and $\phi$:

In [30]:

```
spher.<th,ph> = A.chart(r'theta:(0,pi):\theta phi:(0,2*pi):\phi') ; spher
```

Out[30]:

The specification of the spherical coordinates is completed by providing the transition map with the stereographic chart $(A,(x,y))$:

In [31]:

```
spher_to_stereoN = spher.transition_map(stereoN_A,
(sin(th)*cos(ph)/(1-cos(th)),
sin(th)*sin(ph)/(1-cos(th))))
spher_to_stereoN.display()
```

Out[31]:

We also provide the inverse transition map, asking to check that the provided formulas are indeed correct (argument `verbose=True`

):

In [32]:

```
spher_to_stereoN.set_inverse(2*atan(1/sqrt(x^2+y^2)), atan2(-y,-x)+pi,
verbose=True)
```

The check is passed, modulo some lack of trigonometric simplifications in the first two lines.

In [33]:

```
spher_to_stereoN.inverse().display()
```

Out[33]:

The transition map $(A,(\theta,\phi))\rightarrow (A,(x',y'))$ is obtained by combining the transition maps $(A,(\theta,\phi))\rightarrow (A,(x,y))$ and $(A,(x,y))\rightarrow (A,(x',y'))$:

In [34]:

```
stereoN_to_S_A = stereoN_to_S.restrict(A)
spher_to_stereoS = stereoN_to_S_A * spher_to_stereoN
spher_to_stereoS.display()
```

Out[34]:

Similarly, the transition map $(A,(x',y'))\rightarrow (A,(\theta,\phi))$ is obtained by combining the transition maps $(A,(x',y'))\rightarrow (A,(x,y))$ and $(A,(x,y))\rightarrow (A,(\theta,\phi))$:

In [35]:

```
stereoS_to_N_A = stereoN_to_S.inverse().restrict(A)
stereoS_to_spher = spher_to_stereoN.inverse() * stereoS_to_N_A
stereoS_to_spher.display()
```

Out[35]:

The user atlas of $\mathbb{S}^2$ is now

In [36]:

```
S2.atlas()
```

Out[36]:

Let us draw the grid of spherical coordinates $(\theta,\phi)$ in terms of stereographic coordinates from the North pole $(x,y)$:

In [37]:

```
spher.plot(stereoN, number_values=15, ranges={th: (pi/8,pi)})
```

Out[37]:

Conversly, we may represent the grid of the stereographic coordinates $(x,y)$ restricted to $A$ in terms of the spherical coordinates $(\theta,\phi)$. We limit ourselves to one quarter (cf. the argument ranges):

In [38]:

```
stereoN_A.plot(spher, ranges={x: (0.01,8), y: (0.01,8)}, number_values=20, plot_points=200)
```

Out[38]:

We declare the **North pole** (resp. the **South pole**) as the point of coordinates $(0,0)$ in the chart $(V,(x',y'))$ (resp. in the chart $(U,(x,y))$):

In [39]:

```
N = V.point((0,0), chart=stereoS, name='N') ; print(N)
S = U.point((0,0), chart=stereoN, name='S') ; print(S)
```

Since points are Sage Element's, the corresponding Parent being the manifold subsets, an equivalent writing of the above declarations is

In [40]:

```
N = V((0,0), chart=stereoS, name='N') ; print(N)
S = U((0,0), chart=stereoN, name='S') ; print(S)
```

Moreover, since stereoS in the default chart on $V$ and stereoN is the default one on $U$, their mentions can be omitted, so that the above can be shortened to

In [41]:

```
N = V((0,0), name='N') ; print(N)
S = U((0,0), name='S') ; print(S)
```

In [42]:

```
N.parent()
```

Out[42]:

In [43]:

```
S.parent()
```

Out[43]:

We have of course

In [44]:

```
N in V
```

Out[44]:

In [45]:

```
N in S2
```

Out[45]:

In [46]:

```
N in U
```

Out[46]:

In [47]:

```
N in A
```

Out[47]:

Let us introduce some point at the equator:

In [48]:

```
E = S2((0,1), chart=stereoN, name='E')
```

The point $E$ is in the open subset $A$:

In [49]:

```
E in A
```

Out[49]:

We may then ask for its spherical coordinates $(\theta,\phi)$:

In [50]:

```
E.coord(spher)
```

Out[50]:

which is not possible for the point $N$:

In [51]:

```
try:
N.coord(spher)
except ValueError as exc:
print('Error: ' + str(exc))
```

Let us first declare $\mathbb{R}^3$ as a 3-dimensional manifold covered by a single chart (the so-called** Cartesian coordinates**):

In [52]:

```
R3 = Manifold(3, 'R^3', r'\mathbb{R}^3', start_index=1)
R3.set_calculus_method('sympy')
cart.<X,Y,Z> = R3.chart() ; cart
```

Out[52]:

The embedding of the sphere is defined as a differential mapping $\Phi: \mathbb{S}^2 \rightarrow \mathbb{R}^3$:

In [53]:

```
Phi = S2.diff_map(R3, {(stereoN, cart):
[2*x/(1+x^2+y^2), 2*y/(1+x^2+y^2),
(x^2+y^2-1)/(1+x^2+y^2)],
(stereoS, cart):
[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')
```

In [54]:

```
Phi.display()
```

Out[54]:

In [55]:

```
Phi.parent()
```

Out[55]:

In [56]:

```
print(Phi.parent())
```

In [57]:

```
Phi.parent() is Hom(S2, R3)
```

Out[57]:

$\Phi$ maps points of $\mathbb{S}^2$ to points of $\mathbb{R}^3$:

In [58]:

```
N1 = Phi(N) ; print(N1) ; N1 ; N1.coord()
```

Out[58]:

In [59]:

```
S1 = Phi(S) ; print(S1) ; S1 ; S1.coord()
```

Out[59]:

In [60]:

```
E1 = Phi(E) ; print(E1) ; E1 ; E1.coord()
```

Out[60]:

$\Phi$ has been defined in terms of the stereographic charts $(U,(x,y))$ and $(V,(x',y'))$, but we may ask its expression in terms of spherical coordinates. The latter is then computed by means of the transition map $(A,(x,y))\rightarrow (A,(\theta,\phi))$:

In [61]:

```
Phi.expr(stereoN_A, cart)
```

Out[61]:

In [62]:

```
Phi.expr(spher, cart)
```

Out[62]:

In [63]:

```
Phi.display(spher, cart)
```

Out[63]:

Let us use $\Phi$ to draw the grid of spherical coordinates $(\theta,\phi)$ in terms of the Cartesian coordinates $(X,Y,Z)$ of $\mathbb{R}^3$:

In [64]:

```
graph_spher = spher.plot(chart=cart, mapping=Phi, number_values=11,
color='blue', label_axes=False)
show(graph_spher, viewer=viewer3D)
```

We may also use the embedding $\Phi$ to display the stereographic coordinate grid in terms of the Cartesian coordinates in $\mathbb{R}^3$. First for the stereographic coordinates from the North pole:

In [65]:

```
graph_stereoN = stereoN.plot(chart=cart, mapping=Phi, number_values=25,
label_axes=False)
show(graph_stereoN, viewer=viewer3D)
```

and then have a view with the stereographic coordinates from the South pole superposed (in green):

In [66]:

```
graph_stereoS = stereoS.plot(chart=cart, mapping=Phi, number_values=25,
color='green', label_axes=False)
show(graph_stereoN + graph_stereoS, viewer=viewer3D)
```