[how to?]Runge-Kutta 4 or Implicit Euler Integration

I’ve had suffers searching all around so I’m asking it here.

I have a basic damped spring model in my game. It uses Hooke’s law equation(F = -kx - bv). The stiffness in my case is very high(and so the angular frequency(iirc it’s sqrt(k/m)) is also very high, above 2 or so). This makes the explicit Euler unstable. As it is above 2, the semi-implicit Euler also becomes unstable. I’ve tried them both, indeed.
(explicit Euler)

v += (F/m) * dt
p += v * dt

I’ve heard of 2 methods that are both accurate and stable, however, they’re slower. However, it’s mostly told that they’re not as heavy as they’re stable. Sounds worth to try. Those are Runge-Kutta 4 and Implicit Euler. However, I have no idea how to implement them in my system. I need to do 2 steps - update velocity and update position while knowing the force. So how can I do this with each of those integrations?

Is all I saw there. I wish I actually understood half of what you wrote.
I wish I’d studied harder at math in school. :frowning:

But… A quick google check bring up this:
http://www.codeproject.com/Tips/792927/Fourth-Order-Runge-Kutta-Method-in-Python
I guess you will be doing it in 3 dimensions so you should probably use the second implementation there.

Hahaha…

Thanks for the link.
It may help. However, there is one problem: f(x,t) or f(y,t) is different approach, often used in curve generation. I must somehow manage to replace those with the physics integration all the time, but it’s not that simple.

I have an implementation at home I’ll dig it out

Nice! It’d be really, really great!:slight_smile:

By the way - there is just one problem for me at this case.

Where it does f(t + h/2, x + k1/2) I can’t understand what to place here is f(t, x) is equal to acc in my case. Is it even possible to do? Appearantly it is, cause some people has already done it somehow.

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 &gt; 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).

Ive got a headache

Yeah… My system even went unstable at 500 Hz and 1% real-time(50 kHz per game second) with explicit euler. That’s, of course, is because those springs are for softbody simulation of very, very, very stiff object(stiffness goes around 200-1000 kN/m) at very small, almost diverging springs(their length is around 5-15 cm). This adds the instability. The wile loop seems interesting. Also the fact that the function uses spring is interesting. However, in this case it seems that the spring must be integrated seperately from the rest? Or should I actually calculate the acceleration from force in k1, k2, k3, k4, than integrate velocity from it? Oh, and in the k2 it’s calculating like if it would have moved by half an Euler step? So that’s how it works. You are slowly removing my headache - thanks!:slight_smile:

My system even went unstable at 500 Hz and 1% real-time(50 kHz per game second)

500Hz should be plenty for just about anything. What are you simulating?
I really think you have a mistake somewhere if your system is unstable at 500Hz. I recently simulated some thermal stuff with a rise time of time of 2 seconds at 4Hz (so about 4 times niquist).
One thing to note is that high levels of spring-force and high levels of damping do go unstable when solved iteratively. This is simply the nature of them. If you have a spring of 1000000000kNm and a damping of 100000000kNms, you’ll find they are always unstable unless you simulate in nano-second timesteps. This is because when the spring is disturbed slightly, a huge response force appears, and a huge resulting displacement appears in a single timestep - before the (stabilizing) damping appears. But at that level of springiness, a rigid body will act exactly the same. Using a convergent method (such as backwards euler) may even be unstable, or at least take a long time to stabilize. I am not sure how RK4 would respond as I’ve never implemented it. Ask Agoose I guess?
Consider that you can scale everything down, which may help stabilize things: is springs of 20n/m, but make sure masses are 100,000 times lighter. This may stabilize things, it may not.

So what are you simulating?

Or should I actually calculate the acceleration from force in k1, k2, k3, k4, than integrate velocity from it? Oh, and in the k2 it’s calculating like if it would have moved by half an Euler step? So that’s how it works. You are slowly removing my headache - thanks!

If you work in vectors, all this comes out in the wash. Calculate your displacement and forces in vector form and then you can just sum and integrate them all at the same time. This should all be rolled into the function f(t, position).


It sounds to me like you are getting into the world of Finite Element Analysis. I did a course on that last year, and it was a pretty deadly course. The math was way over my head, but the techniques for simulating large arrays of simple systems may be useful. You’ll have to find someone else to help you with it though.

The basics of FEA is around connected nodes or volumes, depending on how you decide to make your finite elements. I find nodes conceptually simpler, so I’ll run with that.
Let’s run with 1D for now. If I have a bunch of nodes connected by springs:


Lets number them 0 - 10 from left to right. The force on the first node is simple, F = k (x_0 - x_1). The next node has two forces F = k (x_1 - x_2) - k (x_0 - x_1), and so on.
Turns out, if you grab all of these together and stick them in a vector equation you end up with the simple equation:
F = KX, whereby F is a vector [f0, f1, f2 … f9], X is a vector of positions [x0, x1, … x9], and k is, well, a 10x10 triagonal matrix:


k      -k     0  0   0 ....  
-k    2k   -k   0   0
0      -k   2k  -k  0
0      0      -k   .
.                       .
.                          2k    -k
                            -k    k

Assuming all spring constants are the same. (If you number your nodes right, a the matrix for a 2d array also looks really pretty)
Now we can calculate the forces on all the nodes at the same time. I highly recommend you sit down with paper and prove this to yourself.

And then there is a bunch of gorgeously hideous methods for actually solving these systems. I’m just going to leave it at that before I stuff up. The level of information you want to know is at university level (Euler and numerical integration is first year, FEA is third year), please take some engineering or mathematical modelling courses.

Attachments


Woo, tire simulation. Welcome to a world of hurt.
Why not just use bullet’s softbodies. What features is it missing? Have you considered implementing them in bullet rather than in blender? Just because you can implement it all yourself doesn’t mean you should. I will admit that built-in soft-bodies are pretty bad though.

I compare previous and current compression to get the compression speed(v = (compression - previous_compression) * dt). Than I apply damping.

And this one-frame lag between doing the springiness and damping means that at high stiffness values, things get unstable.

Your admit is perfect. They’re very bad for tires. They just can’t simulate tires. Hm… You say there’s damping lag? So I need to calculate damping implicitily? How can I do it?