This notebook demonstrates some differential geometry capabilities of SageMath on the example of the 3-dimensional sphere, $\mathbb{S}^3$. The corresponding tools have been developed within the SageManifolds project (version 1.3, as included in SageMath 8.3).
Click here to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, via the command sage -n jupyter
NB: a version of SageMath at least equal to 7.5 is required to run this worksheet:
version()
'SageMath version 8.3, Release Date: 2018-08-03'
First we set up the notebook to display mathematical objects using LaTeX formatting:
%display latex
We also define a viewer for 3D plots (use 'threejs'
or 'jmol'
for interactive 3D graphics):
viewer3D = 'threejs' # must be 'threejs', jmol', 'tachyon' or None (default)
To increase the computational speed, we ask for demanding computations to be parallelly performed on 8 cores:
Parallelism().set(nproc=8)
We start by declaring $\mathbb{S}^3$ as a differentiable manifold of dimension 3 over $\mathbb{R}$:
S3 = Manifold(3, 'S^3', latex_name=r'\mathbb{S}^3', start_index=1)
The first argument, 3
, is the dimension of the manifold, while the second argument is the symbol used to label the manifold, with the LaTeX output specified by the argument latex_name
. 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,3\}$; the default value is start_index=0
, yielding to $\{0,1,2\}$.
print(S3)
3-dimensional differentiable manifold S^3
S3
The 3-sphere cannot be covered by a single chart. At least two charts are necessary, for instance the charts associated with the stereographic projections from two distinct points, $N$ and $S$ say,
which we may call the North pole and the South pole respectively. Let us introduce the open subsets covered by these two charts:
$$ U := \mathbb{S}^3\setminus\{N\} $$
$$ V := \mathbb{S}^3\setminus\{S\} $$
U = S3.open_subset('U') ; print(U)
Open subset U of the 3-dimensional differentiable manifold S^3
V = S3.open_subset('V') ; print(V)
Open subset V of the 3-dimensional differentiable manifold S^3
We declare that $\mathbb{S}^3 = U \cup V$:
S3.declare_union(U, V)
Then we introduce the stereographic chart on $U$, denoting by $(x,y,z)$ the coordinates resulting from the stereographic projection from the North pole onto the equatorial plane:
stereoN.<x,y,z> = U.chart()
stereoN
stereoN.coord_range()
Similarly, we introduce on $V$ the coordinates $(x',y',z')$ corresponding to the stereographic projection from the South pole onto the equatorial plane:
stereoS.<xp,yp,zp> = V.chart("xp:x' yp:y' zp:z'")
stereoS
stereoS.coord_range()
We have to specify the transition map between the charts stereoN
= $(U,(x,y,z))$ and stereoS
= $(V,(x',y',z'))$; it is given by the standard inversion formulas:
r2 = x^2+y^2+z^2
stereoN_to_S = stereoN.transition_map(stereoS,
(x/r2, y/r2, z/r2),
intersection_name='W',
restrictions1= x^2+y^2+z^2!=0,
restrictions2= xp^2+yp^2+zp^2!=0)
stereoN_to_S.display()
In the above declaration, 'W'
is the name given to the open subset where the two charts overlap: $W := U\cap V$, the condition $x^2+y^2+z^2\not=0$ defines $W$ as a subset of $U$, and the condition $x'^2+y'^2+z'^2\not=0$ defines $W$ as a subset of $V$.
The inverse coordinate transformation is computed by means of the method inverse()
:
stereoS_to_N = stereoN_to_S.inverse()
stereoS_to_N.display()
Note that the situation is of course perfectly symmetric regarding the coordinates $(x,y,z)$ and $(x',y',z')$.
At this stage, the user's atlas has four charts:
S3.atlas()
For future reference, we store $W=U\cap V$ into a Python variable:
W = U.intersection(V)
print(W)
Open subset W of the 3-dimensional differentiable manifold S^3
$N$ is the point of $V$ of stereographic coordinates $(x',y',z')=(0,0,0)$:
N = V((0,0,0), chart=stereoS, name='N')
print(N)
Point N on the 3-dimensional differentiable manifold S^3
while $S$ is the point of $U$ of stereographic coordinates $(x,y,z)=(0,0,0)$:
S = U((0,0,0), chart=stereoN, name='S')
print(S)
Point S on the 3-dimensional differentiable manifold S^3
We have of course
all([N not in U, N in V, S in U, S not in V])
Let us first declare $\mathbb{R}^4$ as a 4-dimensional manifold covered by a single chart (the so-called Cartesian coordinates):
R4 = Manifold(4, 'R^4', r'\mathbb{R}^4')
X4.<T,X,Y,Z> = R4.chart()
X4
The embedding of $\mathbb{S}^3$ into $\mathbb{R}^4$ is then defined by the standard formulas relating the stereographic coordinates to the ambient Cartesian ones when considering a stereographic projection from the point $(-1,0,0,0)$ to the equatorial plane $T=0$:
rp2 = xp^2 + yp^2 + zp^2
Phi = S3.diff_map(R4, {(stereoN, X4):
[(1-r2)/(r2+1), 2*x/(r2+1),
2*y/(r2+1), 2*z/(r2+1)],
(stereoS, X4):
[(rp2-1)/(rp2+1), 2*xp/(rp2+1),
2*yp/(rp2+1), 2*zp/(rp2+1)]},
name='Phi', latex_name=r'\Phi')
Phi.display()
With this choice of stereographic projection, the "North" pole is actually the point of coordinates $(-1,0,0,0)$ in $\mathbb{R}^4$:
X4(Phi(N))
while the "South" pole is the point of coordinates $(1,0,0,0)$:
X4(Phi(S))
The hyperspherical coordinates $(\chi, \theta, \phi)$ generalize the standard spherical coordinates $(\theta, \phi)$ on $\mathbb{S}^2$. They are defined on the open domain $A\subset W \subset \mathbb{S}^3$ that is the complement of the "origin meridian"; since the latter is defined by $y=0$ and $x\geq 0$, we declare:
A = W.open_subset('A', coord_def={stereoN.restrict(W): (y!=0, x<0),
stereoS.restrict(W): (yp!=0, xp<0)})
print(A)
Open subset A of the 3-dimensional differentiable manifold S^3
We then declare the chart $(A,(\chi,\theta,\phi))$ by specifying the intervals spanned by the various coordinates:
spher.<ch,th,ph> = A.chart(r'ch:(0,pi):\chi th:(0,pi):\theta ph:(0,2*pi):\phi')
spher
spher.coord_range()
The specification of the hyperspherical coordinates is completed by providing the transition map to the stereographic chart $(A,(x,y,z))$:
den = 1 + cos(ch)
spher_to_stereoN = spher.transition_map(stereoN.restrict(A),
(sin(ch)*sin(th)*cos(ph)/den,
sin(ch)*sin(th)*sin(ph)/den,
sin(ch)*cos(th)/den))
spher_to_stereoN.display()
We also provide the inverse transition map, asking to check that the provided formulas are indeed correct (argument verbose=True
):
spher_to_stereoN.set_inverse(2*atan(sqrt(x^2+y^2+z^2)),
atan2(sqrt(x^2+y^2), z),
atan2(-y, -x)+pi,
verbose=True)
Check of the inverse coordinate transformation: ch == 2*arctan(sqrt(-cos(ch) + 1)/sqrt(cos(ch) + 1)) th == arctan2(sqrt(-cos(ch) + 1)*sin(th)/sqrt(cos(ch) + 1), cos(th)*sin(ch)/(cos(ch) + 1)) ph == pi - arctan2(sin(ch)*sin(ph)*sin(th)/(cos(ch) + 1), -cos(ph)*sin(ch)*sin(th)/(cos(ch) + 1)) x == x y == y z == z
The check is passed, modulo some lack of trigonometric simplifications in the first three lines.
spher_to_stereoN.inverse().display()
The transition map $(A,(\chi,\theta,\phi))\rightarrow (A,(x',y',z'))$ is obtained by combining the transition maps $(A,(\chi,\theta,\phi))\rightarrow (A,(x,y,z))$ and $(A,(x,y,z))\rightarrow (A,(x',y',z'))$:
spher_to_stereoS = stereoN_to_S.restrict(A) * spher_to_stereoN
spher_to_stereoS.display()
Similarly, the transition map $(A,(x',y',z'))\rightarrow (A,(\chi,\theta,\phi))$ is obtained by combining the transition maps $(A,(x',y',z'))\rightarrow (A,(x,y,z))$ and $(A,(x,y,z))\rightarrow (A,(\chi,\theta,\phi))$:
stereoS_to_spher = spher_to_stereoN.inverse() * stereoS_to_N.restrict(A)
stereoS_to_spher.display()
At this stage, the user atlas of $\mathbb{S}^3$ is
S3.atlas()
Let us get the coordinate expression of the restriction of the embedding $\Phi$ to $A$:
Phi.display(stereoN.restrict(A), X4)
Phi.display(spher, X4)
Phi.display()
We will need some projection operators from (a subset of) $\mathbb{R}^4$ to $\mathbb{S}^3$.
First, let $\mathbb{R}^4_N$ be $\mathbb{R}^4$ minus the hyperplane $T=-1$:
R4N = R4.open_subset('R4N', latex_name=r'\mathbb{R}^4_N',
coord_def={X4: T!=-1})
X4N = X4.restrict(R4N)
and let us consider the following projection $\Pi_N: \mathbb{R}^4_N \to U\subset\mathbb{S}^3$:
ProjN = R4N.diff_map(U, {(X4N, stereoN):
[X/(1+T), Y/(1+T), Z/(1+T)]},
name='P_N', latex_name=r'\Pi_N')
ProjN.display()
Similarly, let $\mathbb{R}^4_S$ be $\mathbb{R}^4$ minus the hyperplane $T=1$ and $\Pi_S$ the following projection $\mathbb{R}_S\to V\subset \mathbb{S}^3$:
R4S = R4.open_subset('R4S', latex_name=r'\mathbb{R}^4_S',
coord_def={X4: T!=1})
X4S = X4.restrict(R4S)
ProjS = R4S.diff_map(V, {(X4S, stereoS):
[X/(1-T), Y/(1-T), Z/(1-T)]},
name='P_S', latex_name=r'\Pi_S')
ProjS.display()
Let us check that once applied to an embedded point of $U\cap V\subset \mathbb{S}^3$, this projection reduces to the identity:
var('a b c', domain='real')
p = S3((1+a^2,b,c), chart=stereoN)
stereoN(p)
all([p in U, p in V])
all([ProjN(Phi(p)) == p, ProjS(Phi(p)) == p])
p = S3((1+a^2,b,c), chart=stereoS)
all([ProjN(Phi(p)) == p, ProjS(Phi(p)) == p])
q = R4((sqrt(3)/2, sin(a)*cos(b)/2, sin(a)*sin(b)/2, cos(a)/2))
X4(q)
all([q in R4N, q in R4S])
all([Phi(ProjN(q)) == q, Phi(ProjS(q)) == q])
We consider the (division) algebra of quaternions $\mathbb{H}$ as $\mathbb{R}^4$ endowed with the following (non-commutative) product:
def qprod(p,q):
if p in R4 and q in R4:
T1, X1, Y1, Z1 = X4(p)
T2, X2, Y2, Z2 = X4(q)
return R4(((T1*T2-X1*X2-Y1*Y2-Z1*Z2).simplify_full(),
(T1*X2+X1*T2+Y1*Z2-Z1*Y2).simplify_full(),
(T1*Y2-X1*Z2+Y1*T2+Z1*X2).simplify_full(),
(T1*Z2+X1*Y2-Y1*X2+Z1*T2).simplify_full()))
if p in S3 and q in S3:
a = qprod(Phi(p),Phi(q))
if X4(a) == (-1,0,0,0):
return N
return ProjN(R4N(a))
raise ValueError("Cannot evaluate qprod of {} and {}".format(p,q))
Note that we have extended the definition of the quaternionic product to $\mathbb{S}^3$ via the embedding $\Phi$.
Let us introduce two special points on $\mathbb{S}^3$: $\mathbf{1}$ and $-\mathbf{1}$.
One = S3((0,0,0), chart=stereoN, name='1', latex_name=r'\mathbf{1}')
X4(Phi(One))
As we can see from the Cartesian coordinates of $\Phi(\mathbf{1})$, the point $\mathbf{1}$ is actually nothing but the "South" pole used to define the stereographic chart $(V,(x',y',z'))$:
One == S
minusOne = S3((0,0,0), chart=stereoS, name='-1', latex_name=r'-\mathbf{1}')
X4(Phi(minusOne))
The point $\mathbf{-1}$ is thus nothing but the "North" pole used to define the stereographic chart $(U,(x,y,z))$:
minusOne == N
Next we introduce the points $\mathbf{i}$, $\mathbf{j}$ and $\mathbf{k}$ on $\mathbb{S}^3$:
I = S3((1,0,0), chart=stereoN, name='i', latex_name=r'\mathbf{i}')
X4(Phi(I))
stereoS(I)
J = S3((0,1,0), chart=stereoN, name='j', latex_name=r'\mathbf{j}')
X4(Phi(J))
stereoS(J)
Since $\mathbf{j}$ lies in $A$, contrary to $\mathbf{i}$, we may ask for its hyperspherical coordinates:
spher(J)
K = S3((0,0,1), chart=stereoN, name='k', latex_name=r'\mathbf{k}')
X4(Phi(K))
stereoS(K)
Hamilton's fundamental relations $$ \mathbf{i} \mathbf{j} \mathbf{k} = \mathbf{-1} $$ $$ \mathbf{i} \mathbf{j} = \mathbf{k},\quad \mathbf{j} \mathbf{k} = \mathbf{i}, \quad \mathbf{k} \mathbf{i} = \mathbf{j}$$ are satisfied:
qprod(I, qprod(J,K)) == minusOne
all([qprod(I,J) == K, qprod(J,K) == I,
qprod(K,I) == J])
These relations imply $\mathbf{i}^2 = \mathbf{-1}$, $\mathbf{j}^2 = \mathbf{-1}$ and $\mathbf{k}^2 = \mathbf{-1}$:
all([qprod(One,One) == One, qprod(I,I) == minusOne,
qprod(J,J) == minusOne, qprod(K,K) == minusOne])
Let us introduce $\mathbf{-i}$, $\mathbf{-j}$ and $\mathbf{-k}$, as points of $\mathbb{S}^3$:
minusI = qprod(minusOne, I)
X4(Phi(minusI))
minusJ = qprod(minusOne, J)
X4(Phi(minusJ))
minusK = qprod(minusOne, K)
X4(Phi(minusK))
In the comments below (but not in the SageMath code), we shall identify $\mathbf{1}\in \mathbb{S}^3$ with $\Phi(\mathbf{1})\in \mathbb{R}^4$, $\mathbf{i}\in \mathbb{S}^3$ with $\Phi(\mathbf{i})\in \mathbb{R}^4$, etc. In particular, we consider $(\mathbf{1}, \mathbf{i}, \mathbf{j},\mathbf{k})$ as a basis of the quaternion algebra $\mathbb{H}$.
The conjugate of a quaternion $q = T + X\mathbf{i} + Y\mathbf{j} + Z\mathbf{k}$ is $\bar{q} = T - X\mathbf{i} - Y\mathbf{j} - Z\mathbf{k}$; hence we define:
def qconj(p):
if p in R4:
T, X, Y, Z = X4(p)
return R4((T, -X, -Y, -Z))
if p in S3:
a = qconj(Phi(p))
if X4(a) == (-1,0,0,0):
return N
return ProjN(a)
raise ValueError("Cannot evaluate qconf of {}".format(p))
In particular, we have $\bar{\mathbf{1}} = \mathbf{1}$, $\bar{\mathbf{i}} = -\mathbf{i}$, $\bar{\mathbf{j}} = -\mathbf{j}$ and $\bar{\mathbf{k}} = -\mathbf{k}$:
all([qconj(One) == One,
qconj(I) == minusI,
qconj(J) == minusJ,
qconj(K) == minusK])
The conjugate of an element of $\mathbb{S}^3$
assume(a != 0) # to ensure that qconj(p) is not N
p = S3((a,b,c), chart=stereoN)
stereoN(qconj(p))
p = S3((a,b,c), chart=stereoS)
stereoS(qconj(p))
forget(a!=0)
The quaternionic norm $\| q\| = \sqrt{q\bar{q}}$ coincide with the Euclidean norm in $\mathbb{R}^4$, so that $\mathbb{S}^3$ can be viewed as the set of unit quaternions; hence we define:
def qnorm(p):
if p in R4:
T, X, Y, Z = X4(p)
return (sqrt(T^2 + X^2 + Y^2 + Z^2)).simplify_full()
if p in S3:
return 1
raise ValueError("Cannot evaluate qnorm of {}".format(p))
var('d', domain='real')
q = R4((a,b,c,d))
qnorm(q)
Let us check that $\| q\|^2 = q\bar{q}$:
R4((qnorm(q)^2,0,0,0)) == qprod(q, qconj(q))
As elements of $\mathbb{S}^3$, $\mathbf{1}$, $\mathbf{i}$, $\mathbf{j}$ and $\mathbf{k}$ have all unit norm:
(qnorm(One), qnorm(I), qnorm(J), qnorm(K)) == (1, 1, 1, 1)
The right translation by $\mathbf{i}$ is the map $\bar{R}_{\mathbf{i}}: p \mapsto p \mathbf{i}$. We define it first at the level of $\mathbb{R}^4$:
p = R4((T,X,Y,Z))
RI_R4 = R4.diff_map(R4, X4(qprod(p, Phi(I))))
RI_R4.display()
Focusing on its action on $\mathbb{S}^3$, we consider then the map ${\bar R}_{\mathbf{i}}\circ\Phi$:
RI_S3_R4 = RI_R4 * Phi
RI_S3_R4.display()
Let $U_{\mathbf{i}} := U \setminus \{\mathbf{i}\}$; since the coordinates of $\mathbf{i}$ in the chart $(U,(x,y,z))$ are $(1,0,0)$, we declare $U_{\mathbf{i}}$ as
UI = U.open_subset('U_I', latex_name=r'U_{\mathbf{i}}',
coord_def={stereoN: (x!=1, y!=0, z!=0)})
If we restrict $R_{\mathbf{i}}\circ\Phi$ to $U_{\mathbf{i}}$ the codomain can be taken to be $\mathbb{R}^4_N$ since $\mathbf{i}$ is the only point of $\mathbb{S}^3$ for which $T(R_{\mathbf{i}}(p)) = -1$. Hence we may apply the operator $\Pi_N$ to define the right translation by $\mathbf{i}$ as a map $U_{\mathbf{i}}\to U$:
RI_UI = ProjN * RI_S3_R4.restrict(UI, subcodomain=R4N)
RI_UI.display()
Similarly, if we restrict $R_{\mathbf{i}}\circ\Phi$ to $V_{-\mathbf{i}} := V \setminus\{-\mathbf{i}\}$, we get a map $V_{-\mathbf{i}} \to \mathbb{R}^4_{S}$, so that composing by $\Pi_S$, the right translation by $\mathbf{i}$ becomes a map $V_{-\mathbf{i}}\to V$:
VmI = V.open_subset('V_mI', latex_name=r'V_{-\mathbf{i}}',
coord_def={stereoS: (xp!=-1, yp!=0, zp!=0)})
RI_VmI = ProjS * RI_S3_R4.restrict(VmI, subcodomain=R4S)
RI_VmI.display()
We note that $\mathbb{S}^3 = U_{\mathbf{i}} \cup V_{-\mathbf{i}}$:
S3.declare_union(UI, VmI)
Consequently, we can define the right translation by $\mathbf{i}$ as a map $R_{\mathbf{i}}: \mathbb{S}^3\to \mathbb{S}^3$ by providing the coordinate expressions obtained above on $U_{\mathbf{i}}$ and $V_{-\mathbf{i}}$:
RI = S3.diff_map(S3, name='R_I', latex_name=r'{R_{\mathbf{i}}}')
RI.add_expression(stereoN.restrict(UI), stereoN,
RI_UI.expr(stereoN.restrict(UI), stereoN))
RI.add_expression(stereoS.restrict(VmI), stereoS,
RI_VmI.expr(stereoS.restrict(VmI), stereoS))
RI.display(stereoN.restrict(UI), stereoN)
RI.display(stereoS.restrict(VmI), stereoS)
Let us check the formulas $R_{\mathbf{i}}(\mathbf{1})=\mathbf{i}$, $R_{\mathbf{i}}(-\mathbf{1})=-\mathbf{i}$, $R_{\mathbf{i}}(\mathbf{i})=-\mathbf{1}$, $R_{\mathbf{i}}(-\mathbf{i})=\mathbf{1}$, $R_{\mathbf{i}}(\mathbf{j})=-\mathbf{k}$, $R_{\mathbf{i}}(-\mathbf{j})=\mathbf{k}$, $R_{\mathbf{i}}(\mathbf{k})=\mathbf{j}$ and $R_{\mathbf{i}}(-\mathbf{k})=-\mathbf{j}$:
all([RI(One)==I, RI(minusOne)==minusI,
RI(I)==minusOne, RI(minusI)==One,
RI(J)==minusK, RI(minusJ)==K,
RI(K)==J, RI(minusK)==minusJ])
Let us recall the expression of the right translation by $\mathbf{i}$ on $\mathbb{R}^4$:
RI_R4.display()
We turn it into a vector field $E_{\mathbf{i}}$ on $\mathbb{R}^4$ by identifying $T_p\mathbb{R}^4$ and $\mathbb{R}^4$ at each point $p\in\mathbb{R}^4$:
EI = R4.vector_field(name='E_I', latex_name=r'E_{\mathbf{i}}')
EI[:] = RI_R4.expression()
EI.display()
The "radial" vector field on $\mathbb{R}^4$ is
r = R4.vector_field(name='r')
r[:] = (T,X,Y,Z)
r.display()
It is clear that $r\cdot E_{\mathbf{i}}=0$, where $\cdot$ denotes the standard Euclidean scalar product in $\mathbb{R}^4$. We can check this property explicitely by introducing the Euclidean metric:
h = R4.metric('h')
h[0,0], h[1,1], h[2,2], h[3, 3] = 1, 1, 1, 1
h.display()
so that $r\cdot E_{\mathbf{i}} = h(r, E_{\mathbf{i}})$:
h(r, EI) == 0
This proves that the vector field $E_{\mathbf{i}}$ is tangent to $\mathbb{S}^3$, or more precisely to the embedded submanifold $\Phi(\mathbb{S}^3)$. Consequently, there exists a vector field $\varepsilon_{\mathbf{i}}$ on $\mathbb{S}^3$, the pushforward of which by $\Phi$ is $E_{\mathbf{i}}$: $$ E_{\mathbf{i}} = \Phi^* \varepsilon_{\mathbf{i}} $$ Let us determine the components of $\varepsilon_{\mathbf{i}}$ in the vector frame $\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right)$ associated with the stereographic coordinates on $U$:
frameN = stereoN.frame()
frameN
The pushforwards by $\Phi$ of the stereographic frame vectors are:
frameN_R4 = [Phi.pushforward(frameN[i]) for i in S3.irange()]
frameN_R4
print(frameN_R4[0])
Vector field Phi^*(d/dx) along the Open subset U of the 3-dimensional differentiable manifold S^3 with values on the 4-dimensional differentiable manifold R^4
The expressions of $\Phi^* \frac{\partial}{\partial x}$, $\Phi^* \frac{\partial}{\partial y}$ and $\Phi^* \frac{\partial}{\partial z}$ in terms of the canonical vector frame $\left(\frac{\partial}{\partial T}, \frac{\partial}{\partial X},\frac{\partial}{\partial Y},\frac{\partial}{\partial Z}\right)$ of $\mathbb{R}^4$ are
frameN_R4[0].display()
frameN_R4[1].display()
frameN_R4[2].display()
Let us denote by $(a,b,c)$ the components of $\left. \varepsilon_{\mathbf{i}} \right|_{U}$ in the stereographic frame: $$ \left. \varepsilon_{\mathbf{i}} \right|_{U} = a \frac{\partial}{\partial x} + b \frac{\partial}{\partial y} + c \frac{\partial}{\partial z} $$ Since $\Phi^*\varepsilon_{\mathbf{i}}=\left. E_{\mathbf{i}} \right| _{\Phi(\mathbb{S}^3)}$, we get the linear system $$ a\, \Phi^* \frac{\partial}{\partial x} + b \, \Phi^* \frac{\partial}{\partial y} + c\, \Phi^* \frac{\partial}{\partial z} = \left. E_{\mathbf{i}} \right| _{\Phi(U)} $$ to be solved in $(a,b,c)$. The right-hand side is
p = U((x,y,z), chart=stereoN, name='p')
EIp = EI.at(Phi(p))
EIp[:]
Hence the system:
eqs = [(a*frameN_R4[0][i] + b*frameN_R4[1][i] + c*frameN_R4[2][i]).expr() == EIp[i]
for i in R4.irange()]
eqs
The unique solution is
sol = solve(eqs, (a,b,c), solution_dict=True)
sol
The expression of $\varepsilon_{\mathbf{i}}$ in terms of the frame associated with the stereographic coordinates on $U$ is thus
epsI = sol[0][a] * frameN[1] + sol[0][b] * frameN[2] + sol[0][c] * frameN[3]
epsI.display()
The vector field $E_{\mathbf{i}}$ never vanishes on $\Phi(\mathbb{S}^3)$, since it vanishes only at $(T,X,Y,Z)=(0,0,0,0)$. It follows that $\varepsilon_{\mathbf{i}}$ never vanishes on $\mathbb{S}^3$. Similarly, starting from the right translation by $\mathbf{j}$ and the right translation by $\mathbf{k}$, we can construct global vector fields $\varepsilon_{\mathbf{j}}$ and $\varepsilon_{\mathbf{k}}$ that are always nonzero on $\mathbb{S}^3$. Moreover, the vector fields $\varepsilon_{\mathbf{i}}$, $\varepsilon_{\mathbf{j}}$ and $\varepsilon_{\mathbf{k}}$ are linearly independent at any point of $\mathbb{S}^3$. They thus form a global vector frame of $\mathbb{S}^3$. This means that, as any Lie group, $\mathbb{S}^3$ is a parallelizable manifold. This contrasts with $\mathbb{S}^2$. Actually the only parallelizable spheres are $\mathbb{S}^1$, $\mathbb{S}^3$ and $\mathbb{S}^7$.
Let us declare the global vector frame $(\varepsilon_{\mathbf{i}}, \varepsilon_{\mathbf{j}}, \varepsilon_{\mathbf{k}})$, using the notations $\varepsilon_1 := \varepsilon_{\mathbf{i}}$, $\varepsilon_2 := \varepsilon_{\mathbf{j}}$ and $\varepsilon_3 := \varepsilon_{\mathbf{k}}$:
E = S3.vector_frame('E', latex_symbol=r'\varepsilon')
E
On $U$, we can set the components of $\varepsilon_1$ in the stereographic frame $\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right)$ to those obtained above:
E[1].restrict(U)[stereoN.frame(),:,stereoN] = (sol[0][a], sol[0][b], sol[0][c])
E[1].display(stereoN.frame())
Let us check that the pushforward of $\varepsilon_1$ by $\Phi$ coincides with $E_{\mathbf{i}}$:
E1U_R4 = Phi.pushforward(E[1].restrict(U))
E1U_R4.display()
all([E1U_R4[i] == EIp[i] for i in R4.irange()])
Let us now determine the components of $\varepsilon_1$ in the vector frame $\left(\frac{\partial}{\partial x'}, \frac{\partial}{\partial y'}, \frac{\partial}{\partial z'}\right)$ associated with the stereographic chart from the South pole:
frameS = stereoS.frame()
frameS
We use the same procedure as for the stereographic frame from the North pole:
frameS_R4 = [Phi.pushforward(frameS[i]) for i in S3.irange()]
frameS_R4
p = V((xp,yp,zp), chart=stereoS, name='p')
EIp = EI.at(Phi(p))
EIp[:]
eqs = [(a*frameS_R4[0][i] + b*frameS_R4[1][i] + c*frameS_R4[2][i]).expr() == EIp[i]
for i in R4.irange()]
eqs
sol = solve(eqs, (a,b,c), solution_dict=True)
sol
E[1].restrict(V)[stereoS.frame(),:, stereoS] = (sol[0][a], sol[0][b], sol[0][c])
E[1].display(stereoS.frame())
Again, we check the correctness by
E1V_R4 = Phi.pushforward(E[1].restrict(V))
all([E1V_R4[i] == EIp[i] for i in R4.irange()])
The vector field $\varepsilon_2 = \varepsilon_{\mathbf{j}}$ is induced by the right translation by $\mathbf{j}$:
p = R4((T,X,Y,Z))
RJ_R4 = R4.diff_map(R4, X4(qprod(p, Phi(J))))
RJ_R4.display()
EJ = R4.vector_field(name='E_J', latex_name=r'E_{\mathbf{j}}')
EJ[:] = RJ_R4.expression()
EJ.display()
We determine the components of $\varepsilon_2$ in the stereographic frame from the North pole as we did for $\varepsilon_1$:
p = U((x,y,z), chart=stereoN)
EJp = EJ.at(Phi(p))
EJp[:]
eqs = [(a*frameN_R4[0][i] + b*frameN_R4[1][i] + c*frameN_R4[2][i]).expr() == EJp[i]
for i in R4.irange()]
sol = solve(eqs, (a,b,c), solution_dict=True)
E[2].restrict(U)[stereoN.frame(),:,stereoN] = (sol[0][a], sol[0][b], sol[0][c])
E[2].display(stereoN.frame())
Check:
E2U_R4 = Phi.pushforward(E[2].restrict(U))
all([E2U_R4[i] == EJp[i] for i in R4.irange()])
We turn now to the stereographic components from the South pole:
p = V((xp,yp,zp), chart=stereoS)
EJp = EJ.at(Phi(p))
EJp[:]
eqs = [(a*frameS_R4[0][i] + b*frameS_R4[1][i] + c*frameS_R4[2][i]).expr() == EJp[i]
for i in R4.irange()]
sol = solve(eqs, (a,b,c), solution_dict=True)
E[2].restrict(V)[stereoS.frame(),:, stereoS] = (sol[0][a], sol[0][b], sol[0][c])
E[2].display(stereoS.frame())
E2V_R4 = Phi.pushforward(E[2].restrict(V))
all([E2V_R4[i] == EJp[i] for i in R4.irange()])
The vector field $\varepsilon_3 = \varepsilon_{\mathbf{k}}$ is induced by the right translation by $\mathbf{k}$:
p = R4((T,X,Y,Z))
RK_R4 = R4.diff_map(R4, X4(qprod(p, Phi(K))))
RK_R4.display()
EK = R4.vector_field(name='E_K', latex_name=r'E_{\mathbf{k}}')
EK[:] = RK_R4.expression()
EK.display()
The components of $\varepsilon_3$ in the two stereographic frames are obtained in the same manner as for $\varepsilon_1$ and $\varepsilon_2$
p = U((x,y,z), chart=stereoN)
EKp = EK.at(Phi(p))
eqs = [(a*frameN_R4[0][i] + b*frameN_R4[1][i] + c*frameN_R4[2][i]).expr() == EKp[i]
for i in R4.irange()]
sol = solve(eqs, (a,b,c), solution_dict=True)
E[3].restrict(U)[stereoN.frame(),:,stereoN] = (sol[0][a], sol[0][b], sol[0][c])
E[3].display(stereoN.frame())
E3U_R4 = Phi.pushforward(E[3].restrict(U))
all([E3U_R4[i] == EKp[i] for i in R4.irange()])
p = V((xp,yp,zp), chart=stereoS)
EKp = EK.at(Phi(p))
eqs = [(a*frameS_R4[0][i] + b*frameS_R4[1][i] + c*frameS_R4[2][i]).expr() == EKp[i]
for i in R4.irange()]
sol = solve(eqs, (a,b,c), solution_dict=True)
E[3].restrict(V)[stereoS.frame(),:, stereoS] = (sol[0][a], sol[0][b], sol[0][c])
E[3].display(stereoS.frame())
E3V_R4 = Phi.pushforward(E[3].restrict(V))
all([E3V_R4[i] == EKp[i] for i in R4.irange()])
The vector fields on $\mathbb{R}^4$ induced by the right translations by respectively $\mathbf{i}$, $\mathbf{j}$ and $\mathbf{k}$ are
show(EI.display())
show(EJ.display())
show(EK.display())
As a check, we note that the above formulas coincide with those given in Problem 8-6 of Lee's textbook Introduction to Smooth Manifolds, 2nd ed., Springer (New York) (2013).
The vector fields $E_{\mathbf{i}}$, $E_{\mathbf{j}}$ and $E_{\mathbf{k}}$ are tangent to $\Phi(\mathbb{S}^3)\subset\mathbb{R}^4$ and therefore induce three vectors fields on $\mathbb{S}^3$, $\varepsilon_1$, $\varepsilon_2$ and $\varepsilon_3$, which form a global vector frame of $\mathbb{S}^3$. The components of these vector fields with respect to the (local) stereographic frames are
for i in S3.irange():
show(E[i].display(stereoN.frame()))
print(" ")
for i in S3.irange():
show(E[i].display(stereoS.frame()))
To complete the connection between the global frame $(\varepsilon_i)$ and the stereographic frames in $U$ and $V$, we shall declare the change-of-basis formulas related the two frames, so that they can be used automatically by SageMath when necessary. Starting with $U$, we first consider the restriction of the global frame to $U$:
E_U = E.restrict(U); E_U
The automorphism $P$ such that on $U$, $\varepsilon_i = P\left(\frac{\partial}{\partial x^i}\right)$ is read directely on the components of $\varepsilon_i$ with respect to the frame $\left(\frac{\partial}{\partial x^i}\right) = \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} \right)$:
P = U.automorphism_field()
for i in S3.irange():
for j in S3.irange():
P[j,i] = E_U[i][j]
P[:]
Check of the formula $\varepsilon_i = P\left(\frac{\partial}{\partial x^i}\right)$:
all([E_U[i] == P(frameN[i]) for i in S3.irange()])
We declare the change of frame by the method set_change_of_frame
:
U.set_change_of_frame(frameN, E_U, P)
The inverse is automatically computed:
U.change_of_frame(E_U, frameN)[:]
We do the same thing on $V$:
E_V = E.restrict(V); E_V
P = V.automorphism_field()
for i in S3.irange():
for j in S3.irange():
P[j,i] = E_V[i][j]
P[:]
all([E_V[i] == P(frameS[i]) for i in S3.irange()])
V.set_change_of_frame(frameS, E_V, P)
V.change_of_frame(E_V, frameS)[:]
Because it admits a global vector frame, $\mathbb{S}^3$ is a parallelizable manifold:
S3.is_manifestly_parallelizable()
Equivalently, the set $\mathfrak{X}(\mathbb{S}^3)$ of vector fields on $\mathbb{S}^3$ is a free module of finite rank over the algebra $C^\infty(\mathbb{S}^3)$ (the algebra of smooth scalar fields on $\mathbb{S}^3$:
XS3 = S3.vector_field_module()
XS3
isinstance(XS3, FiniteRankFreeModule)
print(XS3.category())
Category of finite dimensional modules over Algebra of differentiable scalar fields on the 3-dimensional differentiable manifold S^3
As a free module, $\mathfrak{X}(\mathbb{S}^3)$ admits a basis, which is nothing but the global vector frame constructed above:
XS3.default_basis()
On $U$, we may compute the Lie brackets $[\varepsilon_i, \varepsilon_j]$ of two vectors of the global frame:
E_U[1].bracket(E_U[2]).display(E_U)
E_U[1].bracket(E_U[3]).display(E_U)
E_U[2].bracket(E_U[3]).display(E_U)
Equivalently, the structure coefficients of the frame $\varepsilon$ are
C = E_U.structure_coeff(); C
C.display('C')
By definition, the structure coefficients $C_{kij}$ obey the relation $[\varepsilon_i, \varepsilon_j] = C_{kij} \, \varepsilon_k$, as we can check:
E_U[1].bracket(E_U[2]) == sum(C[[k,1,2]]*E_U[k] for k in S3.irange())
The Lie algebra of $\mathbb{S}^3$ can be identified with the tangent space at the unit element $\mathbf{1}$:
T1 = S3.tangent_space(One)
T1
The values taken by the global frame vectors at $\mathbf{1}$ are
for i in S3.irange():
show(E[i].at(One).display())
Each of these vectors belongs to $T_{\mathbf{1}}\mathbb{S}^3$:
E[1].at(One).parent()
all([E[i].at(One).parent() is T1 for i in S3.irange()])
The simple expressions of $(\varepsilon_1, \varepsilon_2, \varepsilon_3)$ in terms of the basis $\left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\right)$ of $T_{\mathbf{1}}\mathbb{S}^3$ induced by the stereographic coordinates stems from the fact that the $\mathbb{R}^4$ vectors $(E_{\mathbf{i}}, E_{\mathbf{j}}, E_{\mathbf{k}})$ coincide with the vectors $\left(\frac{\partial}{\partial X}, \frac{\partial}{\partial Y}, \frac{\partial}{\partial Z}\right)$ at $\Phi(\mathbf{1})$:
EI.at(Phi(One)).display()
EJ.at(Phi(One)).display()
EK.at(Phi(One)).display()
The extra factor $1/2$, which appears in the above expressions of $(\varepsilon_1, \varepsilon_2, \varepsilon_3)$, arises from the fact that the stereographic coordinates $(x,y,z)$ have been defined with respect to the hyperplane $T=0$ and not to the hyperplane $T=1$ (to which $\Phi(\mathbf{1})$ belongs).
The Hopf coordinates have been introduced in the worksheet 3-sphere: charts, quaternions and the Hopf fribration.
B = U.open_subset('B', coord_def={stereoN.restrict(U):
[x^2+y^2!=0, x^2+y^2+z^2!=1,
(1-x^2-y^2-z^2)*x-2*y*z!=0]})
Hcoord.<eta,alp,bet> = B.chart(r"eta:(0,pi/2):\eta alpha:(0,2*pi):\alpha beta:(0,2*pi):\beta")
Hcoord
Hcoord_to_stereoN = Hcoord.transition_map(
stereoN.restrict(U),
(sin(eta)*cos(alp+bet)/(1+cos(eta)*sin(alp)),
sin(eta)*sin(alp+bet)/(1+cos(eta)*sin(alp)),
cos(eta)*cos(alp)/(1+cos(eta)*sin(alp))))
Hcoord_to_stereoN.display()
Hcoord_to_stereoN.set_inverse(asin(2*sqrt(x^2+y^2)/(1+x^2+y^2+z^2)),
atan2(x^2+y^2+z^2-1, -2*z) + pi,
atan2(-y,-x) - atan2(x^2+y^2+z^2-1, -2*z))
Hcoord_to_stereoN.inverse().display()
Let us express the vectors of the global frame $\varepsilon$ in terms of the vector frame $\left(\frac{\partial}{\partial\eta}, \frac{\partial}{\partial\alpha}, \frac{\partial}{\partial\beta}\right)$ associated with the Hopf coordinates:
E_U[1].display(Hcoord.frame(), Hcoord)
We notice that the components can be simplified to $\varepsilon_1^{\ \, \eta}=\sin(2\alpha+\beta)$, $\varepsilon_1^{\ \, \alpha}=-\cos(2\alpha+\beta)\tan\eta$ and
$\varepsilon_1^{\ \, \beta}=\cos(2\alpha+\beta)/(\cos\eta\sin\eta)$. Hence, we substitute these values via the method add_comp
(NB: the method set_comp
would erase the other sets of components):
E_U[1].add_comp(Hcoord.frame())[1, Hcoord] = sin(2*alp+bet)
E_U[1].add_comp(Hcoord.frame())[2, Hcoord] = -cos(2*alp+bet) * tan(eta)
E_U[1].add_comp(Hcoord.frame())[3, Hcoord] = cos(2*alp+bet) / (cos(eta)*sin(eta))
E_U[1].display(Hcoord.frame(), Hcoord)
Similarly, the components of $\varepsilon_2$ require some simplification "by hand":
E_U[2].display(Hcoord.frame(), Hcoord)
E_U[2].add_comp(Hcoord.frame())[1, Hcoord] = -cos(2*alp+bet)
E_U[2].add_comp(Hcoord.frame())[2, Hcoord] = -sin(2*alp+bet) * tan(eta)
E_U[2].add_comp(Hcoord.frame())[3, Hcoord] = sin(2*alp+bet) / (cos(eta)*sin(eta))
E_U[2].display(Hcoord.frame(), Hcoord)
The expression of $\varepsilon_3$ in terms of the vector frame $\left(\frac{\partial}{\partial\eta}, \frac{\partial}{\partial\alpha}, \frac{\partial}{\partial\beta}\right)$ is particularly simple:
E_U[3].display(Hcoord.frame(), Hcoord)
This is not surprising since $\varepsilon_3$ has been defined as the vector field induced by the right translation by $\mathbf{k}$ and the Hopf fibration considered the worksheet 3-sphere: charts, quaternions and the Hopf fribration has been defined with respect to $\mathbf{k}$ as
$$
\begin{array}{cccc}
h:& \mathbb{S}^3 & \to & \mathbb{S}^2\\
& q & \mapsto & q \mathbf{k} \bar{q}
\end{array}
$$
More precisely, we have shown in the above worksheet that $\frac{\partial}{\partial\alpha}$ is tangent to the $\mathbb{S}^1$ fibers of $h$ and it is easy to see that the field lines of $\varepsilon_3$ coincide with these fibers: let $q\in\mathbb{S}^3$ and $q'\in\mathbb{S}^3$ be a point infinitely close to $q$ on a field line of $\varepsilon_3$:
$$
q' = q + \lambda \left. \varepsilon_3 \right| _q = q + \lambda q \, \mathbf{k},
$$
where $\lambda = o(1)$ is an infinitesimal real number and the second equality stems from the very definition of $\varepsilon_3 = \varepsilon_{\mathbf{k}}$ as the vector field associated with the right translation by $\mathbf{k}$. We have then
$$
\overline{q'} = \overline{q + \lambda q \, \mathbf{k}}
= \bar{q} + \lambda \bar{\mathbf{k}} \bar{q}
= \bar{q} - \lambda \mathbf{k} \bar{q}
$$
so that
$$
h(q') = (q + \lambda q \, \mathbf{k})\mathbf{k}(\bar{q} - \lambda \mathbf{k} \bar{q}) .
$$
Expanding and using $\mathbf{k}^2 = -1$ and $q\bar{q} = 1$, we find
$h(q') = h(q)$ at first order in $\lambda$, which shows that $q$ and $q'$ belong to the same fiber of $h$.
Accordingly the field lines of $\varepsilon_3$ are circles, as we can check on a (projection) plot in the plane $z=0$:
graph_z0 = E[3].plot(chart=stereoN, ambient_coords=(x,y),
fixed_coords={z: 0}, max_range=1, scale=0.5)
show(graph_z0, aspect_ratio=1)
or in the plane $y=0$:
graph_y0 = E[3].plot(chart=stereoN, ambient_coords=(x,z),
fixed_coords={y: 0}, number_values=11,
ranges={x: (-2,2), z: (-1,1)}, scale=0.4)
show(graph_y0, aspect_ratio=1)
A 3D view of $\varepsilon_3$, in terms of the stereographic coordinates $(x,y,z)$:
graph = E[3].plot(chart=stereoN, max_range=1, number_values=7,
scale=0.25, label_axes=False)
show(graph, viewer=viewer3D, online=True, axes_labels=['x','y','z'])