I am trying to understand how to use the improved euler method on MPM simulations.

In the kind of MPM simulation I am doing with forward euler the order of operations is as follows:

  • Write particle material properties (i.e. mass) into the discrete gird cells.
  • Read the material properties from the grid to compute material forces based on physical properties (e.g. use velocity and grid density at the particle to compute a force in the direction that minimizes density the most). Write this force to the grid.
  • Update simulation properties in the grid (i.e. the forces/grid velocities).
  • For each particle, read the forces from the grid and update its position using forward euler, i.e. velocity = velocity + acceleration * time_step and position = position + velocity * time step.

I am trying to understand how to adapt the improved Euler method. My current understanding is that I run the simulation the way I am doing it now, update particle positions but store the prior ones as well as the prior forces.

Run the simulation again with the new particle positions and store the new forces.

Update the prior particle positions as:

velocity = velocity + (prior_vel + current_vel) / 2

position = position + (prior_pos + current_pos) / 2

But this is not giving me the best results so I must have messed it up.

  • $\begingroup$ You did apply the timestep in the second order scheme? $\endgroup$ Commented Jul 5, 2023 at 7:14
  • $\begingroup$ Could you write the form of your evolution equation(s) in the continuous setting? $\endgroup$
    – lightxbulb
    Commented Jul 5, 2023 at 12:48

1 Answer 1


I haven't studied the specifics of MPM and I can't tell the exact problem you're solving since you didn't formulate an evolution equation explicitly, but I can write the general form of the standard Euler methods, so hopefully that can help you.

Euler Methods

Consider the simplest case where you just have $d_tu(t) = f(t,u(t))$, with an initial condition $u(0) = u^0$. You want to discretise time on a grid $t_j = j\tau$ with time step $\tau$. That is, you are looking for a sequence $\tilde{u} = (\tilde{u}^0, \tilde{u}^1,\tilde{u}^2,\ldots)$, where $\tilde{u}^0 = u^0 = u(0)$, and $\tilde{u}^j$ aims to approximate $u(t_j) = u(j\tau)$, then the Euler discretisations are: \begin{alignat}{3} \tilde{u}^{j+1} &= \tilde{u}^j + \tau f(t_j, \tilde{u}^j), &\quad& (\textrm{explicit Euler}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \tau f(t_{j+1}, \tilde{u}^{j+1}), &\quad& (\textrm{implicit Euler}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \frac{\tau}{2}\bigl( f(t_{j}, \tilde{u}^{j}) + f(t_{j+1}, \tilde{u}^{j+1}) \bigr), &\quad& (\textrm{Crank-Nicolson}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \frac{\tau}{2}\bigl( f(t_{j}, \tilde{u}^{j}) + f(t_{j+1}, \tilde{u}^{j}+\tau f(t_j,\tilde{u}^j)) \bigr), &\quad& (\textrm{improved Euler}). \end{alignat} The simplest discretisation is explicit Euler, since it just involves evaluating $f$ for the known arguments $(t_j, \tilde{u}^j)$ in order to compute $\tilde{u}^{j+1}$. In contrast, implicit Euler requires solving a potentially nonlinear equation for $\tilde{u}^{j+1}$ (the latter is usually implemented with a fixed-point iteration or Newton-Rhapson). The main benefit of implicit Euler is that it is L-stable, while explicit Euler is not even A-stable. This generally means that you have severe restrictions on the step size $\tau$ in explicit Euler, while this isn't the case for implicit Euler. For example if we had a heat diffusion equation, $\tau$ can be taken arbitrarily large for implicit diffusion, but it must be less than $\frac{h^2}{2}$ for explicit Euler in 1-D with a standard approximation in space of $\partial_{xx}u$ with spatial step size $h$. Both explicit and implicit Euler are of first order in time (i.e.\ the error scales as $O(\tau)$). That's where Crank-Nicolson comes in which is second order in time $O(\tau^2)$. Since it involves a term with $\tilde{u}^{j+1}$ on the right-hand side, you have to solve an equation for $\tilde{u}^{j+1}$ similar to the implicit method. Since this may be undesirable, one can consider the improved Euler method which also has order $O(\tau^2)$ but doesn't require solving such an equation. This is achieved by replacing the $\tilde{u}^{j+1}$ term on the right-hand side by its approximation by an explicit Euler step: $\tilde{u}^{j+1} \to \tilde{u}^j + \tau f(t_j, \tilde{u}^j).$

Euler Methods For Higher Derivatives

In your question you mention that you update both the position and the velocity, which to me suggests that you have at least something involving $d_{tt}$: $d_{tt} u(t) = f(t, u(t), d_t u(t))$ with two initial conditions $u(0) = u^0$ and $d_t u(0) = v^0$ (e.g. initial position and velocity). Let $v(t)=d_t u(t)$, then we may rewrite the problem as a system of two coupled ODEs: $$d_t u(t) = v(t), \quad d_{tt} u(t) = d_t v(t) = f(t, u(t), v(t)).$$ Discretising with a time grid $t_j = j\tau$ we are looking for the approximations of $u$ and $v$ at these points. The approximations will be two sequences $\tilde{u} = (\tilde{u}^0, \tilde{u}^1, \ldots)$ and $\tilde{v} = (\tilde{v}^0, \tilde{v}^1, \ldots)$, where the initial conditions give you the first terms in the sequences: $\tilde{u}^0 = u^0 = u(0)$, $\tilde{v}^0 = v^0 = v(0)$. The Euler discretisations now take the form:

\begin{alignat}{3} \tilde{u}^{j+1} &= \tilde{u}^j + \tau\tilde{v}^j, &\quad& (\textrm{explicit Euler}), \\ \tilde{v}^{j+1} &= \tilde{v}^j + \tau f(t_j, \tilde{u}^j, \tilde{v}^j), &\quad& (\textrm{explicit Euler}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \tau \tilde{v}^{j+1}, &\quad& (\textrm{implicit Euler}), \\ \tilde{v}^{j+1} &= \tilde{v}^j + \tau f(t_{j+1}, \tilde{u}^{j+1}, \tilde{v}^{j+1}), &\quad& (\textrm{implicit Euler}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \frac{\tau}{2}( \tilde{v}^j + \tilde{v}^{j+1}), &\quad& (\textrm{Crank-Nicolson}), \\ \tilde{v}^{j+1} &= \tilde{v}^j + \frac{\tau}{2}\bigl( f(t_{j}, \tilde{u}^{j}, \tilde{v}^{j}) + f(t_{j+1}, \tilde{u}^{j+1}, \tilde{v}^{j+1}) \bigr), &\quad& (\textrm{Crank-Nicolson}), \\ \tilde{u}^{j+1} &= \tilde{u}^j + \frac{\tau}{2}\bigl( \tilde{v}^{j} + \tilde{v}^{j}+\tau f(t_{j}, \tilde{u}^{j},\tilde{v}^j) \bigr), &\quad& (\textrm{improved Euler}), \\ \tilde{v}^{j+1} &= \tilde{v}^j + \frac{\tau}{2}\bigl( f(t_{j}, \tilde{u}^{j}, \tilde{v}^j) + f(t_{j+1}, \tilde{u}^{j}+\tau \tilde{v}^j, \tilde{v}^{j}+\tau f(t_j,\tilde{u}^j, \tilde{v}^j)) \bigr), &\quad& (\textrm{improved Euler}). \end{alignat}

And this generalises to higher order PDEs in the same manner.

Your Problem

From the description of the implementation of your problem my guess is that you have something like this:

$$d_{tt} \vec{p}_i(t) = f_i(t, \{\vec{p}_k(t)\}_{k=1}^N, \{d_t\vec{p}_k(t)\}_{k=1}^N), \, 1\leq i \leq N.$$

Where $\vec{p}_i(t)$ is the position of the $i$-th particle at time $t$. Here I assume that all particles are allowed to interact, and thus the right-hand side contains all of those. So essentially you have $N\times d$ equations, where $N$ is the number of particles, and $d$ is the dimension of the particles. As before we can split this into two pieces $d_t \vec{p}_i(t) = \vec{v}_i(t)$ and $d_t \vec{v}_i(t) = f_i(t, \{\vec{p}_k(t)\}_{k=1}^N, \{\vec{v}_k(t)\}_{k=1}^N)$. For simplicity I will suppress the vectorial notation, but it should be clear that $\tilde{p}^j_k$ and $\tilde{v}^j_k$ are $d$-dimensional vectors in the following. Explicit Euler is given as: $$\tilde{p}^{j+1}_i = \tilde{p}^j_i + \tau \tilde{v}^j_i, \quad \tilde{v}^{j+1}_i = \tilde{v}^j_i + \tau f_i(t_j, \{\tilde{p}^j_k\}_{k=1}^N, \{\tilde{v}^j_k\}_{k=1}^N).$$

Here the $f$ subsumes your computations of the forces on the grid, by using the positions $\{\tilde{p}^j_k\}_{k=1}^N$ and velocities $\{\tilde{v}^j_k\}_{k=1}^N$ at time $t_j = j\tau$. Now consider the improved Euler scheme for the same:

\begin{align} \tilde{p}^{j+1/2}_i &= \tilde{p}^j_i + \tau \tilde{v}^j_i, \\ \tilde{v}^{j+1/2}_i &= \tilde{v}^j_i + \tau f_i(t_j, \{\tilde{p}^j_k\}_{k=1}^N, \{\tilde{v}^j_k\}_{k=1}^N), \\ \tilde{p}^{j+1}_i &= \tilde{p}^{j}_i + \frac{\tau}{2}\bigl(\tilde{v}^j_i + \tilde{v}^{j+1/2}_i\bigr), \\ \tilde{v}^{j+1}_i &= \tilde{v}^j_i + \frac{\tau}{2}\bigl( f_i(t_j, \{\tilde{p}^j_k\}_{k=1}^N, \{\tilde{v}^j_k\}_{k=1}^N)+f_i(t_{j+1}, \{\tilde{p}^{j+1/2}_k\}_{k=1}^N, \{\tilde{v}^{j+1/2}_k\}_{k=1}^N)\bigr). \end{align}

Note that now you have to evaluate $f$ two times at every step with different position and velocities. Once at the positions and velocities $\{\tilde{p}^j_k\}_{k=1}^N$, $\{\tilde{v}^j_k\}_{k=1}^N$; and once at $\{\tilde{p}^{j+1/2}_k\}_{k=1}^N$, $\{\tilde{v}^{j+1/2}_k\}_{k=1}^N$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.