Numerical Integration

There are more accurate integration method than Euler''s, which can handle larger timesteps, from Gaussian Quadrature to Fourth Order Runge-Kutta (for Differential eqns).

Usually, the frame rate fixes the timestep, not the other way around, because you measure the elapsed time between 2 frames and use that as your timestep.

Fruny sums up my thoughts pretty well. I would just add that if you have enough CPU, then you can do multiple Euler steps for each frame, so you can get the "accuracy" and still maintain frame rates.

Don''t know what you''re doing, but you should know that Euler integration is highly susceptable to numerical instability, and if you''re doing physics this can show up depending how you implement collision response, rigid body joints, or if you have springs.

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Unless I misunderstand the question you do it in a loop so perhaps delta_t is the frame duration divided by 10 and you do it 10 times. Integrating an arbitrary function:

float step_size = (end_time - start_time) / num_of_steps; float value = 0.0f; for (float curr_time=start_time; curr_time < end_time; curr_time += step_size) value += f(curr_time)*step_size; //or for midpoint value += f(curr_time+step_size/2)*step_size; //or trapazoid value += (f(curr_time+step_size)-f(curr_time))/2*step_size+f(curr_time)*step_size;

Actually, that all sounds pretty good ... I think you all three have answered my question ... thank you very much ...

Well, now, how exactly do I implement "Gaussian Quadrature" in my program? The closest information I found on it was having to do with multiple integrals ...

I HATE our ISP.asdfas.dfsajfd;lksafjaslkjdfsa, .

I typed in a very nice reply, hit post reply....web timeout. message lost...this happens all the time...ISP sucks.

Anyway, to summarize.

You DON'T want to use Gaussian Quadrature for time marching problems such as simulation of physics and motion, !!!!! . GQ is not designed for this type of problem, and I should have commented on this before after I read Fruny's post.

The Runge-Kutta method is appropriate for this type of problem, but you probably don't need 4th order. Try a 2nd order Runge-Kutta. It should be just fine, and faster than 4th order. (People always seem to like 4th order RK, but the extra precision may not be that useful. You do get some numerical stability benefits with the higher order RK, but those benefits are insignificant compared to changing to an implicit method.) Alternatively, use a Verlet integrator, but do NOT use the "leapfrog" version of Verlet. Leapfrog is just as problematic as explicit Euler. (You are using explicit Euler currently. There is another Euler technique, called implicit Euler. Its way better, probably even better than the Runge-Kutta, but its much harder to implement---implicit methods always are.)

I had some very nice links, but I didn't save them in my favorites and now, thanks to our lovely ISP, all that is gone. I'm too pissed to look them up again. Try my favorite search engine, www.google.com.

(It might not be the ISP. It might be our Corporate intranet server, which is problematic, .)

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Edited by - grhodes_at_work on December 18, 2001 11:21:04 AM

and what about Monte-Carlo integration? ...you may find something at http://mathworld.wolfram.com

Dave007
dave007@volny.cz

Thanks everyone (sorry about the ISP Graham, I would have loved to read your post). What I find most frustrating is that I have a site with the equations for 4th order Runge Kutta, yet no idea as to how I implement it. In LilBudyWizer''s code, what is "f?" If I was integrating velocity for position then "f" would obviously be the velocity, but how do I find the velocity at (time + step/2)? Seems that I would have to integrate acceleration all over again for velocity ... maybe this is correct ... ?

Basically, how do I transform my crappy way of integrating to this new magical way.

quote:

---

Original post by Dave007
and what about Monte-Carlo integration?

---

Monte-Carlo isn't really a type of integration. Its a type of simulation that deals with uncertainties and randomness. Basically, you loop over a calculation tens or hundreds or thousands of times, changing the calculation slightly with random numbers to find out how the randomness affects the solution. It would sort of wrap around a numerical integration, if you were interested in representing uncertainties. In the engineering world, it is used to design things (cars, airplanes) that are more robust, that react well even when things are random such as air turbulence, small manufacturing flaws, etc. (Actually, my company develops software called ProFES which implements Monte-Carlo simulation, . http://www.profes.com. web site is outdated).

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Edited by - grhodes_at_work on December 18, 2001 5:15:32 PM

i played around recently with a simple simulation of a ball bouncing on a plane (no squashing of ball, no loss of energy, constant gravity). I was using a dynamic collision detection routine, and also simple Eular integration. the height of the bounce reduced with time until the ball was still on the plane. I also did a similar thing, but where force downward grew inversely proportional to distance to plane. In this case bounce height fluctuated.
I was just wondering if these problems were to do with the Numerical Instabililty you mentioned? What would be the usual way to get round these problems in games?

thanks

Cameron

quote:

---

Original post by kamrann
i played around recently with a simple simulation of a ball bouncing on a plane (no squashing of ball, no loss of energy, constant gravity). I was using a dynamic collision detection routine, and also simple Eular integration. the height of the bounce reduced with time until the ball was still on the plane. I also did a similar thing, but where force downward grew inversely proportional to distance to plane. In this case bounce height fluctuated.
I was just wondering if these problems were to do with the Numerical Instabililty you mentioned? What would be the usual way to get round these problems in games?

---

Cameron,

Your unstable case is a perfect example of the inherent instability of the explicit (simple) Euler integration. The way you did your collision response in the latter case is exactly to implement the reaction force like a spring. Spring forces cause a sinusoidal response and the Euler method is unstable whenever there is a sinusoidal response. No matter how weak that response is. I imagine your first approach, the one that did not show the instabilities, calculated the after-collision velocities directly rather than apply a reaction force over a few time steps to simulate the collision process? The closed-form collision response approach is wonderful when you can do it, and I recommend it.

The reduction of the height of the ball even when you did not mean to reduce the total energy is also a result of the Euler method. All numerical methods (Euler, RK, Verlet, implicit, everything else) have truncation errors as a result of their being derived from truncated Taylor series expansions. When the truncation errors are nonzero (and they usually are), you will either artificially add/remove energy to/from the system. (This is known as "artificial dissipation".) Only in very special situations (known as "perfect shift" conditions in which all truncation errors are exactly zero) will the total energy truly remain constant. There is no real way to avoid this, but you can reduce the energy loss by choosing a different integrator. You don't want an energy gain since that means instability. Higher order methods result in lower magnitude energy fluctuations (when they are stable) but higher order methods also are less stable in general (due to reasons beyond the scope of this post!). You just want to do the best you can to minimize the energy loss.

There are a couple of ways to get around the instability issues. The first is to use closed form motion solutions, which are simple and even fairly trivial for simple projectiles and closed form momentum exchange for collisions. Its all the projectile and billiard ball collision stuff from high school physics. If you want to include continuous but nonconstant forces such as aerodynamic loading/drag/lift, A variation on the full closed form simulation is to do what you probably did for the first case---numerically simulate motion until a collision occurs and then calculate after-collision velocities directly using the closed form equations rather than numerically simulate the collision process itself.

The second approach, suitable for full numerical simulation of everything, is to choose a better integration method, such as the Runge-Kutta or Verlet integrators (but never leapfrog), or to go to implicit integrators (which are difficult to implement).

It turns out that many people actually implement explicit Euler wrong somehow. This is a result of the fact that you end up integrating two equations (equation for velocity and equation for position), and a failure to pay attention to details (like making sure the right-hand-side is evaluated entirely at the current time step and not the next time step). When people get it wrong by chance, what they actually end up implementing is perhaps closer to a verlet integrator which is more stable than the Euler integrator.

So long answer, but I hope it offers some useful insight.

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Edited by - grhodes_at_work on December 19, 2001 4:27:07 PM

lost me in places, but very interesting still.

Where you say about implementing Eular wrongly, are you saying that the code at the start of this thread is wrong and should really be:


position += velocity;
velocity += acceleration;


??

thankyou

Cameron

Edited by - kamrann on December 19, 2001 5:14:20 PM

Yeah, that part of your post confused me too Graham. Is what Kamrann posted the correct way to do Euler integration?

quote:

---

Original post by kamrann
lost me in places, but very interesting still.

Where you say about implementing Eular wrongly, are you saying that the code at the start of this thread is wrong and should really be:


position += velocity;
velocity += acceleration;


---

Yes, the code at the beginning of the article is wrong (e.g., it is not simple Euler). The implementation you provide is correct! The reason why is this. The explicit Euler method says that the right-hand-side of the equation uses the values from the current/previous time step and the left-hand-side is the next time step.

So in your first equation you're calculating position at time = t+dt and using velocity at time = t.

Same thing in your second equation. You're calculating velocity at time = t+dt and using acceleration at time = t.

Since you did explicit Euler correctly, your simulations are subject to the instabilities it can incur.

Here is something that is not true explicit Euler. It is the implementation from jag_oes.


velocity += acceleration;
position += velocity;


Looks almost the same but it is not. The velocity update is explicit Euler, but the position update is not explicit Euler. Why? Because the position update is calculated using velocity at time = t + dt, as calculated from the first equation. This subtle little change makes the second equation something different than explicit Euler. Happily, it may actually make the problem stable when true explicit Euler is not stable.

The thing that people may want to realize is that sometimes they think they've done explicit (or simple) Euler when they actually did this other integration. As a result, sometimes people defend explicit Euler as being better than it actually is---since their results are better than explicit Euler can produce.

Does that clarify my previous post a bit?

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.


Edited by - grhodes_at_work on December 21, 2001 10:58:45 AM

quote:

---

Original post by jag_oes
Yeah, that part of your post confused me too Graham. Is what Kamrann posted the correct way to do Euler integration?

---

Yes, Kamrann''s post is the correct way to do Euler integration. See my detailed response to his email. But I would say that your implementation is more stable than his. (Remember, he had discovered the Euler instability in his bouncing ball example.) You might want to experiment, though, just to see the difference.

Graham Rhodes
Senior Scientist
Applied Research Associates, Inc.

Clarified. Thanks Graham.

Cameron

as far as i know it is like this:

p(t+dt) ~= p(t) + v(t)dt + .5a(t)dt^2

so your code should be like this:

position += velocity*time + .5*acceleration*time*time;

oh, and you dont even need velocity at all, just put it out..

p(t+dt) ~= 2p(t) - p(t-dt) + a(t)dt^2

its even a little more accurate as far as i can see because you dont need velocity, wich you dont HAVE but integrate with error yourself..

we wanna play, not watch the pictures