Mass spring algorithm:
This 1D example shows the motion of a mass spring system that consists of a static and a dynamic particle linked by a spring. In each simulation step the following steps are performed:
- compute spring forces
- time integration to get new particle positions and velocities
On the left the simulation is shown while on the right we see the plot of the motion in x-direction compared to the analytic solution.
1. Compute spring forces
In this simple 1D example we have a spring with a stiffness $k$ and a rest length of 0. Therefore, the spring force is detertmined as
$$f(x) = -k x.$$
2. Time integration
Analytic solution
The analytic solution for the position $x$ of a particle which is linked by a spring to the origin is
$$
x(t) = x(0) \cos\left(\sqrt{\frac{k}{m}} t \right ),
$$
where $m$ is the mass of the particle and $x(0)$ the initial position.
Explicit Euler
The explicit Euler method is a first-order accurate method for solving ordinary differential equations (ODEs):
$$\begin{align*}
x(t + \Delta t) &= x(t) + \Delta t v(t) \\
v(t + \Delta t) &= v(t) + \Delta t \frac{f(t)}{m}.
\end{align*}$$
Symplectic Euler
The symplectic Euler method is also known as semi-implicit Euler since first the new velocity is determined and then the new position is computed using the new velocity value:
$$\begin{align*}
v(t + \Delta t) &= v(t) + \Delta t \frac{f(t)}{m} \\
x(t + \Delta t) &= x(t) + \Delta t v(t + \Delta t).
\end{align*}$$
The method is also first-order accurate but yields better results than the explicit Euler since it is a symplectic method.
Runge-Kutta 2
Typically numerical integration methods are defined for a differential equation determined by a function $f(t, s(t))$, where $s$ defines the state of the system at time $t$. In our 1D example the state is defined by the current position and velocity of the particle and the function is defined as:
$$\begin{equation*}
f(t, s(t)) = \begin{pmatrix} \dot{x} \\ \dot{v} \end{pmatrix}, \quad\quad s(t) = \begin{pmatrix} x(t) \\ v(t) \end{pmatrix}
\end{equation*}$$
Using this function the second-order Runge-Kutta integration is defined as
$$\begin{align*}
k_1 &= \Delta t f(t, s(t)) \\
k_2 &= \Delta t f(t + \frac12 \Delta t, s(t) + \frac12 k_1) \\
s(t + \Delta t) &= s(t) + k_2.
\end{align*}$$
Note that the Runge-Kutta method has multiple stages (in our case 2) to achieve a higher-order accuracy. However, this approach is more expensive than the Euler methods since the function has to be evaluated once per stage.
Implicit Euler
The implicit Euler method is a first-order accurate method for solving ordinary differential equations (ODEs):
$$\begin{align*}
x(t + \Delta t) &= x(t) + \Delta t v(t + \Delta t) \\
v(t + \Delta t) &= v(t) + \Delta t \frac{f(x(t + \Delta t))}{m}.
\end{align*}$$
The method is unconditionally stable but in general a non-linear system has to be solved.
In our simple 1D example we linearize the spring force $f(x)$ using a first-order Taylor approximation:
$$
f(x(t + \Delta t)) = f(x + \Delta x) = f(x) + \frac{\partial f(x)}{\partial x} \Delta x,
$$
where $\frac{\partial f(x)}{\partial x} = -k$ and we write $x=x(t)$ and $x + \Delta x = x(t + \Delta t)$ for simplicity. Using this linearization and transforming the above equations yields the linear equation
$$
\left (1 + \frac{\Delta t^2}{m}k \right ) \Delta v = \frac{\Delta t}{m} (f(x) - \Delta t k v),
$$
which can be simply solved for the velocity change $\Delta v$. Finally, the position and velocity are updated as
$$\begin{align*}
v(t+\Delta t) &= v(t) + \Delta v \\
x(t+\Delta t) &= x(t) + \Delta t v(t+\Delta t).
\end{align*}
$$
|