I have a computer simulation system of bodies connected by springs, so their movement is governed by:

$x_{n+1} = x_n-\Delta tk(x_n-r)$

Where $r$ is the idea distance between every two bodies, and $\Delta t$ is the timestep.

Using the Backward Euler method, I end up with the equation:

$x_{n+1} = x_n-\Delta tk(x_{n+1}-r)$

which simplifies to:

$x_{n+1} = \frac{x_n+\Delta tkr}{1+\Delta tk}$

I then find that $\Delta tk(x_{n+1}-r)= x_n-x_{n+1} $ so I know the magnitude of my movement, and I can easily calculate a direction vector to apply to it.

I calculate such interactions between each two sets of bodies that are joined by a spring (I have about 28 bodies), and then add up the movement vectors for each body and apply the movement.

This is where I run into a problem. The numbers I get just keep growing and growing through the iterations until the reach infinity and my simulation explodes, unless I set a small timestep $\Delta t$. I thought Backward Euler should be able to avoid this problem and be unconditionally stable.

Is the problem on my implementation, or is there a limit to Backward Euler's stability too that is governed by the time step $\Delta t$

  • $\begingroup$ What you call the governing equation is not; it's simply the forward Euler discretization. $\endgroup$ Commented Apr 8, 2016 at 4:05
  • $\begingroup$ okay, I'm rusty on the definitions and wording, kinda new to this field. Hope that down-vote made you satisfied, you are officially smarter than someone who's been at it for a couple of months. $\endgroup$
    – WhiteTiger
    Commented Apr 8, 2016 at 7:03

1 Answer 1


The implicit Euler method is unconditionally stable alright, but what you are doing is not the implicit Euler method. Rather, what you do is compute where the particle would be at the end of the time step using only 2-particle interactions, compute the location update, and then sum these updates up for all particle interactions. But the implicit Euler method would compute all of these updates where you use as $x_{n+1}$ the final position using all interactions, not just the one interaction you are currently considering.

The method you describe is called "operator splitting", and it is not unconditionally stable even if each building block (each interaction separately) is discretized in a way that is unconditionally stable.

  • $\begingroup$ So for each particle I have to get the final position using an equation that has a term for each other particle? For example, if I have a particle that is attached to two other ones, my equation is $\frac{\Delta\vec{x}}{\Delta t} = -k(\vec{x_{01}}-\vec{r_{01}}) - k(\vec{x_{02}} - \vec{r_{02}})$? $\endgroup$
    – WhiteTiger
    Commented Apr 7, 2016 at 21:46
  • $\begingroup$ Yes, both $x_{01}$ and $x_{02}$ both need to be computed with the final position of the particles for this to be a completely implicit scheme. $\endgroup$ Commented Apr 8, 2016 at 18:57

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.