Suppose I am at a point $p$ in a (smooth, connected, complete) Riemannian manifold $(M, g)$ and want to find a minimizing path to point $q \in M$.
Obviously I should travel along a geodetic and at least locally I know how geodetics look like, I just have to solve some ODEs. (For the sake of this question let us assume that is easily done.)
If $p, q \in U$ for a chart $(\kappa, U)$, then I am done, I simply choose a ”correct“ geodetic.
However, if $q$ is “farther away”, I do not know how to proceed. I could follow some (randomly chosen?) geodetic for a while. After leaving $U$ I can extend the geodetic by solving another ODE and then travel further, repeating this process if necessary. Problem: Even if I reach $q$ at some point, I may not have traveled a minimizing path.
Worse, there are uncountable many geodetics to choose from. Thus, the probability that this approach succeeds should be $0\%$.
Is there any terminating algorithm to find a minimizing path from $p$ to $q$, if one only knows $(M, g)$ (as well as $p$ and $q$, of course)? For any given point I am able to calculate locally all geodetics going through that point (for the cost of a single operation), but for non-local statements I need to repeat that process.
Does this task become easier for compact manifolds? Then there is at least a finite atlas …
This answer is mostly an attempt at writing down an arbitrary chunk of the theory of variational calculus on the fly while learning it. The two major references this answer is based on are (1) Milnor, "Morse theory" (part III) (2) Oancea, "Morse theory, closed geodesics, and the homology of free loop spaces" (chapter 2-3). Thanks to @MikeMiller for pointing out to the second reference which has the relevant theorems and constructions to make gradient flow work, and thanks to @0celo7 for giving a concise and clear introduction to the theory of Sobolev spaces that is essential to the second half of the answer.
Suppose $(M, g)$ is the Riemannian manifold and $p, q$ two points on $M$. Denote by $\Omega_{p, q}M$ to be the space of smooth paths $\gamma : [0, 1] \to M$ such that $\gamma(0) = p$ and $\gamma(1) = q$. This should be thought as an infinite dimensional manifold; if $\gamma \in \Omega_{p, q}M$ is a point one can define the tangent space $T_\gamma \Omega_{p, q}M $ to be the space of vector fields $X$ along $\gamma$ such that $X(p) = X(q) = 0$. The intuition is, if we have such a vector field $X$ along a path $\gamma$ in $M$, one can ask for solutions $\alpha : (-\epsilon, \epsilon) \times [0, 1] \to M$ of the "variational equation" given by $$\dfrac{d}{ds}\alpha(0, t) = X(\gamma(t))$$ with the initial conditions $\alpha(s, 0) = p$ and $\alpha(s, 1) = q$ for all $s$ (this happens because $X$ is zero at the endpoints of $\gamma$) and $\alpha(0, t) = \gamma(t)$ for all $t$. Therefore $\gamma_s = \alpha(s, -)$ are simply variations of $\gamma$ for $s\neq 0$ obtained from "flowing a little to the sideways along $X$". There is always a unqiue such variation for any given vector field $X \in T_\gamma \Omega_{p, q}M$ by the Picard-Lindelof theorem.
There is a function $\mathcal{E} : \Omega_{p, q}M \to \Bbb R$ known as the "energy functional", which is defined by $$\mathcal{E}(\gamma) = \displaystyle \int_0^1 g(\gamma'(t), \gamma'(t))dt$$ (note that this is similar to the formula for arclength) The "first derivative" $d\mathcal{E}$ of this functional is $d\mathcal{E}(\gamma)(X) = d/ds\,\mathcal{E}(\gamma_s)|_{s=0}$ for any given $\gamma$ in $\Omega_{p, q}M$ and $X$ in $T_\gamma \Omega_{p, q}M$ where $\gamma_s(t) = \alpha(s, t)$ is the variation of $\gamma$ defined by $X$.
Theorem: A path $\gamma$ in $\Omega_{p, q}M$ is a critical point of the energy functional (i.e., $d\mathcal{E}(\gamma)(X) = 0$ for all $X$ in $T_\gamma \Omega_{p, q}M$) if and only if $\gamma$ is a geodesic of $M$
Proof: Let's denote $\partial_s$ to be the derivative in the direction of $X$ and $\partial_t$ to be the derivative in the direction of $\gamma'$ for brevity's sake. Notice that $[\partial_t, \partial_s] = 0$. By differentiating under the integral sign, $$\begin{align} d\mathcal{E}(\gamma)(X) = \int_0^1 \partial_s g(\partial_t \gamma_s(t), \partial_t \gamma_s(t))dt & = 2\int_0^1 g(\nabla_{\partial_s} \partial_t \gamma_s(t), \partial_t \gamma_s(t))dt \\ &= 2\int_0^1 g(\nabla_{\partial_t} \partial_s \gamma_s(t), \partial_t \gamma_s(t)) dt\\&= 2\int_\gamma g(\nabla_{\gamma'} X, \gamma') \end{align}$$ where we have used symmetry of the Levi-Civita connection to conclude $\nabla_{\partial_t} \partial_s = \nabla_{\partial_s} \partial_t$. Using integration by parts (and using $X(\gamma(0)) = X(\gamma(1)) = 0$), we conclude $$d\mathcal{E}(\gamma)(X) = - 2\int_\gamma g(X, \nabla_{\gamma'} \gamma')$$ This formula implies that $d\mathcal{E}(\gamma)(X) = 0$ if and only if $\nabla_{\gamma'} \gamma' = 0$ (i.e. $\gamma$ is a geodesic). $\blacksquare$
Recall that if $f : M \to \Bbb R$ is a smooth function on a Riemannian manifold $(M, g)$ then the gradient vector field $\text{grad} \, f$ is defined by $g(\text{grad}\, f, X) = X(f)$. By Cauchy-Schwarz inequality, $g(\text{grad}\, f, X) \leq \|\text{grad}\, f\|_g \|X\|_g$ with equality iff $\text{grad} \, f = c X$ for some scalar function $c$. In other words, the maximum direction of ascent for $f$ is given by the gradient. Therefore if one looks at the negative gradient flow $\phi : \Bbb R \to \text{Diffeo}(M)$ given by $\partial_t \phi = -\text{grad}\, f$ the flowlines $\{\phi_t(x_0)\}$ follow the direction of maximal descent of $f$ (this algorithm is therefore known as the method of steepest descent). If $\|\text{grad} \, f\|_g$ decreases to zero along a certain flowline which does not escape to infinity, it must converge to a critical point of $f$.
To apply similar ideas to find the critical points of $\mathcal{E}$, we upgrade $\Omega_{p, q}M$ to a larger space $\widetilde{\Omega_{p, q} M}$. First, the setup: $H^1([0, 1], \Bbb R^n) \subset L^2([0, 1], \Bbb R^n)$ be the subspace of paths $\gamma : [0, 1] \to \Bbb R^n$ admitting a weak first derivative $\gamma'$ such that the Sobolev norm $\|\gamma\|_{H^1} := \left(\int \|\gamma\|^2 + \|\gamma'\|^2 \right)^{1/2}$ is finite. This is a Hilbert space with inner product $\langle \gamma_0, \gamma_1 \rangle_{H^1} = \int (\langle \gamma_0, \gamma_1 \rangle + \langle \gamma'_0, \gamma'_1 \rangle)$ inducing the $H^1$-norm. By Nash embedding theorem, realize $M$ as an isometrically embedded submanifold of $\Bbb R^n$ for some sufficiently large $n$. Then define $\widetilde{\Omega_{p, q}M}$ as a subspace of $H^1([0, 1], \Bbb R^n)$ consisting of $H^1$-paths $\gamma : [0, 1] \to \Bbb R^n$ with $\gamma(t) \in M$ for almost all $t$ such that $\gamma(0) = p$ and $\gamma(1) = q$. This turns out to be a Hilbert manifold (I don't see how to prove this) Tangent space $T_\gamma \widetilde{\Omega_{p, q} M}$ at a point $\gamma$ is defined to be the space of $H^1$-vector fields along $\gamma$ which are zero at the endpoints as before. The energy functional $\mathcal{E}$ is a $C^1$-functional on $\widetilde{\Omega_{p, q} M}$ and the previous formula for $d\mathcal{E}$ holds. The critical points of $d\mathcal{E}$ are now $H^1$-solutions to $\nabla_{\gamma'} \gamma' = 0$, which is an elliptic PDE of one variable , hence by elliptic regularity $H^1$-solutions are automatically $C^\infty$, therefore are exactly the geodesics in $M$. $\widetilde{\Omega_{p, q} M}$ has the structure of a Hilbert manifold with Riemannian metric induced from the $H^1$-inner product $$\tilde{g}(X, Y) = \int_\gamma g(X, Y) + \int_\gamma g(\nabla_{\gamma'} X, \nabla_{\gamma'} Y)$$ We can now define the gradient $\text{grad}\, \mathcal{E}$ so that $\text{grad}\, \mathcal{E}(\gamma)$ is the unique vector field along $\gamma$ so that $\tilde{g}(X, \text{grad}\, \mathcal{E}(\gamma)) = d\mathcal{E}(\gamma)(X)$ for all $X \in T_\gamma \widetilde{\Omega_{p, q} M}$. Comparing this with the previous formulas, $$\int_\gamma g(X, \text{grad}\, \mathcal{E}(\gamma)) + \int_\gamma g(\nabla_{\gamma'} X, \nabla_{\gamma'} \text{grad}\, \mathcal{E}(\gamma)) = -2\int_\gamma g(X, \nabla_{\gamma'} \gamma')$$ using integration by part on the middle integral along with $X|_{\partial \gamma} = 0$ we derive $$\int_\gamma g(X, \text{grad}\, \mathcal{E}(\gamma)) - \int_\gamma g(X, \nabla^2_{\gamma'} \text{grad}\, \mathcal{E}(\gamma)) = -2\int_\gamma g(X, \nabla_{\gamma'} \gamma')$$ Which implies $\text{grad}\, \mathcal{E}(\gamma)$ is the unique solution to the ODE $Y - \nabla_{\gamma'}^2 Y + 2 \nabla_{\gamma'} \gamma' = 0$ with initial condition $Y|_{\partial \gamma} = 0$ (or rather $Y(p) = Y(q) = 0)$. This determines $\text{grad}\, \mathcal{E}$ on every smooth path completely independently of the Hilbert manifold formalism of $\widetilde{\Omega_{p, q} M}$.
Theorem: The energy functional $\mathcal{E} : \widetilde{\Omega_{p, q} M} \to \Bbb R$ satisfies the Condition C, namely, if $(\gamma_k)$ is a sequence of paths in $\widetilde{\Omega_{p, q} M}$ such that $(\mathcal{E}(\gamma_k))$ is bounded and $\|\text{grad}\, \mathcal{E}(\gamma_k)\|_{\tilde{g}} \to 0$ then every limit point $\gamma$ of $(\gamma_k)$ is a critical point of $\mathcal{E}$, i.e., $d\mathcal{E}(\gamma)(X) = 0$ for all $X$ in $T_\gamma \widetilde{\Omega_{p, q} M}$.
(as a note here, it's interesting to understand why we took the energy functional and not the arclength functional $\mathcal{L} : \Omega_{p, q} M \to \Bbb R$ defined by $\mathcal{L}(\gamma) = \int_\gamma \sqrt{g(\gamma', \gamma')}$. The issue is that $\Omega_{p, q} M$ has a natural "gauge symmetry" given by action of the orientation preserving diffeomorphism group $\text{Diff}_{+}(I)$ acting by reparametrization on the paths in $\Omega_{p, q} M$. $\mathcal{L}$ preserves the gauge symmetry, i.e., $\mathcal{L}(\gamma \circ f) = \mathcal{L}(\gamma)$ for any $f \in \text{Diff}_+(I)$ ("arclength is independent of parametrization"). This means critical points of $\mathcal{L}$ gives rise to critical $\text{Diff}_+(I)$-orbits in $\Omega_{p, q} M$. But then there's no hope of having the above theorem hold for $\mathcal{L}$ as Condition C immediately implies the critical sets are compact. Passing to the energy functional $\mathcal{E}$ implicitly mods out by the gauge action as $\Omega_{p, q} M /\text{Diff}_+(I)$ can be identified with the subspace $\mathcal{A} \subset \Omega_{p, q} M $ of arclength parametrized curves between $p$ and $q$ in $M$, where the two functionals $\mathcal{E}$ and $\mathcal{L}$ agree. Simultaneously $\mathcal{E}$ as a functional on $\Omega_{p, q} M$ is not $\text{Diff}_+(I)$-invariant, so it "breaks the gauge symmetry".)
Therefore if $\gamma_s$ is a variation of $\gamma \in \widetilde{\Omega_{p, q} M}$ such that $\partial_s \gamma_s = -\text{grad}\, \mathcal{E}(\gamma_s)$, i.e., $\gamma_s$ is a flowline of the "negative gradient flow on $\widetilde{\Omega_{p, q}}$" starting at $\gamma$, then under a suitable choice of $\gamma$ we could hope for $\gamma_s$ to converge to a critical point of $\mathcal{E}$ by the previous theorem. Moreover if $\gamma_s$ converges to a minimum of $\mathcal{E}$, then it'd be a minimal geodesic as Cauchy-Schwarz implies $\mathcal{L}(\gamma)^2 \leq \mathcal{E}(\gamma)$. This gives an algorithm analogous to the steepest descent for finding (minimal) geodesics between two points $p, q$ of $M$ by flowing along the negative gradient in the space of paths.