Yes, simulating physics using any iterative method requires a small enough timestep otherwise it goes unstable - as you have discovered.
Bullet physics solves this using substeps, and things like constraints can be set to simulate many times per frame, however, performing physics operations in python that need high accuracy, once per frame simply isn’t enough.
However, I would think that if you’re hitting instability at 60FPS, there is likely something else wrong. I don’t think you’re doing things at massive speeds, or massively changing forces/accelerations. But then, if your spring is really high valued, you will be. Have you looked at the resonant frequency of a car recently? If you push down on one side, it takes maybe a second to do a full osscilation. 60FPS of physics should be plenty for car suspension.
For those who don’t understand what’s going on here. (Note that F0 is the force at time zero, F1 is the force at time one etc.)
Eulers method:
Iteration <b>1</b>:
Calculate force F<b>1</b>
F1 = f(t=1, position_at_time_0) # Some function of time and position last frame
Apply force F<b>1</b> to an object by:
accel = (F1/m)
velocity = accel*dt
position_at_time_1 = velocity * dt
Increment dt
Iteration 2:
Calculate force F2:
F2 = f(t=2, position)
Apply force F2 to an object:
accel = (F2/m)
....
Backwards Euler’s method:
Iteration <b>1</b>:
Calculate force F<b>2</b> by:
F2 = f(t=<b>2</b>, position_at_time_1) # Oh help, we need the position from the future to calculate the force now!
Apply force F<b>2</b> to an object:
accel = (F2/m)
velocity = accel*dt
position_at_time_1 = velocity * dt
Increment dt
Iteration <b>2</b>:
Calculate force F<b>2</b> by:
F2 = f(t=<b>2</b>, position_at_time_2)
Apply force F<b>2</b> to an object:
accel = (F2/m)
velocity = accel*dt
position_at_time_2 = velocity * dt
Increment dt
Notice that the backwards eulers method required information from the future, ie, the position of the current timestep.
So how do we calculate this position_at_time_2 before we have the force at time 2 which we need to calculate it?
Well, we hope they converge from the previous location, and run it iteratively:
For a single timestep:
while converged > 0.01:
Calculate_force:
F = f(t=1, guess_at_position)
new_guess_at_position = F/m*dt*dt
converged = abs(guess_at_position - new_guess_at_position)
guess_at_position = new_guess_at_position #For the next iteration, use a better guess
The idea is that every time that loop runs, the position will tend towards a certain location. This is ‘fixed point iteration’
However, let’s stop and consider your situation. You have a spring. You calculate it’s force discretely, as a function of position. So where does this go? It goes into the function f(t, position). You’ll notice there are a few caveats about implicit euler:
- Will it converges so it doesn’t infinite loop?
- How long it will take to converge? (ie how much computation it will take)
- How do you do the first guess_at_position anyway?
Well, to answer those there are many courses on system identification, modelling and so on at your local university. I haven’t taken any studying them in-depth, so I can’t really help. (Or if I have, I’ve forgotten them already)
I’m not familiar with the Runge-Kutta method, though I’m sure I’ve come across the name before, but here goes, after some quick research
The wikipedia article tells me to start with an initial value problem, which is:
velocity = function(position)
velocity = -k * position / m * dt #For a simple spring, note the single euler integration though
Fo solve for position, you do:
def f(position):
return -k * position / m * dt #For a simple spring
k1 = f(position_last_frame)
k2 = f(position_last_frame + k1 * dt / 2)
k3 = f(position_last_frame + k2 * dt / 2)
k4 = f(position_last_frame + dt*k3)
position = position_last_frame + (k1 + 2k2 + 2k3 + k4)*dt / 6
But we still have a single euler integration to convert from acceleration to velocity, so we have to clobber that somehow:
We have the second order:
Utt = -k * U #Spring
Ut = Z #One first order equation....
Zt = -k * U #Hey look, another first order equation, because we seperated out the second integration
So either we can make do with one euler integration, or we can apply the RK4 method twice (I think).
