# 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.

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$$.

• Thank you for taking the time to answer. I mentioned I don't have an ending point. I'd like to solve for the ending point and velocity. Ending point in this case is on the ground, so x=?, y=?, z=0. Can you comment on why direct minimization of the action functional starting with a position and velocity doesn't seem to be common compared to numerical integration of the equations of motion? Apr 2, 2020 at 20:12
• @user55937 I am not sure what you are trying to achieve then. First of all, for ballistic trajectories, if you are considering just constant vertical gravity, then everything is explicit, it's simply a parabolic trajectory. If you are talking about a more complicated situation, the minimization in physics is not with respect to the initial and final point (they are kept fixed, as are the initial and the final moment). The minimization is with respect to all trajectories connecting these two fixed points (and time moments) . You do not vary the endpoints, but you vary the curves. Apr 2, 2020 at 21:39
• @user55937 In general, solving initial value problems (starting position and velocity) is easier than solving boundary value problems (two fixed points). If you have a starting point and velocity, then you should work with initial value problem and track the trajectory until it reaches the surface you want. Apr 2, 2020 at 21:44
• direct minimization of the action functional starting with a position and velocity doesn't seem to be common, because the Lagrange minimization problem is not formulated this way. It is formulated with a starting and ending point, not with a starting point and velocity. The global minimization problem leads to the local Euler-LAgrange equations, which are easier to handle with initial values, rather than endpoints. Apr 2, 2020 at 22:11