A second order mechanical system can be described with:

$$ y + k_1\dot{y} + k_2\ddot{y} = x + k_3\dot{x} $$

Where:

• y is position
• $\dot{y}$ is velocity
• $\ddot{y}$ is acceleration
• $k_1$, $k_2$, $k_3$ are system describing constants

To solve this, we use the Semi-implicit Euler Method.

$$ y[0] = x_{0} $$

$$ \dot{y}[0] = 0 $$

Given the time between samples (or in this case, frames):

$$ y[n+1] = y[n] + T\dot{y}[n] $$

$$ \dot{y}[n+1] = \dot{y}[n] + T\ddot{y} $$

Solving the equation of motion for the acceleration ($\ddot{y}$), we get:

$$ \ddot{y} = (x + k_3\dot{x} - y - k_1\dot{y})/k_2 $$

We can then substitute our most recent estimations of position and velocity

$$ \ddot{y} = (x[n+1] + k_3\dot{x}[n+1] - y[n+1] - k_1\dot{y}[n])/k_2 $$

Substituting this into out velocity updating function, we get:

$$ \dot{y}[n+1] = \dot{y}[n] + T(x[n+1] + k_3\dot{x}[n+1] - y[n+1] - k_1\dot{y}[n])/k_2 $$

(Notice that we use $y[n+1]$ - updated - instead of the previous, $y[n]$)

If $[n+1]$ is unknown, it can be estimated as an average between new and last time sample:

$$ \dot{x}[n+1] = (x[n+1] - x[n])/T $$

As a result, we get the following equations:

$$ \dot{x}[n+1] = (x[n+1] - x[n])/T $$

$$ y[n+1] = y[n] + T\dot{y}[n] $$

$$ \dot{y}[n+1] = \dot{y}[n] + T(x[n+1] + k_3\dot{x}[n+1] - y[n+1] - k_1\dot{y}[n])/k_2 $$

Constants $k_1$, $k_2$ and $k_3$ can be described by more intuitive factors:

• $f = \frac{1}{2\pi\sqrt{k_2}}$
• $\zeta = \frac{k_1}{2\sqrt{k_2}}$
• $r = \frac{2k_3}{k_1}$

By substitution, we get the following definitions for $k_1$, $k_2$, $k_3$:

$k_1 = \frac{\zeta}{\pi f}$

$k_2 = \frac{1}{(2\pi f)^2}$

$k_3 = \frac{r\zeta}{2\pi f}$

If the resonant frequency $f$ is too high relative to the framerate (dt), the system will become unstable. Just keeping $f$ low is not enough, as a single lag spike may be enough to break the simulation.

A simple way of looking at it is looking at a first-ordere system:

$$ y[n+1] = A y[n]$$

Here, if $|A|> 1$, the next iteration will be larger than the previous. For this type of system to be stable, it should converge towards 0.

In the case of a second-order system, by substituting the update of position into the update of acceleration, we get a system of linear equations, which can be written in this form:

$$[\frac{y}{\dot{y}}]^{[n+1]} = A[\frac{y}{\dot{y}}]^{[n]} + B[\frac{x}{\dot{x}}]^{[n+1]} $$

Where A and B are matrices based on the k-constants and T.

Whether A is "greater than 1" (which doesn't really make sense for a matrix) is defined by its eigen values. However, for this implementation we'll clamp either T or k_2 instead:

$$ T < \sqrt{4k_2+k_1²} - k_1 $$

For $T&gt;$, the computation can be divided, however that will require and extra iteration of computation.

$$ k_2 > \frac{T^2}{4} + \frac{Tk_1}{2} $$

Clamping $k_2$ to be above this value is not physically correct, but will stabilize the system.

Jitter can also happen (caused by negative eigen values), and the following clamping of $k_2$ can remove this:

$$ k_2 > \frac{T^2}{2} + \frac{Tk_1}{2} \\ \text{and} \\ k_2 > T k_1 $$

• Eigen values and critical time step:
  • https://www.youtube.com/watch?v=KPoeNZZ6H4s
• Pole-zero matching for more stable and precise simulation.