Find geodesics given two points

Suppose I want to find the geodesic on a paraboloid going through two points located on the paraboloid, using Mathematica 12.

I chose the paraboloid parametrization as follows:

paraboloid[{u_, v_}] := {Sqrt[u] Sin[v], Sqrt[u] Cos[v], u}


Load the required package to do variational calculus:

Needs["VariationalMethods"]


Then I set up the Euler-Lagrange equations like so:

eq = EulerEquations[Sqrt[Total[D[paraboloid[{u, v[u]}], u]^2]], v[u],u]


The output looks quite involved, so I was surprised that DSolve could obtain an exact solution using

DSolve[eq, v[u], u]


The two solutions contains ArcTanh, ArcSinh and quite a few roots (one solution below):

v[u_] := (u (ArcTanh[(2 Sqrt[1 + 4 u])/Sqrt[4 - u C[1]]] Sqrt[C[1]] Sqrt[4 - u C[1]] -
4 ArcSinh[(Sqrt[1 + 4 u] Sqrt[C[1]])/Sqrt[-16 - C[1]]] Sqrt[-16 -
C[1]] Sqrt[(4 - u C[1])/(16 + C[1])]))/(Sqrt[C[1]] Sqrt[u^2 (-4 + u C[1])]) + C[2]


Since I am interested in a specific geodesic passing through two given points, e.g. $$p_1=(u_1,v_1)=(1,1)$$ and $$p_2=(u_2,v_2)=(-1,-1)$$ on the paraboloid, I tried to determine the values of C[1] and C[2] like so:

Solve[{v[1] == 1, v[-1] == -1}, {C[1], C[2]}]


Neither Solve nor NSolve would return for at least an hour which is when I stopped the kernel, so how would I go about finding the specific solution making the geodesic go through $$(u_1,v_1)$$ and $$(u_2,v_2)$$?

EDIT: To clarify, I would like to know whether solving for C[1] and C[2] the way I tried is the proper way to obtain the geodesic connecting the points on the given surface. If the principle is sound, how can I let MMA arrive at a (numeric) solution?

UPDATE: I realized I made a mistake: The paraboloid parametrization at the top makes use of cylinder coordinates, therefore $$u$$ and $$v$$ represent the radius and the azimuth angle, respectively. So in order to trace out a complete, non-degenerate paraboloid, we need $$u > 0$$ and $$v \in [0,2\pi]$$. That means the second point through which the specific geodesic should pass cannot, by definition, be $$p_2=(-1,-1)$$ as above. $$(2,\frac{3\pi}{4})$$ is inside the intervals for $$u$$ and $$v$$, so that would be OK.

Several posts suggested using FindRoot to solve (systems of) transcendental equations like this one, so I took the solution to the E-L equation from above and modified the parameter list to include the integration constants $$C_1$$ and $$C_2$$:

v[u_, C1_, C2_] := (u (ArcTanh[(2 Sqrt[1 + 4 u])/Sqrt[4 - uC1]]Sqrt[C1]Sqrt[4 - u C1] - 4 ArcSinh[(Sqrt[1 + 4 u] Sqrt[C1])/Sqrt[-16 - C1]] Sqrt[-16 - C1] Sqrt[(4 - u C1)/(16 + C1)]))/(Sqrt[C1] Sqrt[u^2 (-4 + u C1)]) + C2


Find $$C_1$$ and $$C_2$$ to get the geodesic through $$p_1$$ and $$p_2$$: Let $$p_1=(1,0)$$, $$p_2=(2,\frac{3\pi}{4})$$, so

FindRoot[{v[1, C1, C2] == 0, v[2, C1, C2] == (3 \[Pi])/4}, {{C1, 1}, {C2, 1}}]


FindRoot returns {C1 -> -21.7912 - 1.99114*10^-14 I, C2 -> -1.52118 - 0.824159 I}. I suppose I may safely assume the imaginary part of $$C_1$$ to be zero, but not the one in $$C_2$$. Also, it says

FindRoot: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to ...


I looked at quite a number of network posts that reported the same issue, most of them seem to have specific workarounds which don't attack the problem simply by increasing WorkingPrecision - how can I get a correct solution?

• Between two points there are many local geodesics on the surface of a paraboloid. However there exists a unique path among them that globally minimizes path length. Oct 3 '20 at 16:37
• Another way. (r,z) are known. in a boundary value problem. For zero geodesic curvature after going through first given start point, by "shoot through iteration " .... of changing arc length until the geodesic passes through the second point of given (r,z).... can fix the unknown $\theta$ Oct 3 '20 at 22:02

Looks like you will only find solutions for $$(u_1,v_1)=(-1,1)$$ and $$(u_2,v_2)=(-1,-1)$$.

Clear[c, d, u, v]
paraboloid[{u_, v_}] := {Sqrt[u] Sin[v], Sqrt[u] Cos[v], u}

Needs["VariationalMethods"]

eq = EulerEquations[Sqrt[Total[D[paraboloid[{u, v[u]}], u]^2]], v[u], u];
sol = DSolve[eq, v[u], u] /. {C[1] -> c, C[2] -> d};
v[u_, c_, d_] := Evaluate[sol[[1, 1, 2]]]

Manipulate[Plot[v[u, c, d], {u, -3, 3},
PlotRange -> {{-3, 3}, {-3, 3}}], {c, -10, 10}, {d, -10, 10}]


sol2 = Solve[With[{u = -1}, v[u, c, d] == -1], d];
c = 0.18;
d /. sol2

{-0.324666 + 0. I}


In other words, -0.324666

Likewise for

v[u_, c_, d_] := Evaluate[sol[[2, 1, 2]]];
sol2 = Solve[With[{u = -1}, v[u, c, d] == 1], d]
c = 0.18;
d /. sol2

{0.324666 + 0. I}


Also for

c = 1;
d /. sol2

{0.363122 + 0. I}


and so on.

The imaginary component becomes significant when u = -1 cannot be reached, e.g.

c = -6;
d /. sol2

{4.55771 - 0.944697 I}