Finding energy minimizing path I'm trying to find an approximation for the optimal path for a material point, minimizing the integral associated with the total energy.
I managed to write the exact formula for the energy along a path covering the points $(x,y(x))$, which should be:
$$ E(y) = \int_{x_0}^{x_f} \sqrt{a-by(x)}\sqrt{1+(y^\prime(x))^2} \ \ \mathrm{d}x $$
Where $a$ and $b$ are constants related to mass, initial position and speed, and gravitational constant.
If I knew the path, numerically computing the integral would be trivial, but I don't know any methods to estimate
$$ y = \arg\min_{z} J(z) $$
I suppose that $y$ should be at least a $\mathcal{C}^1$ function (possibly $\mathcal{C}^\infty$?), but the only way to approach the problem in a general way that I can think about is to use a (linear) spline function, replacing differentiation with incremental ratios.
I could start with a straight line path, split it in $n$ parts and find a greedy optimization for each segment of the spline, moving the node $(x_1,y_1)$ upwards or downwards along the $y$-axis until an approximate minimum is found, then repeating the same with the node $(x_2,y_2)$, but this is obviously not going to work: .
I could instead try to move all the nodes at once (except the starting and ending ones, fixed at $(x_0,y_0)$ and $(x_f,y_f)$), but this seems too computationally intensive to pursue when using a finer grid.
Is there an efficient method to find an optimal path?
 A: Under the assumption that your solution is smooth enough, you can use the calculus of variations to turn your minimization problem into a differential equation. The differential equation is then given by the Euler-Lagrange equation, and can be solved by standard methods for numerical simulation of differential equations.
A: It is not hard to show, assuming that $b>0$, that the smooth solutions of the Euler-Lagrange equation are given by the formula
$$
y(x) = \frac{a}{b} - c_1^2 - \left(\frac{x-c_0}{2c_1}\right)^2 < \frac{a}{b}
$$
where $c_0$ and $c_1\not=0$ are constants.  For example, see my answer to the MO question Riemannian surfaces with an explicit distance function?.  Along such a path, it is easy to compute the energy using your integration.
However, you should note that, when $4(a-b y_0)(a-b y_f) < b^2(x_0-x_f)^2$, the actual minimizing path is to go straight from $(x_0,y_0)$ to $(x_0,a/b)$, then go along the line $y=a/b$ to $(x_f,a/b)$, and then go straight to $(x_f,y_f)$.  For an explanation of this phenomenon, see the above mentioned answer.
When $b<0$, one has to flip some signs, but essentially the same analysis holds.  When $b=0$, the minimizing curves are straight line segments.
A: The following iterative algorithm may provide an efficient way of iteratively generating estimates with increasing precision for the optimal path:
As the optimal path must be contained in a box $[x_0-c \lt x\lt x_f+c]\times[-\infty\lt y_{min} \le y\le \frac{a}{b}]$ for sufficiently large $c$ and sufficiently small $y_{min}$, we can restrict our considerations to that strip and proceed as follows:  


*

*estimate the values of $c$ and $y_{min}$  

*choose a $2D$-distribution that can be truncating to an arbitrary $2D$-interval.  

*Start: sample the distribution until the point-density inside $[x_0-c,\ x_f+c]\times[y_{min},\ \frac{a}{b}]$ is high enough and add points $(x_0,y_0)$ and $(x_f,y_f)$.  

*Repeat until the geometry of the calculated shortest path is precise enough:  


*

*construct the point set's Delaunay Triangulation and replace each edge by a pair of antiparallel arcs  

*calculate the arc weights via the formula supplied in the statement of the problem.  

*calculate the minimum-weight path in the modified Delaunay Triangulation from $(x_0,y_0)$ to $(x_f,y_f)$  

*cover the (undirected) edges of the shortest path with sufficiently small rectangular boxes that can be e.g., determined with an algorithm of raster graphics.  

*let the path-covering boxes be $\lbrace B_1,\ \dots,\ B_k\rbrace$, to which weights $\lbrace beta_1,\ \dots,\beta_k\rbrace$ are assigned, that equal the integral of the distribution's density function over the corresponding box.  

*assign to box $B_i$ the right-open range $\frac{\left[\sum_{j\lt i}\beta_j,\ \sum_{j\le i}\beta_j\right)}{\sum_{j\le k}\beta_j} $  

*repeat until the point-density in the path-covering boxes is high enough: 


*

*determine a random point from $\left[0,1\right)$ by sampling an equal distribution over that range; let $B_j$ be the path-covering box, that corresponds to the containing interval.  

*sample the point-generating distribution's restriction to $B_j$   
One should however be cautious to use a shortest-path algorithm, that works with negative edge-weights, such as the Bellman-Ford algorithm;  there is however no need to deal with negative cycles, because of the physical background of the arc-weights.
