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}

(Source: https://mathworld.wolfram.com/ParaboloidGeodesic.html)

Load the required package to do variational calculus:


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?

  • $\begingroup$ 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. $\endgroup$
    – Narasimham
    Commented Oct 3, 2020 at 16:37
  • $\begingroup$ 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$ $\endgroup$
    – Narasimham
    Commented Oct 3, 2020 at 22:02

1 Answer 1


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}


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}]

enter image description here

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}

enter image description here

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}

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.