x[ϕ_,ψ_]:=Cos[ϕ]Sin[ψ];
y[ϕ_,ψ_]:=-Cos[ϕ]Cos[ψ];
z[ϕ_,ψ_]:=Sin[ϕ];
nothing =Graphics3D[{Opacity->0,Point[{0,0,0}]}];

PlotDiagram[vertex_List,edge_List,args___]:=Module[{v,vl,vp,e,el,es},v=Range[Length@vertex];
vl=Table[vertex[[i,1]],{i,1,Length@vertex}];
vp=Table[vertex[[i,2]],{i,1,Length@vertex}];
e=Table[edge[[i,1]],{i,1,Length@edge}];
el=Table[edge[[i,2]],{i,1,Length@edge}];
es=Table[edge[[i,3]],{i,1,Length@edge}];
Graph[e,VertexCoordinates->vp,VertexLabels->Table[v[[i]]->Placed[vl[[i]],{0.5,0.5}],{i,1,Length[v]}],EdgeLabels->Table[e[[i]]->el[[i]],{i,1,Length[e]}],EdgeShapeFunction->Table[e[[i]]->es[[i]],{i,1,Length[e]}],args,VertexLabelStyle->Directive[Italic,18],EdgeLabelStyle->Directive[Italic,15],VertexSize->0,VertexStyle->Directive[EdgeForm[],White]]];

s1[al_,gam0_,ϕ_,θ_,t1_,cc0_]:=NDSolve[{q0'[t]==(q1[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q2[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q1'[t]==(-q0[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q3[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q2'[t]==(-q3[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q0[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q3'[t]==(q2[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q1[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,gam'[t]== cc[t],cc'[t]==-Sin[al]Sin[2 gam[t]],q0[0]==1,q1[0]==0,q2[0]==0,q3[0]==0,gam[0]==gam0,cc[0]==cc0
},{q0,q1,q2,q3,gam,cc},{t,t1}];

s2[al_,gam0_,ϕ_,θ_,t1_]:=NDSolve[{q0'[t]==(q1[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q2[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q1'[t]==(-q0[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q3[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q2'[t]==(-q3[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q0[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q3'[t]==(q2[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q1[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,gam'[t]== cc[t],cc'[t]==-Sin[al]Sin[2 gam[t]],q0[0]==1,q1[0]==0,q2[0]==0,q3[0]==0,gam[0]==gam0,cc[0]==x[θ,ϕ]*Sqrt[Cos[al]+Sin[al]]Sin[gam0]/z[θ,ϕ]-y[θ,ϕ]*Sqrt[Cos[al]-Sin[al]]Cos[gam0]/z[θ,ϕ]
},{q0,q1,q2,q3,gam,cc},{t,t1}];

s3[al_,gam0_,ϕ_,θ_,t1_,s_]:=NDSolve[{q0'[t]==(q1[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q2[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q1'[t]==(-q0[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]-q3[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q2'[t]==(-q3[t]Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q0[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,
q3'[t]==(q2[t] Sin[gam[t]] Sqrt[Cos[al]-Sin[al]]+q1[t] Cos[gam[t]] Sqrt[Cos[al]+Sin[al]])/2,gam'[t]== cc[t],cc'[t]==-Sin[al]Sin[2 gam[t]],q0[0]==1,q1[0]==0,q2[0]==0,q3[0]==0,gam[0]==(2s+1)Pi/2,cc[0]==Tan[Pi*cc0/2]
},{q0,q1,q2,q3,gam,cc},{t,t1}];

eps=0.0001;

f3[t_]:=Sqrt[M/(a^2(1-a^2))](EllipticPi[(a^2*k^2)/(a^2-1),JacobiAmplitude[a*t+EllipticK[k^2],k^2],k^2]-EllipticPi[(a^2*k^2)/(a^2-1),JacobiAmplitude[EllipticK[k^2],k^2],k^2]);

R[t_]:={{p1[0]Sqrt[1-a^2]/Sqrt[M-p0[0]^2],p2[0]/Sqrt[M-p0[0]^2],0},{-p2[0]/Sqrt[M-p0[0]^2],p1[0]Sqrt[1-a^2]/Sqrt[M-p0[0]^2],0},{0,0,1}}.

{{1,0,0},{0,p0[0]/Sqrt[M],Sqrt[(M-p0[0]^2)/M]},{0,-Sqrt[(M-p0[0]^2)/M],p0[0]/Sqrt[M]}}.

{{Cos[f3[t]],-Sin[f3[t]],0},{Sin[f3[t]],Cos[f3[t]],0},{0,0,1}}.

{{1,0,0},{0,p0[t]/Sqrt[M],-Sqrt[(M-p0[t]^2)/M]},{0,Sqrt[(M-p0[t]^2)/M],p0[t]/Sqrt[M]}}.

{{p1[t]Sqrt[1-a^2]/Sqrt[M-p0[t]^2],-p2[t]/Sqrt[M-p0[t]^2],0},{p2[t]/Sqrt[M-p0[t]^2],p1[t]Sqrt[1-a^2]/Sqrt[M-p0[t]^2],0},{0,0,1}}/.M->1-a^2(1-k^2);

p1[t_]:=-JacobiDN[a*t+EllipticK[k^2],k^2];
p2[t_]:=k*JacobiSN[a*t+EllipticK[k^2],k^2];
p0[t_]:=a*k*JacobiCN[a*t+EllipticK[k^2],k^2];

-(Sqrt[1-(a^2)]) (Sqrt[M]) Sin[f3[t]] p1[t]-Cos[f3[t]] p0[t] p2[t]/.M->1-a^2(1-k^2)//TraditionalForm

-(Sqrt[1-(a^2)]) Cos[f3[t]] p0[t] p1[t]+(Sqrt[M]) Sin[f3[t]] p2[t]/.M->1-a^2(1-k^2)//TraditionalForm

Almost-Riemmanian geometry 
of the two-sphere

I. Yu. Beschastnyi

Program Systems institute of RAS

International Conference “Geometry & control”
Moscow, 14-18 April, 2014

Almost-Riemannian structure on (S^2)

Let (OverVector[x])∈(S^2) be a point on a two-dimensional sphere and let (Subscript[(OverVector[b]), 1])∈(ℝ^3), (Subscript[(OverVector[b]), 2])∈(ℝ^3) be two noncollinear vectors.  Consider two vector fields (Subscript[X, 1]) and (Subscript[X, 2])∈Vec((S^2)) which are just infinitesimal rotations around this two vectors:

({{(Subscript[X, 1])(x)=(OverVector[x])×(Subscript[(OverVector[b]), 1]),          (Subscript[X, 2])(x)=(OverVector[x])×(Subscript[(OverVector[b]), 2]),            ‖x‖=1.}})

These vector fields span a distribution Δ:

({{({{(Subscript[Δ, x])=span}}){(Subscript[X, 1])(x),(Subscript[X, 2])(x) }.}})

Assume also that we have a scalar product g(·,·) on  Δ. Then informally the triple ((S^2),Δ,g) defines an almost-Riemannian structure on (S^2).

Why is it called almost-Riemannian?

Manipulate[Show[

Graphics3D[{Opacity->1,Lighter[White],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,1,0}]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,-1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{-1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,-1}]}],

If[XLines,ParametricPlot3D[Table[{Cos[a],Sin[a]*Cos[t],Sin[a]*Sin[t]},{a,0,Pi,Pi/10}],{t,0,2Pi},PlotStyle->{Blue}],nothing],

If[VFrame,Graphics3D[{Blue,Thick,Arrow[{{x[ϕ,ψ],y[ϕ,ψ],z[ϕ,ψ]},{x[ϕ,ψ],y[ϕ,ψ]+z[ϕ,ψ],z[ϕ,ψ]-y[ϕ,ψ]}}]}],nothing],

If[YLines,ParametricPlot3D[Table[{Sin[a]*Cos[t],Cos[a],Sin[a]*Sin[t]},{a,0,Pi,Pi/10}],{t,0,2Pi},PlotStyle->{Red}],nothing],

If[VFrame,Graphics3D[{Red,Thick,Arrow[{{x[ϕ,ψ],y[ϕ,ψ],z[ϕ,ψ]},{x[ϕ,ψ]-z[ϕ,ψ],y[ϕ,ψ],z[ϕ,ψ]+x[ϕ,ψ]}}]}],nothing],

If[TangentPlane,ParametricPlot3D[{x[ϕ,ψ]-u*z[ϕ,ψ],y[ϕ,ψ],z[ϕ,ψ]+u*x[ϕ,ψ]}+{0,v*z[ϕ,ψ],-v*y[ϕ,ψ]},{u,0,1},{v,0,1},Mesh->0,PlotStyle->{Opacity[0.2],Blue}],nothing],

If[Equator,ParametricPlot3D[RotationMatrix[t,{0,0,1}].{1,0,0},{t,0,2Pi},PlotStyle->{Thick}],nothing],

Boxed->False,Axes->False,PlotRange->All,ViewPoint->{1,-1,0.5},ImageSize->{650,550}],
Style["Point coordinates",Bold],{{ϕ,Pi/2},-Pi/2,Pi/2},{{ψ,0},-Pi,Pi},
Delimiter,
Style["Frames and Lines",Bold],Column[{Control[{{XLines,False,"Lines of (Subscript[X, 1])"},{False,True}}],Control[{{YLines,False,"Lines of (Subscript[X, 2])"},{False,True}}],Control[{{VFrame,False,"Show frame"},{False,True}}],Control[{{TangentPlane,False,Style["Tangent plane"]},{False,True}}],Control[{{Equator,False,"Show equator"},{False,True}}]},Right,BaselinePosition->Top],

TrackedSymbols:>{ϕ,ψ,XLines,YLines,VFrame,TangentPlane,Equator}, ControlPlacement->Left,LabelStyle->Large]

Outline of the talk

Slide

Formal definitions

Parameterization of geodesics

Neccesary optimality conditions

Formal definition of an almost-Riemannian structure on (S^2)

arrowShape=({Black,Arrowheads[0.1],Arrow[#,{20,20}]}&);
vertex={{Style["E",48],{0,140}},{Style["(TS^2)",48],{140,140}},{Style["(S^2)",48],{140,10}}};
edge={{1\[DirectedEdge]2,Placed[Style["f",24],{.5,{.5,-0.2}}],arrowShape},{2\[DirectedEdge]3,Placed[Style["π",24],{.5,{-1.2,.5}}],arrowShape},{1\[DirectedEdge]3,Placed[Style["(Subscript[π, E])",24],{.5,{1.2,1.2}}],arrowShape}};
graph=PlotDiagram[vertex,edge,AspectRatio->1,ImageSize->2*72,ImagePadding->{{15,40},{0,40}}]

E is a vector bundle of a rank 2 with a scalar product (Subscript[g, x])(· ,·) defined on (Subscript[E, x]) that depends smoothly on x∈(S^2);

f is a bundle morphism that induces an injective mapping from the set Γ(E) of all smooth sections of E to Vec((S^2)). We denote Δ=(SubStar[f])Γ(E)

The triple (E,f,g) defines an almost-Riemannian structure on (S^2).

The distance function

Let (Subscript[X, u])∈(Subscript[Δ, x]). Then we define a norm in (Subscript[Δ, x]) by

({{({{‖(Subscript[X, u])‖=min{(Subscript[g, x])(u,u), u∈(Subscript[E, q]) s.t. (Subscript[X, u])=f(u,x)}}})}})

A Lipschitz curve x(t) is said to be admissible if there exists a function u(t)∈(L^∞)([0,T],(Subscript[U, x(t)])) s.t.

({{({{(OverDot[x])(t)=f(u(t),x(t))}}).}})

The length of an admissible curve:

({{({{l(x(t))=(∫_0^T)‖(OverDot[x])(t)‖dt}})}})

The distance function:

({{({{d((Subscript[x, 0]),(Subscript[x, T]))=inf{l(x), x is admissible, x(0)=(Subscript[x, 0]),  x(T)=(Subscript[x, T])}}})}})

Optimal control problem

({{(OverDot[x])=x×ω,   ω=(Subscript[u, 1])(Subscript[b, 1])+(Subscript[u, 2])(Subscript[b, 2])}, {((Subscript[u, 1]),(Subscript[u, 2]))∈(ℝ^2);       x,(Subscript[b, 1]),(Subscript[b, 2])∈(ℝ^3);        ‖x‖=1;}, {x(0)=(Subscript[x, 0]),     x(T)=(Subscript[x, T]),}, {(∫_0^T)‖x×ω‖dt->min.}})

We’re going to use the isomorphism between (ℝ^3) and so(3). Namely if b∈(ℝ^3) and B∈so(3) are isomorphic then

({{({{b=(Subscript[b, 1])(Subscript[e, 1])+(Subscript[b, 2])(Subscript[e, 2])+(Subscript[b, 3])(Subscript[e, 3])≃(Subscript[b, 1])(Subscript[A, 1])+(Subscript[b, 2])(Subscript[A, 2])+(Subscript[b, 3])(Subscript[A, 3])}})=B}})

History of the problem

Slide

V.Grushin studied one of the first examples of almost-Riemannian structures now called the Grushin plane;

A.Agrachev, U.Boscain, G.Charlot, M.Sigalotti, R.Ghezzi (Ph.D. thesis), B.Bonnard, G.Janin, J.-B.Caillau studied almost-Riemannian structures on two-dimensional manifolds;

The almost-Riemannian structure (S^2) was explicitly studied in papers of U.Boscain, G.Charlot, M.Sigalotti and in papers of B.Bonnard, J.-B.Caillau.

Recent results in “B.Bonnard, O.Cots, J.-B.Pomet,            N.Sherbakova. Riemannian metrics on 2D manifolds related to the Euler-Poinsot rigid body motion”.

The lifted optimal control problem

Slide

Let Ω∈so(3) be isomorphic to ω∈ (ℝ^3). Then (Subscript[x, t])=(Subscript[R, t]^(-1))(Subscript[x, 0]) is a solution of (OverDot[x])=x×ω, where (Subscript[R, t]) satisfies

({{({{(OverDot[R])=RΩ=}})R((Subscript[u, 1])(Subscript[B, 1])+(Subscript[u, 2])(Subscript[B, 2])),}, {R∈SO(3),      Ω∈so(3),      ((Subscript[u, 1]),(Subscript[u, 2]))∈(ℝ^2),}, {R(0)=(e^(Subscript[βX, 0])),     R(T)=(e^(Subscript[βX, 0]))(Subscript[R, T]),}, {(∫_0^T)(Sqrt[g(Ω,Ω)])dt->min.}})

where (e^(Subscript[βX, 0])) is a rotation on angle βaround (Subscript[(OverVector[x]), 0]) and (Subscript[B, 1]),(Subscript[B, 2]) skew-symmetric matrices isomorphic to (Subscript[(OverVector[b]), 1]),(Subscript[(OverVector[b]), 2]).

Classification of sub-Riemannian structures on 3D Lie groups (A.Agrachev, D.Barilari)

Slide

Import["/ab_classification.png"]

Red points indicate sub-Riemannian structures, where minimal length curves are completely studied.

A family of optimal control problems

Slide

We normalize in such a way that κ+χ=1 and let 2χ=(a^2). Then we have a family of optimal control problems:

({{({{(OverDot[R])=RΩ=}})R((Subscript[u, 1])(Subscript[A, 2])+(Subscript[u, 2])(Sqrt[1-(a^2)])(Subscript[A, 1])),}, {R∈SO(3),     ((Subscript[u, 1]),(Subscript[u, 2]))∈(ℝ^2),      a∈(0,1)}, {R(0)=(e^(Subscript[βX, 0])),     R(T)=(e^(Subscript[βX, 0]))(Subscript[R, T]),}, {(∫_0^T)(Sqrt[(Subscript[u, 1]^2)+(Subscript[u, 2]^2)])dt->min.}})

If a=0, we get a surface of revolution. This case was studied in papers of U.Boscain, G.Charlot, M.Sigalotti and B.Bonnard, J.B.Caillau, O.Cots, J.-B.Pomet, N.Sherbakova.

Geodesic equations

Slide

After applying some first order necessary conditions, we get the following geodesics equations

({{({{(OverDot[R])=}})R((Subscript[p, 1])(Subscript[A, 2])+(Subscript[p, 2])(Sqrt[1-(a^2)])(Subscript[A, 1]))=RΩ,}, {(Subscript[(OverDot[p]), 1])=(Subscript[p, 2])(Subscript[p, 3]),}, {(Subscript[(OverDot[p]), 2])=-(Subscript[p, 1])(Subscript[p, 3]),
(Subscript[(OverDot[p]), 3])=(a^2)(Subscript[p, 1])(Subscript[p, 2]).}})

The vertical subsystem can be written in Lax form

({{({{(OverDot[P])=[P,Ω]}})}})

where P=(Subscript[p, 2])(Subscript[A, 1])+(Sqrt[1-(a^2)])(Subscript[p, 1])(Subscript[A, 2])+(Subscript[p, 3])(Subscript[A, 3])

Integration of the vertical subsystem

Slide

Substitution

({{({{(Subscript[p, 1])=cosγ,     (Subscript[p, 2])=-sinγ}})}})

reduces vertical subsystem to the equations of a mathematical pendulum

({{({{(OverDot[γ])=(Subscript[p, 3]),}, {(Subscript[(OverDot[p]), 3])=-((a^2)/2)sin 2γ.}})}})

Phase portrait of mathematical pendulum

Slide

Manipulate[StreamPlot[{y,-a^2Sin[2x]/2},{x,-Pi/2,3Pi/2},{y,-1.5,1.5},ImageSize->{480,480}],{{a,1/2},0,0.99999},LabelStyle->Large]

Integration of the horizontal part

Slide

Parametrization of the geodesics on SO(3) can be represented as follows

({{({{R=(e^(-(Subscript[ϕ, 1])(0)(Subscript[A, 3])))(e^(-(Subscript[ϕ, 2])(0)(Subscript[A, 1])))}})(e^((Subscript[ϕ, 3])(t)(Subscript[A, 3])))(e^((Subscript[ϕ, 2])(t)(Subscript[A, 1])))(e^((Subscript[ϕ, 1])(t)(Subscript[A, 3]))).}})

A classical approach that is very often used in classical mechanics is to get expressions for cos (Subscript[ϕ, 1]), sin (Subscript[ϕ, 1]), cos (Subscript[ϕ, 2]), sin (Subscript[ϕ, 2]) in terms of (Subscript[p, i])

({{({{cos (Subscript[ϕ, 2]) = ((Subscript[p, 3])/(‖p‖)),     }})                ({{sin (Subscript[ϕ, 2]) = ((Sqrt[((‖p‖)^2)-(Subscript[p, 3]^2)])/(‖p‖)),}})
cos (Subscript[ϕ, 1]) = (((Subscript[p, 1])(Sqrt[1-(a^2)]))/(Sqrt[((‖p‖)^2)-(Subscript[p, 0]^2)])),       sin (Subscript[ϕ, 1]) = ((Subscript[p, 2])/(Sqrt[((‖p‖)^2)-(Subscript[p, 0]^2)]))}})

and then solve the differential equations for  (Subscript[ϕ, 3]).

Transversality conditions

Slide

One can proof that, if the transversality conditions are satisfied in the source, then they are automatically satisfied in the target. Since the transversality conditions are also left-invariant we only need to check them at the identity element. So:

({{({{⟨(Subscript[p, 0]),T((e^(Subscript[βX, 0])))⟩=0}})⟺-sin (Subscript[γ, 0])(Subscript[x, 0])+cos (Subscript[γ, 0])(Subscript[y, 0])(Sqrt[1-(a^2)])+(Subscript[p, 3])(0)(Subscript[z, 0])=0}})

Transversality conditions

Slide

Manipulate[Row[{Show[StreamPlot[{y,-a^2Sin[2x]/2},{x,-Pi/2,3Pi/2},{y,-1.5,1.5}],ContourPlot[-Sin[γ]x[coor[[1]],coor[[2]]]+Cos[γ]y[coor[[1]],coor[[2]]]Sqrt[1-a^2]+p3*z[coor[[1]],coor[[2]]]==0,{γ,-Pi/2,3Pi/2},{p3,-1.5,1.5},ContourStyle->{Red,Thick}],ImageSize->{370,370}],

Show[Graphics3D[{Opacity[0.2],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,-1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
Graphics3D[{PointSize[Large],Point[{x[coor[[1]],coor[[2]]],y[coor[[1]],coor[[2]]],z[coor[[1]],coor[[2]]]}]}],Boxed->False,ImageSize->{420,420}]}],{{a,0.5},0,0.9999},{{coor,{0,-Pi/2},Style["Coordinates"]},{0,-Pi/2},{2Pi,Pi/2}},LabelStyle->Large]

Examples of geodesics

Slide

Manipulate[Row[{Show[StreamPlot[{c,-Sin[al]Sin[2gam]},{gam,-Pi/2,3Pi/2},{c,-3,3}],Plot[x[coor[[1]],coor[[2]]]*Sqrt[Cos[al]+Sin[al]]Sin[gam0]/z[coor[[1]],coor[[2]]]-y[coor[[1]],coor[[2]]]*Sqrt[Cos[al]-Sin[al]]Cos[gam0]/z[coor[[1]],coor[[2]]],{gam0,-Pi/2,3Pi/2},PlotStyle->{Thick,Red}],
Graphics[{PointSize[Large],Point[{gam,x[coor[[1]],coor[[2]]]*Sqrt[Cos[al]+Sin[al]]Sin[gam]/z[coor[[1]],coor[[2]]]-y[coor[[1]],coor[[2]]]*Sqrt[Cos[al]-Sin[al]]Cos[gam]/z[coor[[1]],coor[[2]]]}]}],ImageSize->{390,390}],Show[ParametricPlot3D[Evaluate[Table[{x[coor[[1]],coor[[2]]] ((q0[t])^2)+x[coor[[1]],coor[[2]]] ((q1[t])^2)-2 q0[t] (z[coor[[1]],coor[[2]]]q2[t]-y[coor[[1]],coor[[2]]] q3[t])+2 q1[t] (y[coor[[1]],coor[[2]]] q2[t]+z[coor[[1]],coor[[2]]] q3[t])-x[coor[[1]],coor[[2]]] (((q2[t])^2)+((q3[t])^2)),y[coor[[1]],coor[[2]]] ((q0[t])^2)+2 z[coor[[1]],coor[[2]]] q0[t] q1[t]-y[coor[[1]],coor[[2]]] ((q1[t])^2)+2 x[coor[[1]],coor[[2]]] q1[t] q2[t]+y[coor[[1]],coor[[2]]] ((q2[t])^2)+2 (-x[coor[[1]],coor[[2]]] q0[t]+z[coor[[1]],coor[[2]]] q2[t]) q3[t]-y[coor[[1]],coor[[2]]] ((q3[t])^2),z[coor[[1]],coor[[2]]] ((q0[t])^2)-2 q0[t] (y[coor[[1]],coor[[2]]] q1[t]-x[coor[[1]],coor[[2]]] q2[t])-z[coor[[1]],coor[[2]]] (((q1[t])^2)+((q2[t])^2))+2 (x[coor[[1]],coor[[2]]] q1[t]+y[coor[[1]],coor[[2]]] q2[t]) q3[t]+z[coor[[1]],coor[[2]]] ((q3[t])^2)}/.s2[al,gam0,coor[[2]],coor[[1]],t1],{gam0,{gam}}]],{t,0,t1}],
Graphics3D[{Opacity[0.1],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,-1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
Graphics3D[{PointSize[Large],Point[{x[coor[[1]],coor[[2]]],y[coor[[1]],coor[[2]]],z[coor[[1]],coor[[2]]]}]}],
PlotRange->All,Boxed->False,Axes->False,ImageSize->{440,440}]}],Row[{Column[{Control[{{al,Pi/8,Style["α",Italic,15]},0,Pi/4-eps,ImageSize->Medium,Appearance->"Labeled"}],Control[{{t1,2,Style["(Subscript[t, 1])",Italic,15]},eps,50,ImageSize->Medium,Appearance->"Labeled"}],Control[{{gam,0,Style["γ",Italic,15]},-Pi/2,3Pi/2,ImageSize->Medium,Appearance->"Labeled"}]}],Column[{Control[{{coor,{0.1,-Pi/2+0.1},Style["Coordinates"]},{0.01,-Pi/2+eps},{2Pi-0.01,Pi/2-eps}}]}]}],LabelStyle->Large]

Discrete symmetries of the Hamiltonian system

Slide

One can estimate the cut time by finding Maxwell points, i.e. points where two distinct geodesics of the same length meet one another. After a Maxwell point a geodesic is not globally optimal.

We can rewrite the Hamiltonian system in the following way

({{{({{(OverDot[x])=x×ω,}, {(OverDot[p])=p×ω;}})}})

where ω=(Subscript[p, 1])(Subscript[e, 2])+(Subscript[p, 2])(Sqrt[1-(a^2)])(Subscript[e, 1]) and p=(Subscript[p, 2])(Subscript[e, 1])+(Sqrt[1-(a^2)])(Subscript[p, 1])(Subscript[e, 2])+(Subscript[p, 3])(Subscript[e, 3]).

Discrete symmetries of the Hamiltonian system

Slide

Let (Subscript[I, j])=(e^(Subscript[πA, j])). It is easy to check that the following mappings are symmetries of the Hamiltonian system

({{{({{(Subscript[p, s])↦(Subscript[I, j])(Subscript[p, s,])}, {(Subscript[x, s])↦±(Subscript[I, j])(Subscript[x, s]);}})                        {({{(Subscript[p, s])↦-(Subscript[I, j])(Subscript[p, t-s,])}, {(Subscript[x, s])↦±(Subscript[I, j])(Subscript[x, t-s]);}})}})

Import["/reflections.png"]

Discrete symmetries of the exponential map

Slide

These mappings will be symmetries of the exponential map if and only if two conditions are satisfied:

these symmetries must be consistent with the transversality conditions;

these symmetries must fix the initial point.

The first condition is satisfied for all considered symmetries:

({{({{⟨(Subscript[I, j])(Subscript[p, t-s]),±(Subscript[I, j])(Subscript[x, t-s])⟩=}})±⟨(Subscript[p, t-s]),(Subscript[x, t-s])⟩;}, {±(d/ds)⟨(Subscript[p, t-s]),(Subscript[x, t-s])⟩=∓(⟨(Subscript[p, t-s])×(Subscript[ω, t-s]),(Subscript[x, t-s])⟩+⟨(Subscript[p, t-s]),(Subscript[x, t-s])×(Subscript[ω, t-s])⟩)=0;}, {⟨(Subscript[p, 0]),(Subscript[x, 0])⟩=0=>⟨(Subscript[p, t]),(Subscript[x, t])⟩=0.}})

The second condition is satisfied only for some sets that depend on the symmetry itself.

Characterization of Maxwell sets

Slide

Proposition. Maxwell sets for a point (Subscript[x, 0]) are great circles that are obtained by intersecting the coordinate planes with the sphere and which contain (Subscript[x, 0]). If a geodesic crosses a Maxwell set, it is not optimal anymore.

Manipulate[
Show[
Graphics3D[{PointSize[Large],Point[{x[ϕ,ψ]1.01,y[ϕ,ψ]1.01,z[ϕ,ψ]1.01}]}],
If[Abs[x[ϕ,ψ]]<0.05,ParametricPlot3D[RotationMatrix[t,{1,0,0}].{0,1,0},{t,0,2Pi},PlotStyle->{Red,Thick}],ParametricPlot3D[RotationMatrix[t,{1,0,0}].{0,1,0},{t,0,2Pi},PlotStyle->Opacity[0.1]]],
If[Abs[y[ϕ,ψ]]<0.05,ParametricPlot3D[RotationMatrix[t,{0,1,0}].{0,0,1},{t,0,2Pi},PlotStyle->{Red,Thick}],ParametricPlot3D[RotationMatrix[t,{0,1,0}].{0,0,1},{t,0,2Pi},PlotStyle->Opacity[0.1]]],
If[Abs[z[ϕ,ψ]]<0.05,ParametricPlot3D[RotationMatrix[t,{0,0,1}].{1,0,0},{t,0,2Pi},PlotStyle->{Red,Thick}],ParametricPlot3D[RotationMatrix[t,{0,0,1}].{1,0,0},{t,0,2Pi},PlotStyle->Opacity[0.1]]],
Graphics3D[{Opacity[0.1],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,-1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
ImageSize->{500,500},
Boxed->False,Axes->False
],{{ϕ,Pi/4},-Pi/2,Pi/2},{{ψ,0},-Pi,Pi}]

An equation for the Maxwell time

Slide

a (Sqrt[1-(k^2)]) (k^2) (JacobiCD[a t, (k^2)]) (JacobiSD[a t, (k^2)]) cos((Sqrt[((1-(a^2) (1-(k^2)))/((a^2) (1-(a^2))))]) ((EllipticPi3[(((a^2) (k^2))/((a^2)-1)), (JacobiAmplitude[a t+(EllipticK[(k^2)]), (k^2)]), (k^2)])-(EllipticPi[(((a^2) (k^2))/((a^2)-1)), (k^2)])))+(Sqrt[1-(a^2)]) (Sqrt[1-(k^2)]) (Sqrt[1-(a^2) (1-(k^2))]) (JacobiND[a t, (k^2)]) sin((Sqrt[((1-(a^2) (1-(k^2)))/((a^2) (1-(a^2))))]) ((EllipticPi3[(((a^2) (k^2))/((a^2)-1)), (JacobiAmplitude[a t+(EllipticK[(k^2)]), (k^2)]), (k^2)])-(EllipticPi[(((a^2) (k^2))/((a^2)-1)), (k^2)])))=0

k (Sqrt[1-(a^2) (1-(k^2))]) (JacobiCD[a t, (k^2)]) sin((Sqrt[((1-(a^2) (1-(k^2)))/((a^2) (1-(a^2))))]) ((EllipticPi3[(((a^2) (k^2))/((a^2)-1)), (JacobiAmplitude[a t+(EllipticK[(k^2)]), (k^2)]), (k^2)])-(EllipticPi[(((a^2) (k^2))/((a^2)-1)), (k^2)])))-a (Sqrt[1-(a^2)]) k (1-(k^2)) (JacobiND[a t, (k^2)]) (JacobiSD[a t, (k^2)]) cos((Sqrt[((1-(a^2) (1-(k^2)))/((a^2) (1-(a^2))))]) ((EllipticPi3[(((a^2) (k^2))/((a^2)-1)), (JacobiAmplitude[a t+(EllipticK[(k^2)]), (k^2)]), (k^2)])-(EllipticPi[(((a^2) (k^2))/((a^2)-1)), (k^2)])))=0

Case (Subscript[z, 0])=0

Slide

Corollary. A geodesic that starts from a singular set and crosses it at least ones is not optimal. From the parametrization we get the following bounds on the cut time

In (Subscript[C, 1]): (Subscript[t, cut])≤((2K(k))/a);

In (Subscript[C, 2]): (Subscript[t, cut])≤((2k K(k))/a);

In (Subscript[C, 3]): (Subscript[t, cut])≤(π/(Sqrt[1-(a^2)]));

In (Subscript[C, 4]): (Subscript[t, cut])≤π;

In (Subscript[C, 5]): (Subscript[t, cut])=(π/(Sqrt[1-(a^2)])).

Case (Subscript[x, 0])=±1 and (Subscript[y, 0])=±1

Slide

Proposition. If (Subscript[x, 0])=±1 or (Subscript[y, 0])=±1 and the corresponding geodesics do not leave the quarter-spheres limited by the Maxwell sets, then they do not have any cut points that are also Maxwell points.

Manipulate[Row[{Show[
ParametricPlot3D[RotationMatrix[s,{0,0,1}].{0,1,0},{s,0,2Pi},PlotStyle->{Thick}],
If[sym==1,ParametricPlot3D[Evaluate[Table[{x[θ,ϕ] ((q0[t])^2)+x[θ,ϕ] ((q1[t])^2)-2 q0[t] (z[θ,ϕ]q2[t]-y[θ,ϕ] q3[t])+2 q1[t] (y[θ,ϕ] q2[t]+z[θ,ϕ] q3[t])-x[θ,ϕ] (((q2[t])^2)+((q3[t])^2)),y[θ,ϕ] ((q0[t])^2)+2 z[θ,ϕ] q0[t] q1[t]-y[θ,ϕ] ((q1[t])^2)+2 x[θ,ϕ] q1[t] q2[t]+y[θ,ϕ] ((q2[t])^2)+2 (-x[θ,ϕ] q0[t]+z[θ,ϕ] q2[t]) q3[t]-y[θ,ϕ] ((q3[t])^2),z[θ,ϕ] ((q0[t])^2)-2 q0[t] (y[θ,ϕ] q1[t]-x[θ,ϕ] q2[t])-z[θ,ϕ] (((q1[t])^2)+((q2[t])^2))+2 (x[θ,ϕ] q1[t]+y[θ,ϕ] q2[t]) q3[t]+z[θ,ϕ] ((q3[t])^2)}/.s3[al,cc0,ϕ,θ,t1,1],{cc0,-sym*0.9,0.9,0.05}]],{t,0,t1}],nothing],
ParametricPlot3D[Evaluate[Table[{x[θ,ϕ] ((q0[t])^2)+x[θ,ϕ] ((q1[t])^2)-2 q0[t] (z[θ,ϕ]q2[t]-y[θ,ϕ] q3[t])+2 q1[t] (y[θ,ϕ] q2[t]+z[θ,ϕ] q3[t])-x[θ,ϕ] (((q2[t])^2)+((q3[t])^2)),y[θ,ϕ] ((q0[t])^2)+2 z[θ,ϕ] q0[t] q1[t]-y[θ,ϕ] ((q1[t])^2)+2 x[θ,ϕ] q1[t] q2[t]+y[θ,ϕ] ((q2[t])^2)+2 (-x[θ,ϕ] q0[t]+z[θ,ϕ] q2[t]) q3[t]-y[θ,ϕ] ((q3[t])^2),z[θ,ϕ] ((q0[t])^2)-2 q0[t] (y[θ,ϕ] q1[t]-x[θ,ϕ] q2[t])-z[θ,ϕ] (((q1[t])^2)+((q2[t])^2))+2 (x[θ,ϕ] q1[t]+y[θ,ϕ] q2[t]) q3[t]+z[θ,ϕ] ((q3[t])^2)}/.s3[al,cc0,ϕ,θ,t1,0],{cc0,-sym*0.9,0.9,0.05}]],{t,0,t1}],
Graphics3D[{Opacity[0.1],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
PlotRange->All,Boxed->False,Axes->False,ImageSize->{500,500}]}],Row[{Column[{Control[{{al,Pi/8,Style["α",Italic,15]},0,Pi/4-eps,ImageSize->Medium,Appearance->"Labeled"}],Control[{{t1,2,Style["(Subscript[t, 1])",Italic,15]},eps,Pi/Sqrt[1-2Sin[al]],ImageSize->Medium,Appearance->"Labeled"}]}],Column[{Control[{{θ,Pi,Style["θ",Italic,15]},0,2Pi,ImageSize->Medium,Appearance->"Labeled"}],Control[{{ϕ,0,Style["ϕ",Italic,15]},-Pi/2,Pi/2,ImageSize->Medium,Appearance->"Labeled"}],Control[{sym,{1,0}}]}]}]]

Case (Subscript[z, 0])=±1

Slide

Proposition. If  (Subscript[z, 0])=±1, the corresponding geodesics do not leave the quarter-spheres limited by the Maxwell sets and do not cross the singular set, then they do not have any cut points that are also Maxwell points.

Manipulate[Row[{Show[ParametricPlot3D[RotationMatrix[t,{0,0,1}].{1,0,0},{t,0,2Pi},PlotStyle->Thick],ParametricPlot3D[Evaluate[Table[{x[coor[[1]],coor[[2]]] ((q0[t])^2)+x[coor[[1]],coor[[2]]] ((q1[t])^2)-2 q0[t] (z[coor[[1]],coor[[2]]]q2[t]-y[coor[[1]],coor[[2]]] q3[t])+2 q1[t] (y[coor[[1]],coor[[2]]] q2[t]+z[coor[[1]],coor[[2]]] q3[t])-x[coor[[1]],coor[[2]]] (((q2[t])^2)+((q3[t])^2)),y[coor[[1]],coor[[2]]] ((q0[t])^2)+2 z[coor[[1]],coor[[2]]] q0[t] q1[t]-y[coor[[1]],coor[[2]]] ((q1[t])^2)+2 x[coor[[1]],coor[[2]]] q1[t] q2[t]+y[coor[[1]],coor[[2]]] ((q2[t])^2)+2 (-x[coor[[1]],coor[[2]]] q0[t]+z[coor[[1]],coor[[2]]] q2[t]) q3[t]-y[coor[[1]],coor[[2]]] ((q3[t])^2),z[coor[[1]],coor[[2]]] ((q0[t])^2)-2 q0[t] (y[coor[[1]],coor[[2]]] q1[t]-x[coor[[1]],coor[[2]]] q2[t])-z[coor[[1]],coor[[2]]] (((q1[t])^2)+((q2[t])^2))+2 (x[coor[[1]],coor[[2]]] q1[t]+y[coor[[1]],coor[[2]]] q2[t]) q3[t]+z[coor[[1]],coor[[2]]] ((q3[t])^2)}/.s2[al,gam0,coor[[2]],coor[[1]],t1],{gam0,(-1)^sym2*Pi/2,Pi+sym2*Pi/2,0.1}]],{t,0,t1}],
Graphics3D[{Opacity[0.1],Sphere[]}],
Graphics3D[{PointSize[Large],Red,Point[1.01{0,-1,0}]}],
Graphics3D[{PointSize[Large],Blue,Point[1.01{1,0,0}]}],
Graphics3D[{PointSize[Large],Green,Point[1.01{0,0,1}]}],
Graphics3D[{PointSize[Large],Point[{x[coor[[1]],coor[[2]]],y[coor[[1]],coor[[2]]],z[coor[[1]],coor[[2]]]}]}],
PlotRange->All,Boxed->False,Axes->False,ImageSize->{500,500}]}],Row[{Column[{Control[{{al,0.05,Style["α",Italic,15]},0,Pi/4-eps,ImageSize->Medium,Appearance->"Labeled"}],Control[{{t1,2,Style["(Subscript[t, 1])",Italic,15]},eps,6,ImageSize->Medium,Appearance->"Labeled"}]}],Column[{Control[{{coor,{Pi/2,0},Style["Coordinates"]},{0,-Pi/2},{2Pi,Pi/2}}]}],Column[{Control[{{sym2,1,Style["             sym"]},{1,0}}]},Alignment->Right]}],Alignment->Center]

Results

Slides

A full parameterization of geodesics is obtained

Symmetries and Maxwell sets are studied

Bounds for the cut time are obtained

Notes:

In the article of B.Bonnard, O.Cots, J.-B.Pomet and N.Sherbakova a parameterization of geodesics was obtained. They have also proven that a point is always optimal until it crosses the singular set.

Thank you for your attention