Explicit Euler§

  • main idea: use the tangent at the start of the step
  • has issues with stability; simulation will eventually explode

Let $t^{(k)}$ denote the kth timestep, then the slope of a function $y$ at $t^{(k)}$ is given by using the finite difference operator $\Delta \left[y^{(k)}\right]$,

$$ \Delta\left[y^{(k)}\right] = \frac{y^{(k + 1)} - y^{(k)}}{t^{(k + 1)} - t^{(k)}} $$

From here on, we let $h \triangleq t^{(k + 1)} - t^{(k)}.$

Then, rearranging the expression yields

$$ y^{(k + 1)} = y^{(k)} + h \Delta\left[y^{(k)}\right] $$

Our system in the continuous case can be described in terms of a derivative, where $\dot{x}$ denotes the change w.r.t. time.

$$ \begin{equation} \begin{bmatrix} \dot{x}\\ \dot{v} \end{bmatrix} = \begin{bmatrix} v\\ m^{-1}\vec{f}(x, v) \end{bmatrix} \end{equation} $$

Rewriting this expression with discretized timesteps gives

$$ \begin{equation} \begin{bmatrix} x^{(k + 1)}\\ v^{(k + 1)} \end{bmatrix} = \begin{bmatrix} x^{(k)}\\ v^{(k)} \end{bmatrix} + h \begin{bmatrix} v^{(k)}\\ m^{-1}\vec{f}\left(x^{(k)}, v^{(k)}\right) \end{bmatrix} \end{equation} $$

Explicit Euler

for all particles $i$ do:   $v_i^{(k + 1)} \leftarrow v_i^{(k)} + hm^{-1}\vec{f}\left(x^{(k)}, v^{(k)}\right)$   $x_i^{(k + 1)} \leftarrow x_i^{(k)} + hv^{(k)}$

We let the function $\vec{f}$ denote the forces acting upon each particle, for example gravity and drag forces:

$$ \vec{f}(\vec{x}, \vec{y}) = \underbrace{m\vec{g}}_{\mathrm{gravity}} - \underbrace{D\vec{v}}_{\mathrm{drag}}. $$

Symplectic Euler§

  • stable (no numerical dissipation or explosion)
  • very simple to implement (a single change from explicit Euler)

Symplectic Euler

for all particles $i$ do:   $v_i^{(k + 1)} \leftarrow v_i^{(k)} + hm^{-1}\vec{f}\left(x^{(k)}, v^{(k)}\right)$   $x_i^{(k + 1)} \leftarrow x_i^{(k)} + h\underbrace{v^{(k + 1)}}_{change}$

Implicit Euler§

  • main idea: use tangent at the end of the step (instead of start as with explicit)

We redefine the finite difference operation from explicit Euler as follows:

$$ \Delta\left[y^{(k + 1)}\right] = \frac{y^{(k + 1)} - y^{(k)}}{t^{(k + 1)} - t^{(k)}} $$

Then, discretizing Eq. (1) with this redefined operator

$$ \begin{equation} \begin{bmatrix} x^{(k + 1)}\\ v^{(k + 1)} \end{bmatrix} = \begin{bmatrix} x^{(k)}\\ v^{(k)} \end{bmatrix} + h \begin{bmatrix} v^{(k + 1)}\\ m^{-1}\vec{f}\left(x^{(k + 1)}, v^{(k + 1)}\right) \end{bmatrix} \end{equation} $$

Hooke’s Law

The force of a simple spring can be given by Hooke’s law:

$$ f(x, v) = -kx $$

for some spring constant $k.$

Explicit Euler: $$ \begin{bmatrix} x^{(k + 1)}\\ v^{(k + 1)} \end{bmatrix} = \begin{bmatrix} x^{(k)}\\ v^{(k)} \end{bmatrix} + h \begin{bmatrix} v^{(k)}\\ m^{-1}\left(-kx^{(k)}\right) \end{bmatrix} $$

Note: the only difference between explicit and symplectic Euler is that we use $v^{(k + 1)}$ instead of $v^{(k)}$ in computing $x^{(k)}.$

Implicit Euler: $$ \begin{bmatrix} x^{(k + 1)}\\ v^{(k + 1)} \end{bmatrix} = \begin{bmatrix} x^{(k)}\\ v^{(k)} \end{bmatrix} + h \begin{bmatrix} v^{(k + 1)}\\ m^{-1}\left(-kx^{(k + 1)}\right) \end{bmatrix} $$

Next, substitute $x^{(k+1)}$ in the expression for $v^{(k+1)}$

$$ \begin{align*} x^{(k+1)} &= x^{(k)} + hv^{(k+1)}\\ v^{(k+1)} &= v^{(k)} - \frac{hk}{m}x^{(k+1)}\\ &= v^{(k)} - \frac{hk}{m}\left(x^{(k)} + hv^{k + 1}\right)\\ &= v^{(k)} - \frac{hk}{m} x^{(k)} - \frac{h^2k}{m}v^{k + 1} \end{align*} $$

Finally,

$$ \begin{align*} \left(1 + \frac{h^2k}{m}\right)v^{(k+1)} &= v^{(k)} - \frac{hk}{m} x^{(k)}\\ x^{(k + 1)} &= x^{(k)} + hv^{(k + 1)} \end{align*} $$

Note: the only difference between implicit and symplectic Euler is that the coefficient on $v^{(k)}$ is equal to 1 in the symplectic case.

The subsittution and rearrangement allows us to easily compute (or precompute) the coefficient for each step.

Hooke’s Law, Gravity, and Dampening

We can easily combine different forces, for example, spring force (Hooke’s), gravity, and dampening:

$$ f(x, v) = mg - kx - Dv $$

The explicit and symplectic cases are trivial.

Implicit Euler: The velocity update is given by

$$ v^{(k+1)} = v^{(k)} - m^{-1}\left(mg - kx^{(k+1)} - Dv^{(k+1)}\right) $$

Next, we substitute $x^{(k + 1)} = x^{(k)} + hv^{(k+1)}$

$$ \begin{align*} v^{k + 1} &= v^{(k)} + \frac{h}{m} \left(mg - k\left(x^{(k)} + h v^{(k + 1)}\right) - Dv^{(k + 1)}\right)\\ &= v^{(k)} + hg - \frac{h}{m}kx^{(k)} - \frac{h^2}{m}kv^{(k + 1)} - \frac{h}{m}Dv^{(k + 1)} \end{align*} $$

Rearranging terms gives

$$ \left(1 + \frac{h^2}{m}k + \frac{h}{m}D\right) v^{(k + 1)} = v^{(k)} + hg - \frac{h}{m}kx^{(k)} $$

Multiple Interacting Particles§

Two Particle Spring

Suppose we have two particles at each end of a simple spring.

Then, the forces acting on each particle $f_0$ and $f_1$ are equal and opposite, that is

$$ \begin{align*} f_0(x_0, x_1) &= k(x_1 - x_0)\\ f_1(x_0, x_1) &= k(x_0 - x_1)\\ &= -f_0(x_0, x_1) \end{align*} $$

Symplectic Euler:

$$ \begin{align*} v_0^{(k + 1)} &= v_0^{(k)} + hm_0^{-1}k\underbrace{\left(x_1^{(k)} - x_0^{(k)}\right)}_{f_0^{(k)}}\\ v_1^{(k + 1)} &= v_0^{(k)} + hm_1^{-1}k\underbrace{\left(x_0^{(k)} - x_1^{(k)}\right)}_{f_1^{(k)}}\\ x_0^{(k + 1)} &= x_0^{(k)} + hv_0^{(k + 1)}\\ x_1^{(k + 1)} &= x_1^{(k)} + hv_1^{(k + 1)} \end{align*} $$

Implicit Euler: $$ \begin{align*} v_0^{(k + 1)} &= v_0^{(k)} + hm_0^{-1}k\left(x_1^{(k + 1)} - x_0^{(k + 1)}\right)\\ v_1^{(k + 1)} &= v_0^{(k)} + hm_1^{-1}k\left(x_0^{(k + 1)} - x_1^{(k + 1)}\right)\\ x_0^{(k + 1)} &= x_0^{(k)} + hv_0^{(k + 1)}\\ x_1^{(k + 1)} &= x_1^{(k)} + hv_1^{(k + 1)} \end{align*} $$

As with the other implicit Euler examples, we still wish to compute the coefficient term for the above example.

Let $M_i = m_i \mathbf{I}_3$ be the “mass matrix”.

  1. Multiply everything by the mass matrix $$ M_0v_0^{(k + 1)} = M_0v_0^{(k)} + hf_0^{(k + 1)}\left(x_0^{(k + 1)}, x_1^{(k + 1)}\right)\\ M_0v_1^{(k + 1)} = M_0v_1^{(k)} + hf_1^{(k + 1)}\left(x_0^{(k + 1)}, x_1^{(k + 1)}\right) $$

  2. Expand forces $$ \begin{align*} f_0^{(k + 1)} &= k \left(x_1^{k + 1} - x_0^{k + 1}\right)\\ &= k\left(\left(x_1^{(k)} + h v_1^{(k + 1)}\right) - \left(x_0^{(k)} + h v_0^{(k + 1)}\right)\right)\\ &= k\left(x_1^{(k)} - x_1^{(k)}\right) + h \left(v_1^{(k + 1)} - v_0^{(k + 1)}\right)\\ f_1^{(k + 1)} &= -f_0^{(k + 1)} \end{align*} $$

  3. Plug in expanded forces $$ \begin{align*} M_0v_0^{(k + 1)} = M_0v_0^{(k)} + h\left(f_0^{(k)} + hk \left(v_1^{(k + 1)} - v_0^{(k + 1)}\right)\right)\\ M_0v_1^{(k + 1)} = M_0v_1^{(k)} + h\left(f_1^{(k)} + hk \left(v_0^{(k + 1)} - v_1^{(k + 1)}\right)\right) \end{align*} $$

  4. Rewriting in matrix notation $$ \begin{align*} \begin{pmatrix} M_0 & 0\\ 0 & M_1 \end{pmatrix} \begin{pmatrix} v_0^{(k + 1)}\\ v_1^{(k + 1)} \end{pmatrix} &= \begin{pmatrix} M_0 & 0\\ 0 & M_1 \end{pmatrix} \begin{pmatrix} v_0^{(k)}\\ v_1^{(k)} \end{pmatrix} + h \begin{pmatrix} f_0^{(k)}\\ f_1^{(k)} \end{pmatrix} + h^2 k \begin{pmatrix} v_1^{(k + 1)} - v_0^{(k + 1)}\\ v_0^{(k + 1)} - v_1^{(k + 1)} \end{pmatrix}\\ &= \begin{pmatrix} M_0 & 0\\ 0 & M_1 \end{pmatrix} \begin{pmatrix} v_0^{(k)}\\ v_1^{(k)} \end{pmatrix} + h \begin{pmatrix} f_0^{(k)}\\ f_1^{(k)} \end{pmatrix} + h^2 k \begin{pmatrix} -\mathbf{I}_3 & \phantom{-}\mathbf{I}_3\\ \phantom{-}\mathbf{I}_3 & -\mathbf{I}_3 \end{pmatrix} \begin{pmatrix} v_0^{(k + 1)}\\ v_1^{(k + 1)} \end{pmatrix} \end{align*} $$

  5. Rearrangement $$ \underbrace{\left( \begin{pmatrix} M_0 & 0\\ 0 & M_1 \end{pmatrix} - h^2 k \begin{pmatrix} -\mathbf{I}_3 & \phantom{-}\mathbf{I}_3\\ \phantom{-}\mathbf{I}_3 & -\mathbf{I}_3 \end{pmatrix} \right)}_{A} \underbrace{\begin{pmatrix} v_0^{(k + 1)}\\ v_1^{(k + 1)} \end{pmatrix}}_{x} = \underbrace{\begin{pmatrix} M_0 & 0\\ 0 & M_1 \end{pmatrix} \begin{pmatrix} v_0^{(k)}\\ v_1^{(k)} \end{pmatrix} + h \begin{pmatrix} f_0^{(k)}\\ f_1^{(k)} \end{pmatrix}}_{b} $$

    Again, use some library like Eigen to solve the linear system $Ax = b,$ the solution $x = \left(v_0^{(k + 1)}, v_1^{(k + 1)}\right)^T$ contains the velocity update.

In general, the implicit with an $n$ particle system yields the following

$$ \left(\mathfrak{M} - h^2 \mathfrak{K}\right) v^{k + 1} = \mathfrak{M}v^{(k)} + hf^{(k)} $$

where

$$ \mathfrak{M} = \begin{pmatrix} M_0 & 0 & \cdots & 0\\ 0 & M_1 & \cdots & 0\\ 0 & 0 & \ddots & 0\\ 0 & 0 & \cdots &M_n \end{pmatrix} $$

and $\mathfrak{K}$ is the stiffness matrix, which is the derivative of the force w.r.t $x$.

In our 2 particle system with Hooke’s law,

$$ \mathfrak{K} = \frac{\mathrm{d}f}{\mathrm{d}x} = k \begin{pmatrix} -\mathbf{I}_3 & \phantom{-}\mathbf{I}_3\\ \phantom{-}\mathbf{I}_3 & -\mathbf{I}_3 \end{pmatrix}. $$

Again, the symplectic is the same as the implciit without the stiffness.