OK, what's wrong with my Runge-Kutta implementation?

Started by Axter February 04, 2002 06:24 PM

10 comments, last by Axter 23 years ago

Author

122

February 04, 2002 06:24 PM

Hi, I''m trying to implement a 4th-order Runge-Kutta integrator, but I don’t seem to get the correct results. It seems to be relatively close, but I believe my implementation is wrong. I think this is because I don’t understand what the derivative function is supposed to do. I have searched with various search engines, but none of the examples I saw did what I wanted to do. I tried to modify a version that looked relatively close to what I want, but I’m sure I screwed something up. BTW, to test it, I use a loop with many iterations (1000) and simple Euler integration to compare the results. This small-step Euler integrator should give me a very accurate value, right? Basically, I have a simple formula like: force = 1 / (1 – depth) – 1. This formula describes the force applied to the object at a certain depth, where depth can be from 0 to 1. In the example below, I pass in the fTime, for the time step value, as well as the initial fDepth1 and the end fDepth2. Is this the correct parameters I should feed into the formula? Anyway, here is my implementation:

              

float ForceAtDepth(float fDepth1, float fDepth2)
{
	float fForce1 = (1.0f / (1.0f - fDepth1)) - 1.0f;
	float fForce2 = (1.0f / (1.0f - fDepth2)) - 1.0f;
	return (fForce1 + fForce2) * 0.5f;
}



float RungeKutta4(float fTime, float fDepth1, float fDepth2)
{
	float fK1, fK2, fK3, fK4;
	float fT1, fT2, fT3;
	float fForce;
	float fHalfDepth = (fDepth1 + fDepth2) * 0.5f;
	
	fK1 = fTime * ForceAtDepth(fDepth1, fDepth2);
	fT1 = fDepth1 + 0.5f * fK1;
	
	fK2 = fTime * ForceAtDepth(fHalfDepth, fT1);
	fT2 = fDepth1 + 0.5f * fK2;
	
	fK3 = fTime * ForceAtDepth(fHalfDepth, fT2);
	fT3 = fDepth1 + fK3;
	
	fK4 = fTime * ForceAtDepth(fDepth2, fT3);

	fForce = (fK1 + 2.0f * fK2 + 2.0f * fK3 + fK4) / 6.0f;
	return fForce;
}



float Euler(float fTime, float fDepth1, float fDepth2, int nSteps)
{
	float fCurDepth  = fDepth1;
	float fStepDepth = (fDepth2 - fDepth1) / nSteps;
	float fForce     = 0.0f;
	float fStepTime  = fTime / nSteps;

	for(int nLoop = 0; nLoop < nSteps; nLoop++)
	{
		fForce    += ForceAtDepth(fCurDepth, fCurDepth + fStepDepth);
		fCurDepth += fStepDepth;
	}

	fForce *= fStepTime;

	return fForce;
}

Thanks in advance, SS

LilBudyWizer

491

February 04, 2002 11:01 PM

Is your depth a linear function of time? If it isn''t then there is a problem. Even if it is that seems a confusing way to go about it. Either pass a depth or time interval.

Keys to success: Ability, ambition and opportunity.

Axter

Author

122

February 04, 2002 11:16 PM

quote:
Original post by LilBudyWizer
Is your depth a linear function of time? If it isn''t then there is a problem. Even if it is that seems a confusing way to go about it. Either pass a depth or time interval.

Yes, I guess the depth is a linear function of time, because the object moves in a straight line during the timestep. But the force is not linear wrt the depth (as in the ForceAtDepth() function).

The problem is I can''t figure out what I need to change to make it work. I know I need to change something, but all the different implementations I was able to find are different enough that I can''t figure out what I need to change. This is why I asked. I thought my intention would be clear enough so that the smart math guys can tell me where my problem is.

Basically I have a collision with a non-linear depth-to-force ratio. I want to implement a Runge-Kutta integrator to calculate the force that needs to be applied to the object during a specified time-step. I''m sure it''s a pretty basic case.

Thanks,
SS

LilBudyWizer

491

February 05, 2002 03:04 AM

I may be wrong, but I think you are confusing integration of a differential equation with integration of a non-differential equation. You Euler method correctly estimates the area under a force versus time curve. You average the values of the function at the start and end of an interval and multiply by the width of an interval. So you are estimating it with trapazoids. Your Runge-Kutta method is doing an entirely differant thing. It is finding the value of a function y(x) where y''(x,y)=(x/(1-x)+y/(1-y))/2. You have an initial value condition of y(fDepth1)=fDepth2. That is entirely differant than your Euler function which is finding the integral from fDepth1 to fDepth2 of x/(1-x). It may be the same function in the program, but how it is interpetted mathematically is entirely differant. I am quite surprised that they give answers anywhere close to being the same. They are both correct, but they do entirely differant things.

Keys to success: Ability, ambition and opportunity.

Timkin

864

February 05, 2002 09:07 PM

Axter,

I suspect you are going about this problem in a most difficult manner!

Let''s consider an alternative description of your penetration problem. If you decide to continue with your method this will help guide you.

One common method to describe the interpenetration of objects is to imagine that the force applied to an object is proportional to the depth of penetration and in the opposite direction of motion. This is equivalent to solving the problem of a mass on a spring, or Simple Harmonic Motion.

In the SHM formulation the restoration force applied to the mass is given by f=-k x, where k is a positive spring constant. Using Newton''s Second Law, f=ma, we can write

   d2xm  ---   +  k  x = 0    dt2

which has a solution of

x(t) = a cos(w t+b )

where a is the amplitude of the oscillation, w (=k /m) is the frequency of oscillation and b the phase angle (non-zero if the object has a non-zero velocity at t=0).

Here x(t) describes the depth of the mass relative to an equilibrium point (in your problem the contact surface).

Now, you have provided a different formula for force vs depth (x) of the form f = ¹/_(1-x)-1 = ^x/_(1-x). Following the reasoning above we can then write

   d2x        xm  ---   -   ---   = 0    dt2      (1-x)

There may be a closed form (analytic) solution to this ODE but I''m not going to spend the time working it out... you can if you like. Alternatively you can integrate this equation numerically.

My suggestion is stick with the SHM model of penetration... it''s tried and trusted by many!

Cheers,

Timkin

chrisflatley

122

February 06, 2002 05:54 PM

Finally got around to not being an anonymous poster... :-)

My version of RK4, from Baraff's Differential Equation paper taken from his Siggraph97 course say that the formulas should be slightly different. (search google for Baraff's papers you'll find them easily, if you don't already have them)

First you need an initial time and a step time. You've called your step time fTime, but since you function is independent of time it doesn't need to be used. This maybe were you are getting confused.

So you get (I'll try and be consist with your code):

        {// blah blahfK1 = fTime * ForceAtDepth(fDepth1, fDepth2);	fK2 = fTime * ForceAtDepth(fDepth1 + fK1/2.0f, fDepth2);			fK3 = fTime * ForceAtDepth(fDepth1 + fK2/2.0f, fDepth2);		fK4 = fTime * ForceAtDepth(fDepth1 + fK3, fDepth2);	fForce = (fK1 + 2.0f * fK2 + 2.0f * fK3 + fK4) / 6.0f;// blah blah}

I think this mimicks the effect of your euler step but I don't think it's what you want.

You actually need to recalcuate the depth at the different times. So you need to actually work out the effect of the previous calculated force at the intermediate times. This gets a whole lot more messy.

eg what is fDepth2? The depth at the end of the timestep? You don't know that yet - that what the integrator is for.

I'd put some code but Baraff does a much better job. :-)

Chris

Edited by - chrisflatley on February 6, 2002 7:09:40 PM

chrisflatley

122

February 06, 2002 06:15 PM

I don''t think that calculating the forces required to prevent the collision is posibile. I''m not sure - again I point you to Baraff papers, entitled things like Fast Collision Force Computation and so on. You need to solve some serious equation to do this.

Spring style constraints are easy but not very stable - RK4 will help but that k constant will give you nightmares.

A good reference for all of these is www.q12.org where you can download source code (I think) for a real physics engine, with RK4. (Which you could use in your app)

Chris

Axter

Author

122

February 06, 2002 06:15 PM

Chris,

I tried your suggestion, but it doesn't work. I think part of the problem is that in the following lines of code:

    fK1 = fTime * ForceAtDepth(fDepth1, fDepth2);fK2 = fTime * ForceAtDepth(fDepth1 + fK1/2.0f, fDepth2);

...the function converts the depth to a force value (which could have a totally different range), and feeds this value back into the function again as depth in the second line. For instance, the ForceAtDepth function expects values between 0 and 1, but the force at depth of 1 is much larger than 1.

Any more ideas?

Thanks again for the help with this...

SS

Edited by - Axter on February 6, 2002 7:16:36 PM

LilBudyWizer

491

February 07, 2002 01:07 AM

I tried out a modified version of your routine. Runge-Kutta actually does pretty good. Just for comparision I used the left, right, midpoint, trapazoid and simpson''s rules to integrate x/(1-x) from .25 to .75. Since the anti-derivative is easy to come up with, i.e. -(ln(x)+x), I could find how many iterations it took each to get within a given tolerance. To get within .001 of the actual value took 667, 668, 13, 18 and 749 for each in the order listed above. Using the weighted average of the trapazoid rule and midpoint rule with the midpoint rule getting twice the weight got it within tolerance within 3 iterations. To check that it wasn''t a fluke I check the midpoint, trapazoid, and runge-kutta at .0001 yielding (39,55,5) and .00001 yeilding (122,173,8). At .0000001 it took Runge-Kutta only 24 times and at .000000001 only 76 times. That is blowing the other five completely out of the water.

Well, anyway, you just need to call a function returning x/(1-x) instead of the ForceAtDepth function. You don''t need the fT# or second parameters at all assuming you want to integrate x/(1-x).

Keys to success: Ability, ambition and opportunity.

Axter

Author

122

February 07, 2002 01:35 AM

Thanks for looking at this. Yes, Runge-Kutta seems pretty good.

Anyway, I''ll try out what you suggested, and see how that works, thanks!

SS

OK, what's wrong with my Runge-Kutta implementation?

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

OK, what's wrong with my Runge-Kutta implementation?

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines