Numerical Integration for Particles

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.

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)}$$

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} \end{align*}$$

   $$= \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}$$

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.