Theorem: Derivation of Hamilton-Jacobi Equation

Let $\map S {x_0, x_1, \mathbf y} = \map S {x, \mathbf y}$ be the geodetic distance, where $x_0$ is fixed and $x_1=x$.
Let $H$ be Hamiltonian.
Then the following equation holds:
:$\dfrac {\partial S} {\partial x} + \map H {x, \mathbf y, \nabla_{\mathbf y} S} = 0$
and is known as the Hamilton-Jacobi Equation.

Proof:
Consider the increment $\Delta S$:
:$\Delta S = \map S {x + \Delta x, \mathbf y + \Delta \mathbf y} - \map S {x, \mathbf y}$
Note that the change of function $\mathbf y$ denoted by $\Delta \mathbf y$ is dependent on the manner $\Delta x$ is chosen through the definition of geodetic distance.
For sufficiently smooth $S$, $\size {\Delta \mathbf y} \to 0$ as $\size {\Delta \mathbf x} \to 0$.
{{explain|sufficiently}}
By definition of differential, $\Delta S$ can be written as:
:$\map {\Delta S} {x, \mathbf y; \Delta x, \Delta \mathbf y} = \map {\d S} {x, \mathbf y; \Delta x, \Delta \mathbf y} + \epsilon \Delta x + \boldsymbol \epsilon \cdot \Delta \boldsymbol y$
where $\epsilon \to 0$ as $\Delta x \to 0$, and $\size {\mathbf y} \to 0$ as $\size {\Delta \mathbf x} \to 0$.
By definition of the geodetic distance,
:$\Delta S = J \sqbrk {\gamma^*} - J \sqbrk \gamma$
where $\gamma$ and $\gamma^*$ are extremal curves, connecting the fixed initial point with points $\tuple {x, \mathbf y}$ and $\tuple {x + \Delta x, \mathbf y + \mathbf h}$ respectively.
By definition of increment of functional:
:$J \sqbrk {\gamma^*} - J \sqbrk \gamma = \Delta J \sqbrk {\gamma; \Delta \gamma}$
where $\Delta \gamma = \gamma^* - \gamma$.
A differentiable $J$ can be expressed as:
:$\Delta J \sqbrk {\gamma; \Delta \gamma} = \delta J \sqbrk {\gamma; \Delta \gamma} + \epsilon_\gamma \cdot \size {\Delta \gamma}$
where $\epsilon_\gamma \to 0$ as $\size {\Delta \gamma} \to 0$, and $\size {\Delta \gamma} \to 0$ as $\size {\Delta \mathbf x} \to 0$ for sufficiently smooth $S$.
To summarise:
:$\Delta \map S {x, \mathbf y; \Delta x, \Delta \mathbf y} = \Delta J \sqbrk {\gamma; \Delta \gamma}$
{{explain|change in $x$ should be reflected on both sides, but literature omits this}}
Both sides contain terms linear in $\size {\Delta x}$, $\size {\Delta \mathbf y}$, $\size {\Delta \gamma}$ as well terms of higher order.
Higher order terms on both sides are the same.
{{explain|why}}
Hence, the principal parts match:
:$\d S = \delta J$
The variation of extremal $J$ is expressible as
:$\ds \delta J = \sum_{i \mathop = 1}^n p_i \Delta y_i - H \Delta x$
while the differential of $S$ is 
:$\ds \d S = \frac {\partial S} {\partial x} \Delta x + \sum_{i \mathop = 1}^n \frac {\partial S} {\partial y_i} \Delta y_i$
Equivalently:
:$\ds \paren {\frac {\partial S} {\partial x} + H} \Delta x + \sum_{i \mathop = 1}^n \paren {\frac {\partial S} {\partial y_i} - p_i} \Delta y_i = 0$
$\Delta x$ and $\Delta y_i$ are independent variables.
The equation holds only if all the coefficients in front of $\Delta x$ and $\Delta y_i$ vanish simultaneously:
:$\dfrac {\partial S} {\partial x} = -H, \quad \dfrac {\partial S} {\partial y_i} = p_i$
Since $H = \map H {x, \mathbf y, \mathbf p}$, using the second relation to replace $\mathbf p$ together with the first one proves the formula.
{{finish|the proof makes sense, but very technical aspects are missing}}
{{Namedfor|William Rowan Hamilton|name2 = Carl Gustav Jacob Jacobi}}