Implementing semi-implicit backward Euler in a 1-DOF mass-spring system

Question

I have a simple (mass)-spring system wih two points which are connected with a spring. One point is fixed at a ceiling, so I want to calculate the position of the second point using a numerical method. So, basically I get the position of the second point and it's velocity, and want to know how these two value update after one timestep.

The following forces take effect on the point:

• Gravitational force, given by -g * m
• Spring force, given by k * (l - L) with k being the stiffness, l being the current length and L being the initial length
• Damping force, given by -d * v

Summed up, this leads to

• F = -g * m + k * (l - L)
• Fd = -d * v

Applying for example Explicit Euler, one can derive the following:

• newPos = oldPos + dt * oldVelocity
• newVelocity = oldVelocity + dt * (F + Fd) / m, using F = m * a.

However, I now want to use semi-implicit backward Euler, but can't exactly figure out where to derive the Jacobians from etc.

Answer 1

So it's probably easiest to see how this goes from considering the fully implicit method first, then going to the semi-implicit.

Implicit Euler would have (let's call these eqn (1)):

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*(newPos - L) - d*newVelocity)/m

For now let's just measure positions relative to L so we can get rid of that -kL term. Rearranging we end up with

(newPos, newVelocity) - dt * (newVelocity, k/m newPos - d/m newVelocity) = (oldPos, oldVelocity - g*dt)

and putting that into matrix form

((1,-dt),(k/m, 1 - d/m)).(newPos, newVelocity) = (oldPos, oldVelocity -g*dt)

Where you know everything in the matrix, and everything on the RHS, and you just need to solve for the vector (newPos, newVelocity). You can do this with any Ax=b solver (gaussian elimination by hand works in this simple case). But since you mention Jacobians, you're presumably looking to solve this with Newton-Raphson iteration or something similar.

In that case, you're essentially looking to solve the zeros of the equation

((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt) = 0

which is to say, f(newPos, newVelocity) = (0,0). You have a previous value to use as a starting guess, (oldPos, oldVelocity). Now you just want to iterate on

(x,v)ⁿ⁺¹ = (x,v)ⁿ + f((x,v)ⁿ)/f'((x,v)ⁿ)

until you get a sufficiently good answer. Here,

f(newPos,newVel) = ((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt)

and f'(newPos, newVel) is the Jacobian corresponding the matrix

((1,-dt),(k/m, 1-d/m))

Going through the process for semi-implicit is the same, but a little easier - not all of the RHS terms in eqns (1) are new quantities. The way it's usually done is

newPos = oldPos + dt * newVelocity
newVelocity = oldVelocity + dt * (-g * m + k*oldPos - d*newVelocity)/m

eg, the velocity depends on the old time value of the position, and the position on the new time value of the velocity. (This is very similar to "leapfrog" integration..) You should be able to work through the above steps pretty easilly with this slightly different set of equations. Basically, the k/m term in the matrix above drops away.

Answer 2

You wouldn't happen to be an student living in Zürich and taking the Physically based Simulation course at a renowned university? ;-)