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$
Why is this method for simulating a system of springs and masses unstable?
ode;numerical modelling;implicit methods
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.