Just here to prove it to myself while noodling around with some code.
Recall the definition of a taylor series expansion with the normal caveats of smoothness etc. (Also, check out Madhava, Jason Gregory and Newton for historical development; it's also amusing to think about the sheer volume of Stigler's law occurences throughout history; its also heady and humbling to consider that this was initially used to calculate Halley's Comet trajectories and now it's a back of a napkin exercise )
\[ \sum_{n=0}^{\infty} \frac{f^{n}(a)}{n!} \left( t - a \right)^n \]
Where \(a\) is the point around which we're expanding and \(t\) is the variable we're creating the expansion for.
For a simulation, with a discrete time step associated with the physics tick/ cycle, \(\Delta{t}\), we need to be careful about what and where we're expanding.
We want to know what happens at the next time step: \(x\left( t + \Delta{t} \right)\) (This means we're creating an expansion for a function that can talk about that point.)
We know what is happening at the current step, so our expansion is expanding around that point: \(a = t\)
Subsituting into the Taylor series formula: (Note: \(\left( t - a \right) \rightarrow \left( t + \Delta{t} - t \right) = \left( \Delta{t} \right)\) )
\[ \implies \sum_{n=0}^{\infty} \frac{f^{n}(t)}{n!} \left( \Delta{t} \right)^n \]
This is if it's a forward time step. If it's a backward time step, all that would change is that we'd be expanding for \(x\left( t - \Delta{t} \right)\)
\[ \implies \sum_{n=0}^{\infty} \frac{f^{n}(t)}{n!} \left( -\Delta{t} \right)^n \]
The trick or method is to sum (simple verlet) or take the difference (stoermer-verlet) of the two.
Summing allows the even powers in the expansion to cancel which gives the method its stability (\(E \propto O^4\), i.e. the threshold for ignored higher ordered terms).
Expanded, truncated and summed: \[\implies x\left( t + \Delta{t} \right) + x\left( t - \Delta{t} \right) = 2x\left( t \right) + a\left( t \right) {\left( \Delta{t} \right)}^2\]
\[\implies x\left( t + \Delta{t} \right) = 2 x\left(t \right) - x\left( t - \Delta{t} \right) + a\left( t \right) {\left( \Delta{t} \right)}^2\]
We can talk about the position at the next time step if we know the current position, the last position and the current acceleration resolution.
Misc. Resources: More detailed article at the Algorithm Archive Here's a video to give some intuition by reworking the final Verlet equation into something more familiar: https://www.youtube.com/watch?v=lS_qeBy3aQI