Using the Theory of Functional Connections to Solve Boundary Value Geodesic Problems

Abstract

This study provides a least-squares-based numerical approach to estimate the boundary value geodesic trajectory and associated parametric velocity on curved surfaces. The approach is based on the Theory of Functional Connections, an analytical framework to perform functional interpolation. Numerical examples are provided for a set of two-dimensional quadrics, including ellipsoid, elliptic hyperboloid, elliptic paraboloid, hyperbolic paraboloid, torus, one-sheeted hyperboloid, Moëbius strips, as well as on a generic surface. The estimated geodesic solutions for the tested surfaces are obtained with residuals at the machine-error level. In principle, the proposed approach can be applied to solve boundary value problems in more complex scenarios, such as on Riemannian manifolds.

1. Introduction

Geodesic trajectories represent the paths on a curved surface, whose acceleration has no component on the local tangent plane to the surface. In many common scenarios, a geodesic represents the straightest path between two points. In particular, in a Riemannian manifold, geodesics are characterized by the property of having no geodesic curvature [1,2].

Most of the general studies on geodesics’ trajectories are provided as initial value problems, while most of the boundary value problems are related to the biaxial and triaxial ellipsoid. In [3], the boundary value problem on an ellipsoid with boundary (Dirichlet) conditions is replaced by an initial value problem with Dirichlet and Neumann conditions. In particular, the Neumann condition is obtained iteratively by numerically integrating a system of four first-order differential equations. Triaxial ellipsoids are considered in [4,5], while [6] has provided an analytical approach to solve this boundary value problem with symmetry. Reference [7] contains, because of its importance to terrestrial geodesy, a rich literature survey on the geodesic equations for low eccentric ellipsoid in both Cartesian and polar coordinates. Since the problem of fitting ellipsoid is important in geodesy (all geodetic calculations are performed on a reference ellipsoid), Refs. [5,8] analyzes the computational differences in the fitting ellipsoid using biaxial ellipsoid instead of triaxial ellipsoid and by performing least-squares ellipsoid fitting. Noteworthy, Reference [7] has improved (in terms of computational time) the existing solutions using differential equations in Cartesian coordinates and Taylor series expansions by simplifying previous formulations.

In this study, the geodesic boundary value problem is numerically investigated for any curved surface using the general geodesic equations from differential geometry. These geodesic equations are then solved by nonlinear least-squares using the method that the Theory of Functional Connections [9] (TFC) has introduced to solve differential equations [10,11]. The existence of a solution for the geodesic boundary value problem (boundary value problems, in general, may have a unique solution, no solution, or infinite solutions) is guaranteed by the Hopf–Rinow theorem [12,13]:

If a length-metric space $(M,d)$ is complete and compact then any two points, $({\mathit{p}}_1,{\mathit{p}}_2\in M)$, can be connected by a minimizing geodesic.

In this theorem, M indicates a manifold and d the metric distance. Roughly speaking, complete indicates a space where there are no “points missing” from it (inside or at the boundary), while compact indicates a space that is closed (space bounds included) and bounded (distance between any two points limited).

Motivated by this theorem, this study proposes a general numerical and accurate approach to solve boundary value geodesic problems in curved surfaces. First, this study briefly provides the geodesic equations for the general Riemannian spaces, followed by a short background on TFC. Then, it shows how to apply TFC to solve boundary value geodesic problems by nonlinear least-squares. In particular, an ad hoc algorithm to avoid local minima is presented.

The proposed approach is then numerically tested on various quadric surfaces, such as triaxial ellipsoid, elliptic hyperboloid, elliptic paraboloid, hyperbolic paraboloid, torus, one-sheeted hyperboloid, Moëbius strip, and on a generic surface. In all these tests, machine error estimation of the geodesic trajectory is obtained along with the indirect numerical proof that the associated parametric velocity is constant.

2. Background on Geodesics Equations on Riemannian Manifold

A Riemannian manifold is a smooth curved space in which the infinitesimal distance, $\mathrm{d}s$, satisfies

$$\mathrm{d}{s}^2=g_{ij}\mathrm{d}{x}^i\mathrm{d}{x}^j$$

(1)

where $g_{ij}$ is the covariant metric tensor of the space and where Einstein’s notation has been used in Equation (1).

The metric tensor is the fundamental tool used in differential geometry to study curved spaces. Specifically, for Euclidean spaces, the metric tensor is diagonal ($g_{ij}=0$ if $i\ne j$) and, in particular, is the identity, $g_{ij}={\delta }_{ij}$, for the Cartesian metric tensor. The metric tensor, $g_{ij}$, is the matrix composed by all inner products between all partials of the vector defining the Riemannian surface, $\mathit{p}(x_1,x_2,\dots ,x_n)$

$$g_{ij}={\left({\displaystyle \frac{\partial \mathit{p}}{\partial x_i}}\right)}^{\mathrm{T}}{\displaystyle \frac{\partial \mathit{p}}{\partial x_j}}$$

In Riemannian geometry, geodesics are not the same as “shortest curves” between two points, though the two concepts are closely related. The difference is that geodesics are only locally the shortest distance between points. Going the “long way round” on a great circle between two points on a sphere is a geodesic but not the shortest path between the points. In addition, geodesic paths need not be unique.

On a curved surface, the length of a parametric trajectory, defined by the coordinates $x^i=x^i\left(t\right)$ and connecting the position vectors, $\mathit{p}\left(t_0\right)$ and $\mathit{p}\left(t_f\right)$, respectively, is,

$${\displaystyle L={\int }_{t_0}^{t_f}\mathrm{d}s={\int }_{t_0}^{t_f}\sqrt{g_{ij}{\displaystyle \frac{\mathrm{d}{x}^i}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{x}^j}{\mathrm{d}t}}}\mathrm{d}t}$$

(2)

This parametric trajectory has no normal acceleration if it satisfies the following n differential equations [1,2,14],

$${\displaystyle \frac{{\mathrm{d}}^2{x}^i}{\mathrm{d}{t}^2}}+{\mathsf{\Gamma }}_{jk}^i{\displaystyle \frac{\mathrm{d}{x}^j}{\mathrm{d}t}}{\displaystyle \frac{\mathrm{d}{x}^k}{\mathrm{d}t}}=0.$$

(3)

These n second-order ordinary differential equations are the geodesic equations. They define the geodesic trajectories on a manifold with metric tensor, $g_{ij}$. In Equation (3), the ${\mathsf{\Gamma }}_{jk}^i$ terms are the Christoffel symbols of second kind, defined as,

$${\mathsf{\Gamma }}_{jk}^i={\displaystyle \frac{1}{2}}{g}^{mi}\left({\displaystyle \frac{\partial g_{km}}{\partial x^j}}+{\displaystyle \frac{\partial g_{mj}}{\partial x^k}}-{\displaystyle \frac{\partial g_{jk}}{\partial x^m}}\right)$$

(4)

where $g^{mi}$ is the controvariant inverse of the covariant metric tensor (or conjugate or dual metric),

$$g^{mi}={\left(g_{mi}\right)}^{-1}$$

Note that the Christoffel symbols satisfy the relationship, ${\mathsf{\Gamma }}_{jk}^i={\mathsf{\Gamma }}_{kj}^i$ and are related to the Christoffel symbols of the first kind, $[m,jk]$, by the relationship, $[m,jk]=g_{mi}{\mathsf{\Gamma }}_{jk}^i$.

Geodetic (Parametric) Velocity

A curve on a surface is called geodesic if, at every point, the acceleration is either zero or parallel to the normal. A geodesic trajectory, $\mathit{p}\left(t\right)$, on a curved surface has constant parametric speed. The proof is immediate,

$${\displaystyle \frac{\mathrm{d}{\dot{p}}^2}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}\left({\dot{\mathit{p}}}^{\mathrm{T}}\dot{\mathit{p}}\right)}{\mathrm{d}t}}=2{\ddot{\mathit{p}}}^{\mathrm{T}}\dot{\mathit{p}}=0.$$

The word “velocity,” here, is not the instantaneous ratio between distance and time, but the instantaneous ratio between distance and the parameter t selected to describe the trajectory (which can also be the time).

Let us consider a two-dimensional surface in a three-dimensional space. Let the surface be described by the parametic equations,

$$x=f_x(u,v),y=f_y(u,v)\mathrm{and}z=f_z(u,v),$$

then, by computing,

$${\dot{p}}_x={\displaystyle \frac{\partial f_x}{\partial u}}\dot{u}+{\displaystyle \frac{\partial f_x}{\partial v}}\dot{v},{\dot{p}}_y={\displaystyle \frac{\partial f_y}{\partial u}}\dot{u}+{\displaystyle \frac{\partial f_y}{\partial v}}\dot{v},\mathrm{and}{\dot{p}}_z={\displaystyle \frac{\partial f_z}{\partial u}}\dot{u}+{\displaystyle \frac{\partial f_z}{\partial v}}\dot{v},$$

the velocity in a trajectory defined by $\dot{p}(u\left(t\right),v\left(t\right))$, is provided by ${\dot{p}}^2={\dot{p}}_x^2+{\dot{p}}_y^2+{\dot{p}}_z^2$.

To give a trivial example, the parametric description of a unit-radius sphere is provided by,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}sinucosv,& sinusinv,& cosu\end{array}\right\}$$

where $u\in [0,\pi ]$ and $v\in [0,2\pi )$. The covariant metric tensor and the geodesic equations for the sphere are,

$$g_{ij}=\left[\begin{array}{cc}1& 0\\ 0& {sin}^{2}u\end{array}\right]\mathrm{and}\left\{\begin{array}{c}{\mathcal{L}}_u=\ddot{u}-{\dot{v}}^{2}sinucosu=0\hfill \\ {\mathcal{L}}_v=\ddot{v}+2\dot{u}\dot{v}cotu=0\hfill \end{array}\right.,$$

respectively, and the geodetic parametric velocity on the sphere is,

$${\dot{p}}^2=2{\dot{u}}^2{cos}^{2}u+{\dot{v}}^2{sin}^{2}u.$$

Note that $u=0$ and $u=\pi$ make singular the $g_{ij}$ matrix and, consequently, make singular the geodesic equations. This singularity changes location if another set of polar coordinates is selected to describe a sphere.

3. Background on the Theory of Functional Connections

The Theory of Functional Connections (TFC) is a mathematical framework to perform functional interpolation [9,15]. This is performed by analytically deriving some functionals, called constrained expressions, representing all functions satisfying a set of linear constraints in n-dimensional space. This way, constrained optimization problems subject to linear constraints, such as differential equations, are transformed into unconstrained problems. The optimization problem consists then of deriving the expression of an unconstrained free function, $g\left(x\right)$ (Note that $g\left(x\right)$ can be discontinuous, partially defined, and even the Dirac delta function, as long as it is defined on where the constraints are defined.).

In general, for a univariate function, $y\left(x\right)$, subject to n linear constraints (e.g., points, derivatives, integrals, and any linear combination of them), two equivalent definitions of constrained expressions can be introduced [15,16,17],

$$\begin{array}{cc}\hfill y(x,g(x\left)\right)& =g\left(x\right)+\sum _{k=1}^n{\eta }_k(x,g\left(x\right))s_k\left(x\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill y(x,g(x\left)\right)& =g\left(x\right)+\sum _{k=1}^n{\rho }_k(x,g\left(x\right)){\varphi }_k(x,\mathbf{s}\left(x\right))\hfill \end{array}$$

(6)

In the formal definition of Equation (5), the ${\eta }_k(x,g\left(x\right))$ are functional coefficients whose expressions are derived by imposing the n constraints on (5), and the $s_k\left(x\right)$ are a set of n user-defined linearly independent support functions. In Equation (6), the ${\varphi }_k(x,\mathbf{s}\left(x\right))$ are switching functions which imply changing between the constraints, and the ${\rho }_k(x,g\left(x\right))$ are projection functionals, which project the free function $g\left(x\right)$ to the kth constraint. See Ref. [9] for a complete and detailed explanation of these terms.

For example, consider a function $y\left(x\right)$ subject to the constraints

$${\displaystyle \frac{\mathrm{d}y}{\mathrm{d}x}}{|}_{x_1}={\dot{y}}_1\mathrm{and}{\displaystyle \frac{\mathrm{d}y}{\mathrm{d}x}}{|}_{x_2}={\dot{y}}_2.$$

(7)

Using the form given in Equation (5) with $s_1\left(x\right)=x$ and $s_2\left(x\right)=x^2,$ the constraints in Equation (7) can be expressed as

$$y(x,g\left(x\right))=g\left(x\right)+{\eta }_1{s}_1\left(x\right)+{\eta }_2{s}_2\left(x\right)=g\left(x\right)+{\eta }_{1}x+{\eta }_2{x}^2$$

(8)

where the two constants ${\eta }_1$ and ${\eta }_2$ are computed from Equations (7) and (8) as follows:

$$\left\{\begin{array}{c}{\dot{y}}_1-{\dot{g}}_1\\ {\dot{y}}_2-{\dot{g}}_2\end{array}\right\}=\left[\begin{array}{cc}{\dot{s}}_1\left(x_1\right)& {\dot{s}}_2\left(x_1\right)\\ {\dot{s}}_1\left(x_2\right)& {\dot{s}}_2\left(x_2\right)\end{array}\right]\left\{\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right\}\to \left\{\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right\}={\left[\begin{array}{cc}1& 2x_1\\ 1& 2x_2\end{array}\right]}^{-1}\left\{\begin{array}{c}{\dot{y}}_1-{\dot{g}}_1\\ {\dot{y}}_2-{\dot{g}}_2\end{array}\right\}$$

from which,

$$\left\{\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right\}={\displaystyle \frac{1}{x_2-x_1}}\left\{\begin{array}{c}{x}_2({\dot{y}}_1-{\dot{g}}_1)-x_1({\dot{y}}_2-{\dot{g}}_2)\\ -({\dot{y}}_1-{\dot{g}}_1)+({\dot{y}}_2-{\dot{g}}_2)\end{array}\right\}.$$

Substituting the ${\eta }_1$ and ${\eta }_2$ expressions into Equation (8), the following constrained expression,

$$y(x,g\left(x\right))=g\left(x\right)+\underset{{\varphi }_1(x,\mathbf{s}\left(x\right))}{\underbrace{{\displaystyle \frac{x(2x_2-x)}{2(x_2-x_1)}}}}\underset{{\rho }_1(x,g\left(x\right))}{\underbrace{({\dot{y}}_1-{\dot{g}}_1)}}+\underset{{\varphi }_2(x,\mathbf{s}\left(x\right))}{\underbrace{{\displaystyle \frac{x(x-2x_1)}{2(x_2-x_1)}}}}\underset{{\rho }_2(x,g\left(x\right))}{\underbrace{({\dot{y}}_2-{\dot{g}}_2)}},$$

(9)

is obtained for the constraints given in Equation (7). Equation (9) satisfies the constraints (7), no matter what the free function $g\left(x\right)$ is. Moreover, Equation (9) indicates the expressions of the switching functions, ${\varphi }_k(x,\mathit{s}\left(x\right))$, and the projection functionals, ${\rho }_k(x,g\left(x\right))$. Here, the meaning of the switching functions becomes more clear: when the first constraint, $\dot{y}\left(x_1\right)={\dot{y}}_1$, holds then ${\dot{\varphi }}_1\left(x_1\right)=1$ and ${\dot{\varphi }}_2\left(x_1\right)=0;$ when the second constraint, $\dot{y}\left(x_2\right)={\dot{y}}_2$, holds then ${\dot{\varphi }}_1\left(x_2\right)=0$ and ${\dot{\varphi }}_2\left(x_2\right)=1$. The projection functionals are scalars in this univariate case, but they become functionals in the multivariate case (see [9] for full explanation).

The multivariate TFC [9,16] extends the original univariate theory [15] to n dimensions and to any linear (boundary and/or internal) constraints. This extension can be summarized by the expression,

$$y(\mathit{x},g(\mathit{x}\left)\right)=h\left(c\right(\mathit{x}\left)\right)+g\left(\mathit{x}\right)-h\left(g\right(\mathit{x}\left)\right),$$

where $\mathit{x}={(x_1,x_1,\cdots ,x_n)}^{\mathrm{T}}$ is the vector of n coordinates, $c\left(\mathit{x}\right)$ is a function specifying the linear constraints, $h\left(c\right(\mathit{x}\left)\right)$ is any interpolating function satisfying the linear constraints, and $g\left(\mathit{x}\right)$ is the free function. Several examples on how to derive constrained expressions can be found in [9,17,18].

The Theory of Functional Connections has been developed for points, derivatives, integrals, infinite, component constraints, and any linear combination of them for univariate functions [15] as well as multivariate functions [16] in rectangular domains, while in generic domains some initial results have been obtained via domain mapping [19]. The main feature of these functionals (constrained expressions) is that they allow for restricting the whole function space of a constrained optimization problems to just the space of its feasible solutions, fully satisfying the constraints. This way, a large number of constrained optimization problems can be transformed into unconstrained ones, which can be solved by more simple, efficient, robust, reliable, fast, and accurate methods.

The first application of TFC was in solving linear [10] and nonlinear [11] ODEs. This has been conducted by expanding the free function $g\left(x\right)$ in terms of a set of basis functions (e.g., orthogonal polynomials, Fourier, neural networks, etc.). Linear or nonlinear least-squares method is then used to find the coefficients of the expansion. This TFC approach for solving ODEs has many advantages over traditional methods: (1) it consists of a unified framework to solve IVP, BVP, or multi-valued problems, (2) it provides an analytically approximated solution that can be used for subsequent manipulation (derivatives, integrals), (3) the solution is usually obtained in millisecond and at machine error accuracy, (4) the procedure is numerically robust (small condition number), and (5) it can solve differential equations subject to a variety of different constraint types. Additionally, TFC has been also applied to solve other mathematical optimization problems [20] such as in: homotopy continuation for control problems [21], epidemiological models [22], radiative transfer problems [23], rarefied-gas dynamics [24], Timoshenko-Ehrenfest beam [25], hybrid systems [26], machine learning [27,28,29,30], quadratic and nonlinear programming problems subject to linear equality constraints [31], orbit transfer and propagation [18,32,33], optimal control problems via indirect methods, relative motion [34], landing on small and large planetary bodies [35], and intercept problems [30].

4. Solving the Geodesic Equations Using the Theory of Functional Connections

Any trajectory on a two-dimensional surface in three-dimensional space can always be described by two coordinates, $\left[u\right(t),v(t\left)\right]$, depending on a parameter t. Specifically, let us consider a trajectory from an initial point $[u_0,v_0]$ to a final point $[u_f,v_f]$, while the parameter range is $t\in [-1,+1]$.

All possible trajectories, $\left[u\right(t),v(t\left)\right]$, connecting $[u_0,v_0]$ to $[u_f,v_f]$, can be represented by the following two functionals (constrained expressions) [9]

$$\left\{\begin{array}{c}u(t,g_u\left(t\right))=g_u\left(t\right)+{\displaystyle \frac{1-t}{2}}\left(u_0-g_{u0}\right)+{\displaystyle \frac{1+t}{2}}\left(u_f-g_{uf}\right)\hfill \\ v(t,g_v\left(t\right))=g_v\left(t\right)+{\displaystyle \frac{1-t}{2}}\left(v_0-g_{v0}\right)+{\displaystyle \frac{1+t}{2}}\left(v_f-g_{vf}\right)\hfill \end{array}\right.$$

(10)

where $g_u\left(t\right)$ and $g_v\left(t\right)$ are two free functions and where it has been set $g_{u0}=g_u(-1)$, $g_{v0}=g_v(-1)$, $g_{uf}=g_u(+1)$, and $g_{vf}=g_v(+1)$. No matter what the functions, $g_u\left(t\right)$ and $g_u\left(t\right)$ are, the Equation (10) always generate trajectories moving from $[u_0,v_0]$ to $[u_f,v_f]$, as t increases from $t=-1$ to $t=+1$.

The parametric derivatives of Equation (10) are,

$$\left\{\begin{array}{c}\dot{u}(t,g_u\left(t\right))={\dot{g}}_u\left(t\right)-{\displaystyle \frac{1}{2}}\left(u_0-g_{u0}\right)+{\displaystyle \frac{1}{2}}\left(u_f-g_{uf}\right)\hfill \\ \ddot{u}(t,g_u\left(t\right))={\ddot{g}}_u\left(t\right)\hfill \\ \dot{v}(t,g_v\left(t\right))={\dot{g}}_v\left(t\right)-{\displaystyle \frac{1}{2}}\left(v_0-g_{v0}\right)+{\displaystyle \frac{1}{2}}\left(v_f-g_{vf}\right)\hfill \\ \ddot{v}(t,g_v\left(t\right))={\ddot{g}}_v\left(t\right)\hfill \end{array}\right.$$

(11)

Let us express the free functions, $g_u\left(t\right)$ and $g_v\left(t\right)$, as a linear combination of linearly independent basis functions, $\mathit{h}\left(t\right)$, (e.g., orthogonal polynomials). Then,

$$\left\{\begin{array}{c}{g}_u\left(t\right)={\mathbf{\xi }}_u^{\mathrm{T}}\mathit{h}\left(t\right)\hfill \\ g_v\left(t\right)={\mathbf{\xi }}_v^{\mathrm{T}}\mathit{h}\left(t\right),\hfill \end{array}\right.\left\{\begin{array}{c}{\dot{g}}_u\left(t\right)={\mathbf{\xi }}_u^{\mathrm{T}}\dot{\mathit{h}}\left(t\right)\hfill \\ {\dot{g}}_v\left(t\right)={\mathbf{\xi }}_v^{\mathrm{T}}\dot{\mathit{h}}\left(t\right),\hfill \end{array}\right.\mathrm{and}\left\{\begin{array}{c}{\ddot{g}}_u\left(t\right)={\mathbf{\xi }}_u^{\mathrm{T}}\ddot{\mathit{h}}\left(t\right)\hfill \\ {\ddot{g}}_v\left(t\right)={\mathbf{\xi }}_v^{\mathrm{T}}\ddot{\mathit{h}}\left(t\right).\hfill \end{array}\right..$$

(12)

The geodesic equations provided in Equation (3) can be solved using the expressions given in Equations (10) and (11), with the free functions expanded as in Equation (12), and by discretizing the parameter t from $t=-1$ to $t=+1$. (When expressing the free functions in terms orthogonal polynomial, the best discretization of t for the least-squares problem is obtained by the Chebyshev–Gauss–Lobatto points distribution.) The resulting equations are two nonlinear algebraic equations in the unknown coefficient vectors, ${\mathbf{\xi }}_u$ and ${\mathbf{\xi }}_v$, that can be solved by nonlinear least-squares,

$${\left\{\begin{array}{c}{\mathbf{\xi }}_u\\ {\mathbf{\xi }}_v\end{array}\right\}}_{k+1}={\left\{\begin{array}{c}{\mathbf{\xi }}_u\\ {\mathbf{\xi }}_v\end{array}\right\}}_k-{\left({\mathcal{J}}_k^{\mathrm{T}}{\mathcal{J}}_k\right)}^{-1}{\mathcal{J}}_k^{\mathrm{T}}{\left\{\begin{array}{c}{\mathcal{L}}_u\\ {\mathcal{L}}_v\end{array}\right\}}_k\mathrm{where}{\mathcal{J}}_k={\left[\begin{array}{cc}{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}},& {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}\\ {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}},& {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}\end{array}\right]}_k$$

(13)

is the Jacobian of the system.

Note that, the constrained expressions given in Equation (10) are derived using constant and linear terms of support functions. This implies that the set of basis functions adopted in the least-squares process for the free functions, $g_u\left(t\right)$ and $g_v\left(t\right)$, cannot include both the constant and the linear terms. For instance, if Chebyshev orthogonal polynomials are selected as support functions, then ${\mathrm{T}}_0=1$ and ${\mathrm{T}}_1=t$ must be excluded, otherwise the matrix to invert in the least-squares process becomes singular. This is because a least-squares process of a linear combination of functions is singular if not all the functions are linearly independent.

Initial Guess and the Local Minima Problem

Since the constrained expressions provide trajectories always satisfying the constraints (for any expression of the free functions), then the most natural (and simplest) initial guess is to start the iterative least-squares solving Equation (13) by setting ${\mathbf{\xi }}_u={\mathbf{\xi }}_v=\mathbf{0}$. This is equivalent to initially select $g_u\left(t\right)=g_v\left(t\right)=0$ and, consequently—see Equation (10)—, selecting an initial linear variation of the parametric variables, u and v, from the initials ($u_0,v_0$) to the final values ($u_f,v_f$).

Nonlinear least-squares applied to find the geodesic trajectory is, unfortunately, a process affected by local minima. Typically, when an iterative procedure enters into the convergence phase, then each following step is smaller than the previous,

$$L_2\left(\Delta {\mathbf{\xi }}_{k+1}\right)<L_2\left(\Delta {\mathbf{\xi }}_k\right)\mathrm{where}\Delta {\mathbf{\xi }}_k={\left\{\begin{array}{c}\Delta {\mathbf{\xi }}_u\\ \Delta {\mathbf{\xi }}_v\end{array}\right\}}_k$$

(14)

This convergence criteria has two interesting aspects: (1) it can be used to push the convergence to the maximum accuracy and (2) no tolerance is needed to stop the iterations. In fact, when the convergence reaches the maximum accuracy, the procedure cannot anymore improve the estimation of $\Delta {\mathbf{\xi }}_k$ and the inequality, $L_2\left(\Delta {\mathbf{\xi }}_{k+1}\right)>L_2\left(\Delta {\mathbf{\xi }}_k\right)$, happens because of numerical errors (convergence saturation).

Figure 1 shows the flowchart of the algorithm adopted to avoid local minima. Starting with ${\mathbf{\xi }}_0=\mathbf{0}$, the least-squares iterations continue until Equation (14) is verified for N consecutive times (initial convergence loop). During this sequence, if the least-squares diverges, then the algorithm restarts with a new random initial guess. If Equation (14) is verified for N consecutive times, then the algorithm enters into the convergence loop and exits when $L_2\left(\Delta {\mathbf{\xi }}_{k+1}\right)>L_2\left(\Delta {\mathbf{\xi }}_k\right)$ is experienced. Then, the $L_2$ norm of the residuals is computed. The $L_2$ norm of the residuals allows one to discriminate local minima and global minima (associated with the geodesic trajectory) using a small tolerance, such as $\epsilon ={10}^{-15}$. If $L_2\left({\mathit{r}}_k\right)>\epsilon$ then the algorithm restarts with a new random initial guess.

5. Numerical Validations

In this section, some boundary value geodesic problems are numerically solved to validate the proposed methodology on 2-dimensional surfaces, for visualization purposes. The proposed approach has been tested on eight different kind of surfaces: triaxial ellipsoid, elliptic hyperboloid, elliptic paraboloid, hyperbolic paraboloid, one-sheeted hyperboloid, torus, Moëbius strip, and a generic surface. Six of them are shown in Figure 2.

5.1. Triaxial Ellipsoid

A (non-oriented) triaxial ellipsoid is described by,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}asinucosv,& bsinusinv,& ccosu\end{array}\right\}$$

where $u\in [0,\pi ]$ and $v\in [0,2\pi )$, and $(a,b,c)$ are the three ellipsoid semi-major axes. The covariant metric tensor for the ellipsoid is,

$$g_{ij}=\left[\begin{array}{cc}{cos}^{2}u(a^2{cos}^{2}v+b^2{sin}^{2}v)+c^2{sin}^{2}u& -{\displaystyle \frac{sin\left(2u\right)sin\left(2v\right)}{4}}(a^2-b^2)\\ -{\displaystyle \frac{sin\left(2u\right)sin\left(2v\right)}{4}}(a^2-b^2)& {sin}^{2}u\left[{sin}^{2}v(a^2-b^2)+b^2\right]\end{array}\right]$$

Since the ellipsoid is an affine image of the unit sphere, it is affected by the same singularity problem, occurring for $u=0$ and $u=\pi$.

By setting $a=2$, $b=3$, and $c=1$, the geodesic equations for this ellipsoid are [36],

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\ddot{u}+{\displaystyle \frac{(5{cos}^{2}v-32)cosusinu}{{sin}^{2}u\left(5{cos}^{2}v-32\right)+36}}{\dot{u}}^2-{\displaystyle \frac{36cosusinu}{{sin}^{2}u\left(5{cos}^{2}v-32\right)+36}}{\dot{v}}^2=0\hfill \\ \ddot{v}-{\displaystyle \frac{5cosvsinv}{{sin}^{2}u\left(5{cos}^{2}v-32\right)+36}}{\dot{u}}^2+{\displaystyle \frac{2cosu}{sinu}}\dot{u}\dot{v}-{\displaystyle \frac{5{sin}^{2}ucosvsinv}{{sin}^{2}u\left(5{cos}^{2}v-32\right)+36}}{\dot{v}}^2=0\hfill \end{array}\right.\end{array}$$

To solve these differential equations using nonlinear least-squares, the following partials must be evaluated to compute the Jacobian of the system,

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}}={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}·{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}·{\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}·{\displaystyle \frac{\mathrm{d}\ddot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial v}}·{\displaystyle \frac{\mathrm{d}v}{\mathrm{d}{\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}·{\displaystyle \frac{\mathrm{d}\dot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}}={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}·{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}·{\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial v}}·{\displaystyle \frac{\mathrm{d}v}{\mathrm{d}{\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}·{\displaystyle \frac{\mathrm{d}\dot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}·{\displaystyle \frac{\mathrm{d}\ddot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}\hfill \end{array}\right.$$

(15)

Figure 3 shows the numerical results using 40 basis functions (Legendre orthogonal polynomials) to describe the free functions and 200 discretization points. This test validates the approach in terms of fast convergence (just six iterations using a noisy initial guess) and in terms of finding constant the parametric velocity.

The geodetic velocity on a generic triaxial ellipsoid has components,

$$\begin{array}{c}\hfill \dot{\mathit{p}}=\left\{\begin{array}{c}a(\dot{u}cosucosv-\dot{v}sinusinv)\\ b(\dot{u}cosusinv+\dot{v}sinucosv)\\ -c\dot{u}sinu\end{array}\right\}\end{array}$$

5.2. Elliptic Paraboloid

For the elliptic paraboloid, the geodesic equations and the TFC solution procedure are provided. The parametric equations of the elliptic paraboloid, $z=x^2+y^2$, can be described by the vector,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}u,& v,& u^2+v^2\end{array}\right\}$$

where $u\in [u_{min},u_{max}]$ and $v\in [v_{min},v_{max}]$. The partials of the generic point, $\mathit{p}$, are

$${\displaystyle \frac{\partial \mathit{p}}{\partial u}}={\mathit{p}}_u=\left\{\begin{array}{c}1\\ 0\\ 2u\end{array}\right\}\mathrm{and}{\displaystyle \frac{\partial \mathit{p}}{\partial v}}={\mathit{p}}_v=\left\{\begin{array}{c}0\\ 1\\ 2v\end{array}\right\}$$

from which the metric tensor,

$$g_{ij}=\left[\begin{array}{cc}{\mathit{p}}_u^{\mathrm{T}}{\mathit{p}}_u& {\mathit{p}}_u^{\mathrm{T}}{\mathit{p}}_v\\ {\mathit{p}}_v^{\mathrm{T}}{\mathit{p}}_u& {\mathit{p}}_v^{\mathrm{T}}{\mathit{p}}_v\end{array}\right]=\left[\begin{array}{cc}1+4u^2& 4uv\\ 4uv& 1+4v^2\end{array}\right]$$

and its inverse,

$$g^{ij}={\left(g_{ij}\right)}^{-1}={\displaystyle \frac{1}{4(u^2+v^2)+1}}\left[\begin{array}{cc}1+4v^2& -4uv\\ -4uv& 1+4u^2\end{array}\right]$$

are derived. The non-null Christoffel symbols are derived using Equation (4),

$$\begin{array}{c}{\mathsf{\Gamma }}_{11}^1={\displaystyle \frac{4u}{4(u^2+v^2)+1}},{\mathsf{\Gamma }}_{22}^1={\displaystyle \frac{4u}{4(u^2+v^2)+1}},\\ {\mathsf{\Gamma }}_{11}^2={\displaystyle \frac{4v}{4(u^2+v^2)+1}},\mathrm{and}{\mathsf{\Gamma }}_{22}^2={\displaystyle \frac{4v}{4(u^2+v^2)+1}}.\end{array}$$

The geodesic equations, given in Equation (3), become

$$\left\{\begin{array}{c}\ddot{u}+{\mathsf{\Gamma }}_{11}^1{\dot{u}}^2+{\mathsf{\Gamma }}_{22}^1{\dot{v}}^2=0\hfill \\ \ddot{v}+{\mathsf{\Gamma }}_{11}^2{\dot{u}}^2+{\mathsf{\Gamma }}_{22}^2{\dot{v}}^2=0\hfill \end{array}\right.$$

After substituting the Christoffel symbols and rearranging the equations, the geodesic equations for the elliptic paraboloid become,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=(4u^2+4v^2+1)\ddot{u}+4u{\dot{u}}^2+4u{\dot{v}}^2=0\hfill \\ {\mathcal{L}}_v=(4u^2+4v^2+1)\ddot{v}+4v{\dot{u}}^2+4v{\dot{v}}^2=0\hfill \end{array}\right.$$

(16)

To solve these differential equations using nonlinear least-squares, the same structure of partials provided in Equation (15) must be evaluated to compute the Jacobian of the system, where

$$\begin{array}{ccc}{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}=8u\ddot{u}+4({\dot{u}}^2+{\dot{v}}^2),\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}=8u\dot{u},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}=4(u^2+v^2)+1,\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}=8u\dot{v},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}=8u\ddot{v},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}=8v\dot{u},\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial v}}=8v\ddot{v}+4({\dot{u}}^2+{\dot{v}}^2),\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}=8v\dot{v},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial v}}=8v\ddot{u},\hfill \end{array}$$

Note that, the elliptic paraboloid can also be expressed by, $\mathit{p}={\{usinv,ucosv,u^2\}}^{\mathrm{T}}$, where u and v are angles. In this case, the geodesic equations provided in Equation (16) can be rewritten as,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=\left(4u^2+1\right)\ddot{u}+4u{\dot{u}}^2-u{\dot{v}}^2=0\hfill \\ {\mathcal{L}}_v=u\ddot{v}+2\dot{u}\dot{v}=0\hfill \end{array}\right.$$

and the following partials must be computed to solve the boundary value problem,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}·{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}·{\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}·{\displaystyle \frac{\mathrm{d}\ddot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}& ={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}·{\displaystyle \frac{\mathrm{d}\dot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}·{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}{\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}·{\displaystyle \frac{\mathrm{d}\dot{u}}{\mathrm{d}{\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}& ={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}·{\displaystyle \frac{\mathrm{d}\dot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}·{\displaystyle \frac{\mathrm{d}\ddot{v}}{\mathrm{d}{\mathbf{\xi }}_v}}\hfill \end{array}$$

where,

$$\begin{array}{cccc}{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}=8u\ddot{u}+4{\dot{u}}^2-{\dot{v}}^2,\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}=8u\dot{u},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}=4u^2+1,\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}=-2u\dot{v},\hfill \\ {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}=\ddot{v},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}=2\dot{v},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}=2\dot{u},\hfill & {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}=u,\hfill \end{array}$$

and,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}=\mathbf{h}-{\displaystyle \frac{1-t}{2}}{\mathit{h}}_0-{\displaystyle \frac{t+1}{2}}{\mathit{h}}_f,\hfill \\ \hfill {\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}=\dot{\mathit{h}}+{\displaystyle \frac{1}{2}}{\mathbf{h}}_0-{\displaystyle \frac{1}{2}}{\mathit{h}}_f,\hfill \\ \hfill {\displaystyle \frac{\partial \ddot{u}}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial \ddot{v}}{\partial {\mathbf{\xi }}_v}}=\ddot{\mathit{h}},\hfill \end{array}$$

This second parametrization of the elliptic paraboloid is provided because a particular attention must be given to the boundary value of angles to avoid discontinuous angle evolution. For example, if the absolute difference between the two angles is $|v_0-v_f|>\pi$, then the value of the smallest of these two angles is increased by $2\pi$. This avoids the discontinuity at $v=0$ while unchanging the problem. For instance, if $v_0=3\pi /2$ and $v_f=\pi /4$, then the value of $v_f$ is set as $v_f=\pi /4+2\pi =9\pi /4$.

Geodesic (Parametric) Velocity on an Elliptic Paraboloid

The derivative of the position vector is,

$${\dot{\mathit{p}}}^{\mathrm{T}}=\left\{\begin{array}{ccc}{\dot{p}}_x,& {\dot{p}}_y,& {\dot{p}}_z\end{array}\right\}=\left\{\begin{array}{ccc}\dot{u},& \dot{v},& 2u\dot{u}+2v\dot{v}\end{array}\right\}$$

Therefore, the velocity is provided by,

$$\dot{p}\left(t\right)=\sqrt{{\dot{u}}^2+u^2+4{\left(u\dot{u}+v\dot{v}\right)}^2}$$

Since the velocity on a geodesic is constant, then the (instantaneous) value $\dot{p}\left(t\right)$ for a geodesic is not a function of t and, therefore, the expression of $\dot{p}\left(t\right)$ must coincide with the average value of the velocity which, in turn, can be computed using Equation (2),

$${\displaystyle \overline{\dot{p}}={\displaystyle \frac{L}{t_f-t_0}}={\displaystyle \frac{1}{t_f-t_0}}{\int }_{t_0}^{t_f}\sqrt{\left(v^2+1\right){\dot{u}}^2+2uv\dot{u}\dot{v}+\left(u^2+1\right){\dot{v}}^2}\mathrm{d}t}$$

therefore, for a geodesic trajectory on the elliptic paraboloid, we have a closed form expression for the integral,

$${\displaystyle {\int }_{t_0}^{t_f}\sqrt{\left(v^2+1\right){\dot{u}}^2+2uv\dot{u}\dot{v}+\left(u^2+1\right){\dot{v}}^2}\mathrm{d}t=\left(t_f-t_0\right)\sqrt{{\dot{u}}^2+u^2+4(u\dot{u}+v\dot{v})}}$$

5.3. Elliptic Hyperboloid

The parametric description of the elliptic hyperboloid is,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}asinhucosv,& bsinhusinv,& ccoshu\end{array}\right\}$$

and his covariant metric tensor is,

$$g_{ij}=\left[\begin{array}{cc}{cosh}^{2}u(a^2{cos}^{2}v+b^2{sin}^{2}v)+c^2{sinh}^{2}u& -sinh\left(2u\right)sin\left(2v\right){\displaystyle \frac{a^2-b^2}{4}}\\ -sinh\left(2u\right)sin\left(2v\right){\displaystyle \frac{a^2-b^2}{4}}& {sinh}^{2}u((a^2-b^2){sin}^{2}v+b^2)\end{array}\right]$$

By setting $a=1$, $b=2$, and $c=3$, the geodesic equations are,

$$\begin{array}{cc}\hfill \ddot{u}& +{\displaystyle \frac{coshusinhu\left(13+27{cos}^{2}v\right)}{27{cosh}^{2}u{cos}^{2}v+13{cosh}^{2}u-9-27{cos}^{2}v}}{\dot{u}}^2+\hfill \\ \hfill & +{\displaystyle \frac{4coshusinhu}{27{cosh}^{2}u{cos}^{2}v+13{cosh}^{2}u-9-27{cos}^{2}v}}{\dot{v}}^2=0\hfill \\ \hfill \ddot{v}& -{\displaystyle \frac{27cosvsinv}{27{cosh}^{2}u{cos}^{2}v+13{cosh}^{2}u-9-27{cos}^{2}v}}{\dot{u}}^2+\hfill \\ \hfill & +{\displaystyle \frac{2coshu}{sinhu}}\dot{u}\dot{v}-{\displaystyle \frac{27cosvsinv\left({cosh}^{2}u-1\right)}{27{cosh}^{2}u{cos}^{2}v+13{cosh}^{2}u-9-27{cos}^{2}v}}{\dot{v}}^2=0\hfill \end{array}$$

The geodesic velocity components on the elliptic hyperboloid are,

$$\left\{\begin{array}{c}{\dot{p}}_x=a(\dot{u}coshucosv-\dot{v}sinhusinv)\hfill \\ {\dot{p}}_y=b(\dot{u}coshusinv+\dot{v}sinhucosv)\hfill \\ {\dot{p}}_z=c\dot{u}sinhu\hfill \end{array}\right.$$

Figure 4 provides detailed information of the test results using 40 basis functions (for the free functions) and 300 points discretization. In this case, a very noisy initial guess is selected (black trajectory on the left/top figure), instead of the simple, ${\mathbf{\xi }}_0=\mathbf{0}$. The right/top figure shows, for subsequent iterations, the $L_2$ norm of the error. The procedure entered into the convergence phase after about 35 iterations. The convergence loop pushed the accuracy to the limit, when $L_2\left(\Delta {\mathbf{\xi }}_{k+1}\right)>L_2\left(\Delta {\mathbf{\xi }}_k\right)$ occurred. The final solution residuals are shown in the left/bottom figure, while the central/bottom figure shows the initial and final solutions for $u\left(t\right)$ and $v\left(t\right)$. Finally, the right/bottom figure gives the parametric velocity of the initial guess (black) and of the final (red) solution (geodesic trajectory). This validates the parametric velocity with a constant value. For this example, no local minima were ever experienced (convex problem).

5.4. Hyperbolic Paraboloid

The parametric description of the hyperbolic paraboloid is,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}u,& v,& uv\end{array}\right\}$$

the covariant metric tensor for the hyperbolic paraboloid is

$$g_{ij}=\left[\begin{array}{cc}{v}^2+1& uv\\ uv& u^2+1\end{array}\right]$$

and the Christoffel symbols are, ${\mathsf{\Gamma }}_{11}^2={\displaystyle \frac{v}{u^2+v^2+1}}$, ${\mathsf{\Gamma }}_{12}^1={\displaystyle \frac{v}{u^2+v^2+1}}$, ${\mathsf{\Gamma }}_{21}^2={\displaystyle \frac{u}{u^2+v^2+1}}$, and ${\mathsf{\Gamma }}_{22}^1={\displaystyle \frac{u}{u^2+v^2+1}}$, while the geodesic equations are,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=\left(1+u^2+v^2\right)\ddot{u}+2v\dot{u}\dot{v}=0\hfill \\ {\mathcal{L}}_v=\left(1+u^2+v^2\right)\ddot{v}+2u\dot{u}\dot{v}=0\hfill \end{array}\right.$$

The expression of the geodesic velocity on hyperbolic paraboloid is,

$${\dot{p}}^2={\dot{u}}^2+{\dot{v}}^2+{(v\dot{u}+u\dot{v})}^2.$$

Figure 5 contains the results obtained on hyperbolic paraboloid (same meaning than those provided for the elliptic paraboloid), while

Table 1. First eight iterations for the hyperbolic paraboloid.

Iteration k	${\mathit{L}}_2\left(\mathbf{\Delta }{\mathit{\xi }}_{\mathit{k}}\right)$	${\mathit{L}}_2\left({\mathit{r}}_{\mathit{k}}\right)$
0	-	78,019.3597
1	$0.014083$	$682.5215$
2	$0.13725$	$51.5259$
3	$0.10102$	$4.9834$
4	$0.029109$	$0.10211$
5	$0.00097301$	$7.6775\times {10}^{-5}$
6	$9.8617\times {10}^{-7}$	$7.4134\times {10}^{-11}$
7	$1.0705\times {10}^{-12}$	$7.3451\times {10}^{-17}$
8	$2.3949\times {10}^{-19}$	$8.1203\times {10}^{-17}$

shows the $L_2$ norm of the error step, $L_2\left(\Delta {\mathbf{\xi }}_k\right)$, and the associated $L_2$ norm of the differential equations residuals, $L_2\left({\mathbf{r}}_k\right)$, for the first 8 iterations.

5.5. Generic Surface

Let us consider the surface identified by,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}u,& v,& ucosv-vsinu\end{array}\right\}$$

in the range $u,v\in [-5,+5]$. By setting, $\alpha =sinu+usinv$ and $\beta =cosv-vcosu$, the metric tensor and its inverse for this surface are,

$$g_{ij}=\left[\begin{array}{cc}{\beta }^2+1& -\alpha \beta \\ -\alpha \beta & {\alpha }^2+1\end{array}\right]\mathrm{and}{g}^{ij}={\displaystyle \frac{1}{{\alpha }^2+{\beta }^2+1}}\left[\begin{array}{cc}{\alpha }^2+1& \alpha \beta \\ \alpha \beta & {\beta }^2+1\end{array}\right]$$

and the nonzero Christoffel symbols are,

$$\begin{array}{ccc}{\mathsf{\Gamma }}_{11}^1={\displaystyle \frac{vsinu}{{\alpha }^2+{\beta }^2+1}}\beta ,\hfill & {\mathsf{\Gamma }}_{12}^1={\mathsf{\Gamma }}_{21}^1=-{\displaystyle \frac{cosu+sinv}{{\alpha }^2+{\beta }^2+1}}\beta ,\hfill & {\mathsf{\Gamma }}_{22}^1=-{\displaystyle \frac{ucosv}{{\alpha }^2+{\beta }^2+1}}\beta ,\hfill \\ {\mathsf{\Gamma }}_{11}^2=-{\displaystyle \frac{vsinu}{{\alpha }^2+{\beta }^2+1}}\alpha ,\hfill & {\mathsf{\Gamma }}_{12}^2={\mathsf{\Gamma }}_{21}^2={\displaystyle \frac{cosu+sinv}{{\alpha }^2+{\beta }^2+1}}\alpha ,\hfill & {\mathsf{\Gamma }}_{22}^2={\displaystyle \frac{ucosv}{{\alpha }^2+{\beta }^2+1}}\alpha \hfill \end{array}$$

The geodesic equations can be written in the compact form,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=p\ddot{u}+\beta q=0\hfill \\ {\mathcal{L}}_v=p\ddot{v}-\alpha q=0\hfill \end{array}\right.$$

where,

$$p={\alpha }^2+{\beta }^2+1\mathrm{and}q=(vsinu){\dot{u}}^2-2(cosu+sinv)\dot{u}\dot{v}-(ucosv){\dot{v}}^2$$

The partials forming the Jacobian matrix in the least-squares process are,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}}& =\ddot{u}{\displaystyle \frac{\partial p}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}+p{\displaystyle \frac{\partial \ddot{u}}{\partial {\mathbf{\xi }}_u}}+\beta \left({\displaystyle \frac{\partial q}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial q}{\partial \dot{u}}}·{\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_x}}\right)+q{\displaystyle \frac{\partial \beta }{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}& =\ddot{u}{\displaystyle \frac{\partial p}{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}+\beta \left({\displaystyle \frac{\partial q}{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial q}{\partial \dot{v}}}·{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}\right)+q{\displaystyle \frac{\partial \beta }{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}}& =\ddot{v}{\displaystyle \frac{\partial p}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}-\alpha \left({\displaystyle \frac{\partial q}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial q}{\partial \dot{u}}}·{\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_u}}\right)-q{\displaystyle \frac{\partial \alpha }{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}& =\ddot{v}{\displaystyle \frac{\partial p}{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}+p{\displaystyle \frac{\partial \ddot{v}}{\partial {\mathbf{\xi }}_v}}-\alpha \left({\displaystyle \frac{\partial q}{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial q}{\partial \dot{v}}}·{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}\right)-q{\displaystyle \frac{\partial \alpha }{\partial v}}·{\displaystyle \frac{\partial v}{\partial {\mathbf{\xi }}_v}}\hfill \end{array}$$

where,

$${\displaystyle \frac{\partial \alpha }{\partial u}}=cosu+sinv,{\displaystyle \frac{\partial \alpha }{\partial v}}=ucosv,{\displaystyle \frac{\partial \beta }{\partial u}}=vsinu,\mathrm{and}{\displaystyle \frac{\partial \beta }{\partial v}}=-sinv-cosu$$

and,

$${\displaystyle \frac{\partial p}{\partial u}}=2\alpha {\displaystyle \frac{\partial \alpha }{\partial u}}+2\beta {\displaystyle \frac{\partial \beta }{\partial u}}\mathrm{and}{\displaystyle \frac{\partial p}{\partial v}}=2\beta {\displaystyle \frac{\partial \beta }{\partial v}}+2\alpha {\displaystyle \frac{\partial \alpha }{\partial v}}$$

and,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial q}{\partial u}}& =(vcosu){\dot{u}}^2+2(sinu)\dot{u}\dot{v}-(cosv){\dot{v}}^2\hfill \\ \hfill {\displaystyle \frac{\partial q}{\partial v}}& =(sinu){\dot{u}}^2-2(cosv)\dot{u}\dot{v}+(usinv){\dot{v}}^2\hfill \\ \hfill {\displaystyle \frac{\partial q}{\partial \dot{u}}}& =2(vsinu)\dot{u}-2(cosu+sinv)\dot{v}\hfill \\ \hfill {\displaystyle \frac{\partial q}{\partial \dot{v}}}& =-2(cosu+sinv)\dot{u}-2(ucosv)\dot{v}\hfill \end{array}$$

Figure 6 gives the results obtained for this specific generic surface. The convergence is obtained in just 8 iterations, using 30 basis functions for the free functions and a 300 point discretization. The constancy of the parametric velocity is an alternative way validating the solution as a geodesic trajectory.

5.6. One-Sheeted Hyperboloid

This surface is described by the parametric vector,

$${\mathit{p}}^{\mathrm{T}}(u,v)=\left\{\begin{array}{ccc}\sqrt{1+u^2}cosv,& \sqrt{1+u^2}sinv,& 2u\end{array}\right\}$$

The covariant and the contravariant associated metric tensors are,

$$g_{ij}=\left[\begin{array}{cc}{\displaystyle \frac{5u^2+4}{u^2+1}}& 0\\ 0& u^2+1\end{array}\right]\mathrm{and}{g}^{ij}=\left[\begin{array}{cc}{\displaystyle \frac{u^2+1}{5u^2+4}}& 0\\ 0& {\displaystyle \frac{1}{u^2+1}}\end{array}\right]$$

and the nonzero Christoffel symbols are,

$${\mathsf{\Gamma }}_{11}^1={\displaystyle \frac{u}{(5u^2+4)(u^2+1)}},{\mathsf{\Gamma }}_{22}^1=-{\displaystyle \frac{u\left(u^2+1\right)}{5u^2+4}},\mathrm{and}{\mathsf{\Gamma }}_{12}^2={\mathsf{\Gamma }}_{21}^2={\displaystyle \frac{u}{u^2+1}}$$

This implies the geodesic equations for one-sheeted hyperboloid,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=\left(5u^2+4\right)\ddot{u}+{\displaystyle \frac{u}{u^2+1}}{\dot{u}}^2-u\left(u^2+1\right){\dot{v}}^2=0\hfill \\ {\mathcal{L}}_v=\left(u^2+1\right)\ddot{v}+2u\dot{u}\dot{v}=0\hfill \end{array}\right.$$

To apply the nonlinear least-squares, the following partials must be computed,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}·{\displaystyle \frac{\partial \ddot{u}}{\partial {\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}·{\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}& ={\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}·{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}}& ={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}·{\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_u}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}·{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}& ={\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}·{\displaystyle \frac{\partial \ddot{v}}{\partial {\mathbf{\xi }}_v}}+{\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}·{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}\hfill \end{array}$$

The expressions of these partials are,

$$\begin{array}{cccc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \ddot{u}}}& =5u^2+4\hfill & \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{u}}}& ={\displaystyle \frac{2u}{u^2+1}}\dot{u}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial u}}& =10u^2\ddot{u}+{\displaystyle \frac{1-u^2}{{(u^2+1)}^2}}{\dot{u}}^2-\left(3u^2+1\right){\dot{v}}^2\hfill & \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial \dot{v}}}& =-2u\left(u^2+1\right)\dot{v}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \ddot{v}}}& =u^2+1\hfill & \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{v}}}& =2u\dot{u}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial \dot{u}}}& =2u\dot{v}\hfill & \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial u}}& =2u\ddot{v}+2\dot{u}\dot{v}\hfill \end{array}$$

The parametric velocity is given by

$${\dot{p}}^2={\displaystyle \frac{u^2{\dot{u}}^2}{1+u^2}}+4{\dot{u}}^2+(1+u^2)\dot{v}{.}^2$$

The results of the test conducted on the one-sheeted hyperboloid are reported in Figure 7 where the sub-figures provide full information of the least-squares process. The free functions were expanded by 40 basis functions (Legendre orthogonal polynomials) and the discretization was conducted by 300 points.

5.7. Torus

The implicit and parametric equation of the torus are,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}(R+rcosu)cosv,& (R+rcosu)sinv,& rsinu\end{array}\right\}$$

where R and r are the two torus radii. The covariant metric tensor for the torus is,

$$g_{ij}=\left[\begin{array}{cc}{r}^2& 0\\ 0& {(R+rcosu)}^2\end{array}\right]$$

while the geodesic equations are,

$$\left\{\begin{array}{c}{\mathcal{L}}_u=r\ddot{u}+{\dot{v}}^2(R+rcosu)sinu=0\hfill \\ {\mathcal{L}}_v=(R+rcosu)\ddot{v}-2\dot{u}\dot{v}rsinu=0\hfill \end{array}\right.$$

The partials needed to populate the Jacobian are,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_u}}& =r{\displaystyle \frac{\partial \ddot{u}}{\partial {\mathbf{\xi }}_u}}+{\dot{v}}^2\left[(R+rcosu)cosu-r{sin}^{2}u\right]{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_u}{\partial {\mathbf{\xi }}_v}}& =2\dot{v}{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}(R+rcosu)sinu\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_u}}& =-\left[r\ddot{v}sinu+2\dot{u}\dot{v}rcosu\right]{\displaystyle \frac{\partial u}{\partial {\mathbf{\xi }}_u}}-2{\displaystyle \frac{\partial \dot{u}}{\partial {\mathbf{\xi }}_u}}\dot{v}rsinu\hfill \\ \hfill {\displaystyle \frac{\partial {\mathcal{L}}_v}{\partial {\mathbf{\xi }}_v}}& =(R+rcosu){\displaystyle \frac{\partial \ddot{v}}{\partial {\mathbf{\xi }}_v}}-2\dot{u}{\displaystyle \frac{\partial \dot{v}}{\partial {\mathbf{\xi }}_v}}rsinu\hfill \end{array}$$

and the nonlinear least-squares can be performed using Equation (13). The geodetic parametric velocity on the torus is,

$${\dot{p}}^2=r^2{\dot{u}}^2+{(R+rcosu)}^2{\dot{v}}^2$$

The sub-figures given in Figure 8 show the details of the least-squares process to estimate the geodesic trajectory on the torus surface. The number of basis functions were 40 and the number of discretization points were 300. The constancy of the velocity as well as the machine-level $L_2$ norm of the differential equations residuals validate the estimated geodesic trajectory (in red).

5.8. Moëbius Strip

The last test of the proposed least-squares approach is performed on the Moëbius strip, which can be described by,

$${\mathit{p}}^{\mathrm{T}}=\left\{\begin{array}{ccc}\left[2-usin\left({\displaystyle \frac{v}{2}}\right)\right]cos\left(v\right),& \left[2-usin\left({\displaystyle \frac{v}{2}}\right)\right]sin\left(v\right),& ucos\left({\displaystyle \frac{v}{2}}\right)\end{array}\right\}$$

where $u\in [-1,+1]$ and $v\in [0,2\pi )$. The metric tensor and its inverse for this surface can be given in a compact form by setting, $d(u,v)=4usin\left({\displaystyle \frac{v}{2}}\right)\left[usin\left({\displaystyle \frac{v}{2}}\right)-4\right]+u^2+16$,

$$g_{ij}=\left[\begin{array}{cc}1& 0\\ 0& {\displaystyle \frac{d(u,v)}{4}}\end{array}\right]\mathrm{and}{g}^{ij}=\left[\begin{array}{cc}1& 0\\ 0& {\displaystyle \frac{4}{d(u,v)}}\end{array}\right]$$

while the nonzero Christoffel symbols are,

$$\begin{array}{cc}\hfill {\mathsf{\Gamma }}_{22}^1& =sin\left({\displaystyle \frac{v}{2}}\right)\left[2-usin\left({\displaystyle \frac{v}{2}}\right)\right]-{\displaystyle \frac{u}{4}}\hfill \\ \hfill {\mathsf{\Gamma }}_{21}^2& ={\mathsf{\Gamma }}_{12}^2={\displaystyle \frac{u-4sin\left({\displaystyle \frac{v}{2}}\right)\left[2-usin\left({\displaystyle \frac{v}{2}}\right)\right]}{d(u,v)}}\hfill \\ \hfill {\mathsf{\Gamma }}_{22}^2& =-{\displaystyle \frac{2ucos\left({\displaystyle \frac{v}{2}}\right)\left[2-usin\left({\displaystyle \frac{v}{2}}\right)\right]}{d(u,v)}}\hfill \end{array}$$

Hence, the geodesic equations of the Moëbius strip are,

$$\ddot{u}+{\mathsf{\Gamma }}_{22}^1{\dot{v}}^2=0\mathrm{and}\ddot{v}+2{\mathsf{\Gamma }}_{21}^2\dot{u}\dot{v}+{\mathsf{\Gamma }}_{22}^2{\dot{v}}^2=0$$

which can be written as,

$$\left\{\begin{array}{c}\ddot{u}+c_1(u,v){\dot{v}}^2=0\hfill \\ d(u,v)\ddot{v}+c_2(u,v)\dot{u}\dot{v}+c_3(u,v){\dot{v}}^2=0\hfill \end{array}\right.$$

where

$$\left\{\begin{array}{c}{c}_1(u,v)=2sin\left({\displaystyle \frac{v}{2}}\right)-u{sin}^2\left({\displaystyle \frac{v}{2}}\right)-{\displaystyle \frac{u}{4}}\hfill \\ c_2(u,v)=2u-16sin\left({\displaystyle \frac{v}{2}}\right)+8u{sin}^2\left({\displaystyle \frac{v}{2}}\right)\hfill \\ c_3(u,v)=-4ucos\left({\displaystyle \frac{v}{2}}\right)+4u^{2}sinv\hfill \end{array}\right.$$

Figure 9 provide all information about this last test validation case.

Table 2. Moëbius strip: convergence $L_2$ norms terms used in the least-squares approach.

Iteration k	${\mathit{L}}_2\left(\mathbf{\Delta }{\mathit{\xi }}_{\mathit{k}}\right)$	${\mathit{L}}_2\left({\mathit{r}}_{\mathit{k}}\right)$
1	-	$4.4143$
2	2,590,717.6067	$2.7918$
3	125,575.9139	$1.4362$
4	$324.6177$	$0.7541$
5	27,398.8611	$0.87001$
6	$2583.4652$	$0.60306$
7	$2142.7967$	$0.39897$
8	$2330.2236$	$1.2032$
9	$1485.2018$	$0.55737$
10	$1197.0006$	$0.61434$
11	$311.3524$	$0.14478$
12	$487.9178$	$0.012141$
13	$35.4515$	$0.0008965$
14	$1.1049$	$2.1828\times {10}^{-6}$
15	$0.0088866$	$1.433\times {10}^{-7}$
16	$0.00014298$	$3.7817\times {10}^{-9}$
17	$7.4714\times {10}^{-6}$	$1.5985\times {10}^{-10}$
18	$1.5249\times {10}^{-7}$	$6.5714\times {10}^{-12}$
19	$3.3531\times {10}^{-10}$	$2.6535\times {10}^{-13}$
20	$5.5449\times {10}^{-13}$	$1.0555\times {10}^{-14}$
21	$4.0001\times {10}^{-15}$	$4.1414\times {10}^{-16}$
22	$2.8094\times {10}^{-15}$	$4.9073\times {10}^{-17}$
23	$5.1466\times {10}^{-15}$	$3.4713\times {10}^{-17}$

details the iterative results: the $L_2$ norm of the accuracy gain and the $L_2$ norm of the differential equation at each iteration. The iterative process ended because $L_2\left(\Delta {\mathbf{\xi }}_{23}\right)>L_2\left(\Delta {\mathbf{\xi }}_{22}\right)$, due to numerical convergence saturation.

5.9. Discussions

In this article, a new general method to numerical solve boundary value geodesic problems in curved surfaces is presented. The approach takes advantage of the ability to derive special functionals, called constrained expressions, which always satisfy assigned boundary conditions. The proposed approach has been validated by performing numerical tests on several two-dimensional surfaces. In theory, this approach can be extended to manifolds in higher Riemannian spaces. This will require more computational capability than that used for this article and, most likely, optimized and compiled code, instead of the MATLAB interpreter adopted.

The problem of finding the geodesic by least-squares introduces the risk of getting stuck on some local minima. A simple algorithm has been developed to mitigate this problem. Thanks to this algorithm, all boundary geodesic problems considered were quickly solved with machine-error level accuracy, which is an acceptable definition of an exact solution by engineers.

This article restricts the research on finding geodesic trajectories on surfaces that are continuous and differentiable. The problem of finding geodesic-type trajectories on discretized surfaces, which is solved by combinatorial/computational algorithms and methods, is not taken here into consideration. The performed tests have the only purpose to validate the least-squares approach. A complete analysis quantifying the responses of this methodology to various surface shapes, boundary conditions (e.g., singular boundary points for the Ellipsoid), as well as comparisons with competing approaches, will be the subject of future studies.

As for future research directions, it should be natural to derive the real velocity (instead of the parametric velocity) in geodesic problems of a mass particle, and to use the constancy of the velocity—in addition to or instead of—the geodesic equations. Another future research activity is to investigate what the optimal range of basis functions for the free functions is and to investigate the optimal number of discretization points. These two analyses will help the proposed least-squares solution approach to provide optimal performances. Topics more interesting in physics, such as using the Schwarzschild metric (uncharged, non-rotating black holes) or more general space metrics are welcome to be investigated by researchers with good knowledge on general relativity.