ABSTRACT (original article): We consider higher-dimensional generalizations of the -Grushin plane, focusing on the problem of classification of geodesics that minimize length, also known as optimal synthesis. Solving Hamilton’s equations on these spaces using the calculus of generalized trigonometric functions, we obtain explicit conjugate times for geodesics starting at a Riemannian point. From symmetries in the geodesic structure, we propose a conjectured cut time, and prove that it provides an upper bound on the first conjugate time, a key step in the extended Hadamard technique. In the three-dimensional case, we combine this method with a density argument to establish the conjecture. CITATION (original article): Michael Albert, Samuël Borza, Maria Gordina (2025), Geodesics on Grushin spaces, arXiv:2509.03411. https://doi.org/10.48550/arXiv.2509.03411
This notebook contains the code necessary to produce the simulations and graphics for the paper “Geodesics on Grushin Spaces” by Albert, Borza and Gordina, which can be found on the ArXiv. The basic object of the paper is the 3D stratified Grushin space. It is a sub-Riemannian manifold. For an excellent introduction to sub-Riemannian geometry, in particular the Hamiltonian theory that allows computation of geodesics, see the book of Agrachev, Barilari and Boscain titled “A Comprehensive Introduction to Geometry”.
The Grushin space G_(α,β) is simply 3D Euclidean space, but endowed with a special metric, arising by taking the minimum value of the sub-Riemannian length functional over all admissible curves connecting the endpoints. More precisely, we consider X=, Y= , , which form an orthonormal frame of rank 3 everywhere except for the union of planes {y=0} and {x=0}. An admissible curve γ is an absolutely continuous curve whose tangent vectors lie in Span(X,Y,Z) for almost all moments in time. Referring back to X,Y,Z, notice that if we take enough iterated Lie brackets, then we can recover the full tangent space at all points. This is the substance of Hormander’s condition. By virtue of the Chow-Rashevskii Theorem, the Carnot-Cartheodory distance defined by d_{CC}(x,y)=inf{T>0: L(γ), γ is an admissible curve connecting x and y that is arc length parametrized} is well defined. The length minimizers of this distance are called geodesics (although not all geodesics are length minimizing for their full duration), and they are obtained by solving Hamilton’s equations on the Cotangent bundle. In our paper we use the notation (x,y,z) for the spatial variables and (u,v,w) for the momentum (covector) variables on each fiber. The Hamiltonian arising from this structure is H= 1/2 . The solutions to the corresponding Hamiltonian system can be derived exactly in terms of so-called generalized trigonometric functions. These are generalizations of the usual sine and cosine, satisfying a slightly modified second order ODE, that is perfectly suited for Mathematica’s numerical ODE solver capability. Storing the solutions as interpolating functions, we are able to proceed with a suite of analysis and visualization.
∂
x
α
x
∂
y
Z=
α
x
β
y
∂
z
++
2
u
α
x
2
v
α
x
β
y
2
w
Our goal is the so-called “optimal synthesis” of Grushin spaces. This relates to the so-called cut-time of geodesics. The cut time of a geodesic γ is the maximal time t_{cut} on which the geodesic is a length minimizer. The point γ(t_{cut}) is called a cut point to the initial point γ(0). The collection of all cut points is called the cut locus. We will attempt to visualize the structure of the cut locus using Mathematica.
We note that Mathematica does not help one prove anything regarding optimal synthesis, as the existence of symmetrizing geodesics needs to be truly ruled out analytically up to a candidate time. This is done in our paper by utilizing the extended Hadamard technique. However, Mathematica was essential in the early stages for visualizing what the cut locus ended up looking like, as the extended Hadamard technique relies on certain properties of the conjectured cut locus, such as simple connectedness, and in our case the fine properties of the intersecting regions of the conjectured cut locus was also significant in our analysis.
Code
Code
In[]:=
Clear["Global`*"]
This is some root finding technology that we got from a StackExchange post by user Jens. It is well suited for interpolating functions with simple roots.
In[]:=
Clear[findAllRoots]SyntaxInformation[findAllRoots]={"LocalVariables"->{"Plot",{2,2}},"ArgumentsPattern"->{_,_,OptionsPattern[]}};SetAttributes[findAllRoots,HoldAll];Options[findAllRoots]=Join[{"ShowPlot"->False,PlotRange->All},FilterRules[Options[Plot],Except[PlotRange]]];findAllRoots[fn_,{l_,lmin_,lmax_},opts:OptionsPattern[]]:=Module[{pl,p,x,localFunction,brackets},localFunction=ReleaseHold[Hold[fn]/.HoldPattern[l]:>x];If[lmin!=lmax,pl=Plot[localFunction,{x,lmin,lmax},Evaluate@FilterRules[Join[{opts},Options[findAllRoots]],Options[Plot]]];p=Cases[pl,Line[{x__}]:>x,Infinity];If[OptionValue["ShowPlot"],Print[Show[pl,PlotLabel->"Finding roots for this function",ImageSize->200,BaseStyle->{FontSize->8}]]],p={}];brackets=Map[First,Select[(*ThisSplittrickpretendsthattwopointsonthecurveare"equal"ifthefunctionvalueshave_opposite_sign.Pairsofsuchsign-changesformthebracketsforthesubsequentFindRoot*)Split[p,Sign[Last[#2]]==-Sign[Last[#1]]&],Length[#1]==2&],{2}];x/.Apply[FindRoot[localFunction==0,{x,##1}]&,brackets,{1}]/.x->{}](*ThisisacodethatwefoundonStackExchangefromuser"Jens",whichwillhelptofindthezeroesofinterpolatingfunctions*)
Hamiltonian System
Hamiltonian System
Here we define the Hamiltonian system for the 3D Grushin space. Mathematica can easily solve this system and stores the solutions as interpolating functions.
sol[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=NDSolvex'[t]==u[t],y'[t]==v[t],z'[t]==w[t],u'[t]==-a*+,v'[t]==-b*,w'[t]==0,x[0]==x0,y[0]==y0,z[0]==z0,u[0]==u0,v[0]==v0,w[0]==w0,{x,y,z,u,v,w},{t,-15,30};(*ThisistheHamiltoniansystemforgeodesicsinthe3DGrushinspace*)
2a
x[t]
2a
x[t]
2b
y[t]
2a-1
x[t]
2
v[t]
2
w[t]
2b
y[t]
2a
x[t]
2b-1
y[t]
2
w[t]
Here we define the generalized trigonometric function Sin_alpha, which is understood as the unique solution to the parameter dependent second order ODE x’’= , x’(0)=1, x(0)=0. We also define its derivative Cos_α = Sin_α’. Finally, the numerical constant π_α, written here as the function pia[α_] is important for the symmetry properties that these functions enjoy. In the following we will focus on the sub-case when the initial point (x0,y0,z0) is a Riemannian point.
-
2α-1
αx
In[]:=
sol2=ParametricNDSolve[{x''[t]==-ax[t],x[0]==0,x'[0]==1},x,{t,-20,20},{a}];Sina[α_]:=If[α>1,x[α]/.sol2,Sin];Cosa[α_]:=Sina[α]';η[α_,t_]:=t-Sina[α][t]*Cosa[α][t];(*Definingthegeneralizedtrigfunctionsandanauxillaryfunctionη*)pia[a_]:=2NIntegrate,{t,0,1},WorkingPrecision->50;(*Defininingaconstantthatdependsona*)
2a-2
Abs[x[t]]
1
1/2
(1-)
2a
t
The solution for x(t) is given by x(t)=A_1Sin_α(ω_1t+ϕ_1) for some parameters A_1, ω_1, ϕ_1 uniquely determined by the initial data. The solution for y(t) is given by y(t)=A_2Sin_β(ω_2∫x^2α(s)ds+ϕ_2) for some parameters A_2,ω_2,ϕ_2. The integral can be computed using the η function. Note that ∫Sin_α^2α(s)ds= 1/(α+1) η(t)+C. We make an initial conjecture for the cut time as tau_1 = pi_alpha/ |ω_1|. For an arc length parametrized geodesic (meaning that its covector lies on the energy shell H=1/2), the quantity ^2+y0^2βw0^2)^1/2α. One can check is indeed a symmetry for the geodesics using the symmetry properties of the alpha-sine and cosine functions Sin_α(x+π_α)=-Sin_α(x) (similar for cosine).
|ω_1|=(v0
In[]:=
CutTime1[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=If(v0==0)&&(w0==0),0,;(*Thiswillbeoneofthepossiblecuttimes,whichcorrespondstotheonefoundinBorza'spaper*)
pia[a]
1
2a
+
2
v0
2b
y0
2
w0
Another possible cut time is tau_2, which is essentially the first (or second; depending on the initial covector) time that y(t) intersects the plane y=-y0. In the paper it is defined as the minimum solution of a certain non-linear integral equation. This solution is a multiple root, which is why we add a very small offset of .0001. The findAllRoots function is really only proficient for simple roots. This could possibly be avoided by adding a condition on derivatives, but we felt this would slow things down too much.
Our overall conjectured cut time is the minimum of tau_1 and tau_2. We prove in the paper using the extended Hadamard technique that it is the true cut time. See sections 6 and 7.
Our overall conjectured cut time is the minimum of tau_1 and tau_2. We prove in the paper using the extended Hadamard technique that it is the true cut time. See sections 6 and 7.
In[]:=
CutTime2[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=If>(a+1),If[v0>=0,findAllRoots[(y[t]+(y0-.0001))/.(sol[x0,y0,z0,u0,v0,w0,a,b]),{t,0,30},"ShowPlot"->False][[1]],findAllRoots[(y[t]+(y0-.0001))/.(sol[x0,y0,z0,u0,v0,w0,a,b]),{t,0,30},"ShowPlot"->False][[2]]],30;CutTime2[1,1,0,1,-.05,.45,2,2](*Thisistheotherpossiblecuttime*)
2
w0
2b
pia[b]
pia[a]
b+1+b/a
+
2b
y0
2
w0
2
v0
Out[]=
1.38959
In[]:=
CutTime2Sing[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=If>(a+1),findAllRoots[y[t]/.(sol[x0,y0,z0,u0,v0,w0,a,b]),{t,0,30},"ShowPlot"->False][[1]],30;(*Definingthesecondcuttimeinthecasethaty0=0*)
2
w0
2b
pia[b]
pia[a]
b+1+b/a
+
2b
y0
2
w0
2
v0
In[]:=
xcutval[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=First[x[CutTime1[x0,y0,z0,u0,v0,w0,a,b]]/.(sol[x0,y0,z0,u0,v0,w0,a,b])];ycutval[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=First[y[CutTime1[x0,y0,z0,u0,v0,w0,a,b]]/.(sol[x0,y0,z0,u0,v0,w0,a,b])];zcutval[x0_,y0_,z0_,u0_,v0_,w0_,a_,b_]:=First[z[CutTime1[x0,y0,z0,u0,v0,w0,a,b]]/.(sol[x0,y0,z0,u0,v0,w0,a,b])];(*Evaluatingeachgeodesicatthecuttime*)
Inequalities and Type 1/ Type 2 points.
Inequalities and Type 1/ Type 2 points.
In the paper we derive some inequalities that tell us for a given energy shell constrained covector, which of tau_1 or tau_2 is smaller. By using the ImplicitRegion function, we can visualize these regions on the energy shell. We also make sure to shade the conjugate covectors for greater visualization clarity. The red region is where tau_2 is smaller than tau_1, while the light orange region is where tau_1 is smaller than tau_2.
There is a further consideration that needs to be made in terms of type 1 and type 2 initial points. Type 1 points satisfy pi_alpha |x0|^(alpha+1)/ pi_beta(alpha+1) <= |y0|. This is discussed in section 6 of the paper. Essentially this determines whether or not the topological boundaries of the conjectured cut loci within each plane {x=-x0} and {y=-y_0} intersect in a non-trivial way.
There is a further consideration that needs to be made in terms of type 1 and type 2 initial points. Type 1 points satisfy pi_alpha |x0|^(alpha+1)/ pi_beta(alpha+1) <= |y0|. This is discussed in section 6 of the paper. Essentially this determines whether or not the topological boundaries of the conjectured cut loci within each plane {x=-x0} and {y=-y_0} intersect in a non-trivial way.
T1T2[x0_,y0_,z0_,a_,b_]:=ImplicitRegion>=(a+1),{v0,w0}
2
w0
2b
pia[b]
pia[a]
b+1+
b
a
+
2b
y0
2
w0
2
v0
In[]:=
H12[x0_,y0_,z0_,a_,b_]:=Plot3DSqrt1--,v0,,,w0,,,RegionFunction->Function{v0,w0},+<=1,PlotStyle->Opacity[.6];H121[x0_,y0_,z0_,a_,b_]:=Plot3DSqrt1--,{v0,w0}∈T1T2[x0,y0,z0,a,b],PlotRange->Full,PlotStyle->{Red};H122[x0_,y0_,z0_,a_,b_]:=Plot3D-Sqrt1--,{v0,w0}∈T1T2[x0,y0,z0,a,b],PlotRange->Full,PlotStyle->{Blue};(*Describingtheinequalityoncovectorsthatgives>=*)
2a
x0
2
v0
2a
x0
2b
y0
2
w0
-1
a
Abs[x0]
1
a
Abs[x0]
-1
a
Abs[x0]
b
Abs[y0]
1
a
Abs[x0]
b
Abs[y0]
2a
x0
2
v0
2a
x0
2b
y0
2
w0
2a
x0
2
v0
2a
x0
2b
y0
2
w0
2a
x0
2
v0
2a
x0
2b
y0
2
w0
τ
1
τ
2
In[]:=
Show[H121[1,1,3,2,2],H12[1,1,3,2,2],ParametricPlot3D[{Cos[θ]Sin[π/8],Sin[θ]Sin[π/8],Cos[π/8]},{θ,0,2π},PlotStyle->{LightBlue,Dashed}],Graphics3D[{PointSize[.03],Opacity[1],Green,Point[{0,0,1.04}]}],ParametricPlot3D[{Cos[θ]Sin[π/4],Sin[θ]Sin[π/4],Cos[π/4]},{θ,0,2π},PlotStyle->{LightBlue,Dashed}],ParametricPlot3D[{Cos[θ],Sin[θ],0},{θ,0,2π},PlotStyle->{Thickness[.01],Blue}],ParametricPlot3D[{0,Sin[ϕ],Cos[ϕ]},{ϕ,-.495,.495},PlotStyle->{Thickness[.01],Black}],PlotRange->Full]
Out[]=
Here we visualize the non-trivial intersection of the topological boundaries of the conjectured cut loci. This is in the type 2 case, and we take alpha=beta= 1, where the geodesics are given in terms of sines and cosines.
Visualization of the Cut Locus
Visualization of the Cut Locus
Here we are defining functions that will allow us to visualize the topological boundary of the portion of the conjectured cut locus lying in the plane {x=-x0}. It is an egg shaped closed curve of discrete points. One should use this function only in the type 1 case.
Main Graphic of the Paper
Main Graphic of the Paper
This is the code used to create the main graphic for the cut locus of a type 1 Riemannian point. One can adjust alpha and beta or the initial values of x0, y0 to see how it affects the resulting geometry of the cut locus. One should ensure that they remain inside the region of type 1 points though: being sure to check that pi_alpha |x0|^(alpha+1)/ pi_beta(alpha+1) <= |y0|. Observe that R^3 without the conjectured cut locus is simply connected. This is essential for the extended Hadamard technique.
Cut Locus for a singular initial point
Cut Locus for a singular initial point
In the following, we also create graphics for the cut locus starting from singular points, featuring some geodesic pairs and their intersection points. There are 3 types of singular points. Either x0=0, y0 non-zero, or x0 non-zero and y0=0, or both x0, y0=0.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Geodesics on Grushin spaces
by Michael Albert, Samuël Borza, Maria Gordina
Wolfram Community, STAFF PICKS, December 23, 2025
https://community.wolfram.com/groups/-/m/t/3596261
by Michael Albert, Samuël Borza, Maria Gordina
Wolfram Community, STAFF PICKS, December 23, 2025
https://community.wolfram.com/groups/-/m/t/3596261