Can you numerically compute a trajectory by direct minimization of the action functional? Is there a numerical approach to compute simple projectile motion by directly minimizing the action functional? 
I was thinking that the trajectory is essentially a least cost path through phase space where the traversal cost is the Lagrangian and the cost is accumulated over time. If the scenario was something like a catapult where the projectile is launched from ground level at a known velocity and lands at ground level, the boundary conditions would be $y(t_i)=0$, $v(t_i)=v_0$, and $y(t_f)=0$. 
Numerically, I was thinking the trajectory could be represented as a sequence of nodes in phase space (i.e., integration elements). The first node would be fixed at $(y_0, v_0)$ and the final node would be fixed on the ground, but free to move laterally. The trajectory would be described by the locations of the remaining nodes that minimize $\sum_i L(r_i, v_i)$ subject to the constraint that $\Delta r/ v$ is constant across all nodes. The number of nodes in a trajectory would not change. Consequently, the integration step $\Delta t=\Delta r/ v$ might be different for different trajectories or while the trajectory is being minimized, but that's ok so long as $\Delta t$ is uniform across all nodes in a trajectory. 

I haven't been successful performing this optimization because (I think) the constraint represents too small a region of the domain and the node locations are treated completely independently. I also tried representing the constraint as a clique potential to make the space smoother, but the minimization was still unstable. Just FYI, I've been using scipy's optimize and experimenting with different solvers. 
I started thinking about computing the stationary trajectory because I would like to simulate many similar trajectories. I was hoping that, in this approach, the optimal trajectory starting at $r_0, v_0$ would be a good starting point to optimize the trajectory for $r_0, v_0+\delta v$, and therefore there could be some computational gains.  Whereas when performing usual numerical integration, computing one trajectory doesn't buy you anything toward computing an adjacent trajectory. Also, I only want to know where and with what energy the projectile lands, the time dependence is less important.
I haven't been able to find much relevant literature, and so I suspect an approach like this just isn't practical for some reason (although maybe I just don't know the right search terminology). At first I was curious if this could work, but now I'm getting a little obsessed trying to figure out why this won't work given that it's not prevalent. 
 A: So, let us say you are looking for a trajectory $x = x(t) = \big(\, x^1(t), \, x^2(t), ...,\, x^n(t)\, \big)$ that connects two fixed points and times $x_0, t_0$ and $x_1, t_1$. Such a trajectory is an optimal trajectory for the action:
$$S[x] =  \int_{t_0}^{t_1}  L\Big(\,x, \,\frac{dx}{dt}\,\Big) \, dt = \int_{t_0}^{t_1}  L\Big(x(t), \,\frac{dx}{dt}(t)\Big)\, dt .$$ 
Where $L(x, \dot{x})$ is the lagrangian of the system.
Then, according to the calculus of variations, the optimizing curves of the functional $S$ provide the motion of the system with respect to time and are the solutions of the Euler-Lagrange equations $$\frac{d}{dt}\Big(\,  \nabla_{\dot{x}}L\Big(\,x, \,\frac{dx}{dt}\,\Big) \, \Big) = \nabla_{{x}}L\Big(\,x, \,\frac{dx}{dt}\,\Big) \, \Big)$$ 
One can go for a straight-forward discretization of the formalizm outlined above as follows: replace the derivative $\frac{dx}{dt}$ by a difference $(\tilde{x} - x)/h$ (or something even better if you wish) and consider the discrete Lagrangian 
$$L_{h}(x,\,\tilde{x}) := L\left( x, \, \frac{\tilde{x}-x}{h} \,
\right)$$
By analogy with the continuous case outlined above, the action in
the discrete case is
$$S_{h}[\hat{x}] = \sum_{k=0}^{N}
L_h\big(x_k,\,x_{k+1} \big)\,h$$  where $x_0$ is your initial point
and $x_{N+1}$ is your final point and the points that describe the
discrete trajectory form the multi-point $$\hat{x} = \big(\,x_1,
\, x_2,\, x_3,\, ..., x_k,\, ...,\, x_N\,\big) \, \in \, \mathbb{R}^{nN}$$
(recall, $x_k = \big(\,x_k^1,\, x_k^2,\, ..., x_k^n\,\big) \, \in \mathbb{R}^n$). Then, the critical discrete trajectory should simply be
the solution to the zero-gradient equations $\nabla S_{h}[\hat{x}]
= 0$ which componentwise leads basically to the discrete version
of the Euler-Lagrange equations
$$\nabla_{x_{k}} L_h\big(x_{k-1},\,
x_k\big) + \nabla_{x_{k}} L_h \big(x_{k},\, x_{k+1}\big) = 0
\,\,\, \text{ for } \,\,\, k=1,...,N$$
$$x_0 = \text{ fixed initial point, } \,\, x_{N+1} = \text{ fixed finial point. }$$
To avoid confusion, I am going to denote by $\nabla_1$  the
gradient derivatives of the lagrangian $L_h\big(x ,\,
\tilde{x}\big)$ with respect to the first set of variables $x$ and
by $\nabla_2$ gradient derivatives of $L\big(x ,\, \tilde{x}\big)$
with respect to the first set of variables $\tilde{x}$. Thus, the
discrete Euler-Lagrange equations become:
$$\nabla_{2} L_h\big(x_{k-1},\,
x_k\big) + \nabla_{1} L_h \big(x_{k},\, x_{k+1}\big) = 0 \,\,\,
\text{ for } \,\,\, k=1,...,N$$
$$x_0 = \text{ fixed initial point, } \,\, x_{N+1} = \text{ fixed finial point. }$$
The latter is a system of algebraic equations, something like $n
N$ equations and variables. The solution is a sequence of points
$\,\, x_0, \, x_1, \, x_2, \, ... \, , \, x_{N+1} \,$ which should
approximate the time-parametrized trajectory of the system between
$x_0$ and $x_{N+1}$. You can interpolate between consecutive
points if you want to get a smoother curve.
Now, I am going to deviate, and although you are not asking about
initial value problems, I will nevertheless discuss them, as I
believe they are conceptually important and may give you some
better insight. In order to solve discrete initial value problems,
i.e. to find discrete trajectories with a stating point and a
direction vector, one can look locally at the discrete
Euler-Lagrange equations and write them down using a simplified
superscript notation (I basically drop the index $k$ and keep only
the increments $-1$ and $1$).
$$\nabla_2 L_{h}\big(x_{(-1)},\, x\big)
+ \nabla_1 L_{h}\big(x,\, x_{1}\big) = 0$$
Introduce the variable $p = \nabla_1L_h \big(x,\, x_{1}\big)$.
Then the discrete Euler-Lagrange equations become
$$\nabla_2L_h\big(x_{(-1)}, \,x\big)\,
+\,p \,= \, 0$$ which shifted by one subscript, turn into
$\nabla_2L_h\big(x,\, x_1\big) + p_1= 0$. Thus we obtain the
equations
\begin{align*}
&p = \nabla_1L_h\big(x,\, x_{1}\big)\\
&p_1 = -\,\nabla_2L_h\big(x,\, x_1\big)
\end{align*}
If one can express $x_1$ as a function of $(x,\,p)$ from the first
equation, then the second equation also gives us $p_1$ as a
function of $(x,\,p)$. Thus, we can obtain a map $\Phi_h : (x,\,p)
\mapsto (x_1,p_1)$. Observe that this is a map $\Phi_h : T^*\mathbb{R}^n
\to T^*\mathbb{R}^n$, which turns out to be symplectic (a local
symplectomorphism), because the Lagrangian $L_h$ is in fact a
generating function of the symplectomorphism $\Phi_h$.
Then, given an initial point $x_0$ and a direction vector $v_0$,
take let's say something like $x_1 =x_0 + h \, v_0 $ and obtain
$p_0 = \nabla_1L_h\big(x,\, x_1\big)$. As a result, starting from
$x = x_0$ and $p=p_0$, iterate the map $\Phi_h$:
$$\big(x_{k+1},\, p_{k+1}\big) = \Phi_h\big(x_{k},\, p_{k}\big)$$
obtaining again a sequence $\,\, x_0, \, x_1, \, x_2, \, ... \, ,
\, x_{k}, \,\, ... ,\, x_{N+1}$ which approximates the
time-parametrized trajectory of the system (here you have another
sequence $\,\, p_0, \, p_1, \, p_2, \, ... \, , \, p_{k}, \,\,
...,\, p_{N+1}\,$ which is a sequence of momenta, dually related
to the tangent velocity vectors.
Coming back to the the discrete Euler-Lagrange equations from
before, looks like you are after the solutions of the
zero-gradient system of equations $$\nabla S_h[\hat{x}] = \nabla
S_h( x_1, x_2, ..., x_k, ... x_{N}) = 0$$ i.e. you are looking for
the multi-point $\hat{x} = ( x_1, x_2, ..., x_k, ... x_{N}) \, \in
\, \mathbb{R}^{nN}$ for which the gradient of $S_h\big( x_1,\, x_2,\,
...,\, x_k,\, ...,\, x_{N})$ is zero. Therefore, in the case of
given initial and final point for the geodesic, a first line
computational approach could be a gradient descent method (or a
version of Newton's method). In the case of an initial point and a
direction vector, simply iterate the map $\Phi_h$. Even if the map
is not explicit, the problem is still slightly simpler because one
goes step by step, every time solving only the system involving
the variables $x, x_1, p, p_1$.
To make the gradient descent method a bit more efficient, one can
choose a smart initial guess for the discrete trajectory $x_0,\,
x_1, \, ..., \, x_{N+1}$. We have $x_0$ and $x_{N+1}$ fixed. These
are the initial point and the end point. Denote by $\hat{x}(m) =
\big(x_0, \, x_1(m),\, x_2(m), \, ...,\, x_{N}(m), \, x_{N+1}
\big)$ a sequence of $N+2$ points at iteration $m$. Then say we
use some sort of gradient descent $$\hat{x}(m+1) = \hat{x}(m) -
\alpha_m \, \nabla S_h \big[\hat{x}(m)\big]$$ with starting point
$\bar{x}(0) = \big(x_0, \, x_1(0),\, x_2(0), \, ...,\, x_{N}(0),
\, x_{N+1}\big)$. If $\hat{x}(0)$ is picked carefully, then the
gradient descent could make less iterations before arriving at a
very good approximate solution to the problem $\nabla \,
S_h[\hat{x}] = 0$.
If you have some extra information, you can use $x_0$ as a
starting point and a guess for a second point $x_1$, then run the
initial value iteration, using $\Phi_h$, to create a discrete
geodesic you can use as an initial guess for your gradient descent
algorithm (or Newton's method or whatever numerical scheme you use
to solve the discrete variational problem). That's why I have also
discussed the initial value problem because one can use it as an
auxiliary tool. Or, you could combination of the two methods: one
(or more) step with the gradient descend and one trajectory
generated by iterating $\Phi_h$.
