# Why is this method for simulating a system of springs and masses unstable?

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$

• What you call the governing equation is not; it's simply the forward Euler discretization. – David Ketcheson Apr 8 '16 at 4:05
• 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. – WhiteTiger Apr 8 '16 at 7:03

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.
• 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}})$? – WhiteTiger Apr 7 '16 at 21:46
• 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. – Wolfgang Bangerth Apr 8 '16 at 18:57