Skip to main content
gamedev.net
🔒 Locked

RK4 with multiple objects

Started by daniel_i_l Aug 14, 2006 at 12:15 PM 6 replies 1.1k views
Original Post
daniel_i_l
daniel_i_l
I know that in RK4 you take 4 derivatives at different stages of the time interval. My question is, if I have 2 objects connected with a spring. Do I do all the stages on the first one, changing the position of the spring in each different stage only on the side of the first on and keep the end of the second one constant - and then reverse this for the second one? Or do I keep the second one constant and do the first stage with the first one, and then keep the first one constant and do the first stage on the second one, and keep on doing this for all the stages? Or do I do something else? Thanks
"We've all heard that a million monkeys banging on a million typewriters will eventually reproduce the entire works of Shakespeare. Now, thanks to the internet, we know this is not true." -- Professor Robert Silensky
Trap
Trap
You do it exactly the same way as with a single object. Just put everything you have into a single vector and redefine everything to work on vectors. (vectors meaning mathematical vectors or c++ vectors)
MrRowl
MrRowl
There's an example of how to do it in the biblical Baraff siggraph97 notes (Rigid Body unconstrained dynamics).

The short answer can also be found in this thread....
daniel_i_l
daniel_i_l
Thanks,
I don't understand how can I put all the bodies into one mathematical vector?
And if I put them in a C++ vector like this: vector Bodies;
then for every step I loop through the bodies? - for example instead of:
Deravitive GetDerivative(Body &body)
{
calculate deravitive on body;
}
I'd do:
Deravitive GetDerivative(vector &bodies)
{
for all bodies
calculate deravitive on current body;
}

and the same for all the other functions I have? Or did you mean something else? Could someone explain?
Thanks!
"We've all heard that a million monkeys banging on a million typewriters will eventually reproduce the entire works of Shakespeare. Now, thanks to the internet, we know this is not true." -- Professor Robert Silensky
Trap
Trap
That's already pretty close, conceptually you need something like this:

vector GetDerivative(vector &bodies)
{
for all bodies
calculate deravitive on current body;
}

At places where you use + and * in the normal algorithm you just apply them compenentwise.

This adds some overhead, if you know that some bodies don't interact with any other body it might be faster to integrate them seperately.
daniel_i_l
daniel_i_l
Ok, so basiclly I do everything together? I'll try this out and post what happens!
"We've all heard that a million monkeys banging on a million typewriters will eventually reproduce the entire works of Shakespeare. Now, thanks to the internet, we know this is not true." -- Professor Robert Silensky
daniel_i_l
daniel_i_l
Is this right?
Simulation.h
class CSimulation{private:    double dt;    	//get acceleration from current state	Vector2D Acceleration(RigidBody &bodies);    //get angular acceleration from state	double   AngularAcceleration(RigidBody &bodies);	//overload evaluation	vector<Derivative> Evaluate(vector<RigidBody> &initials);	//evaluate derivative	vector<Derivative> Evaluate(vector<RigidBody> &initials, double dt, vector<Derivative> &d);	//integrate with RK4	void Integrate(vector<RigidBody> &bodies, double dt);public:	CSimulation();	//do one iteration	void Update(vector<RigidBody> &currents, double DeltaTime);};


Simulation.cpp
CSimulation::CSimulation(){	dt = 0.01;}Vector2D CSimulation::Acceleration(RigidBody &body){	Vector2D TotalForce = body.TotalForce;	return  TotalForce/body.m_MyState.Mass;}double   CSimulation::AngularAcceleration(RigidBody &body){    double   Torque       = body.TotalTorque;	double   AngularAccel = Torque / body.m_MyState.MomentOfInertia;	return  AngularAccel;}vector<Derivative> CSimulation::Evaluate(vector<RigidBody> &initials){	Derivative output;	vector<Derivative> outputs;	   for(unsigned int i=0; i<initials.size(); ++i)   {	initials.ApplyForces();	output.dPosition = initials.m_MyState.Velocity;	output.dVelocity = Acceleration(initials);		output.dRotation = initials.m_MyState.AngularVelocity;	output.dAngularVelocity = AngularAcceleration(initials);	outputs.push_back(output);   }	return outputs;}vector<Derivative> CSimulation::Evaluate(vector<RigidBody> &initials, double dt, vector<Derivative> &d){	RigidBody  B;	vector<Derivative> outputs;	Derivative output;       for(unsigned int i=0; i<initials.size(); ++i)   {	B = initials;    	B.m_MyState.Position = initials.m_MyState.Position + d.dPosition * dt;	B.m_MyState.Velocity = initials.m_MyState.Velocity + d.dVelocity * dt;	B.m_MyState.Rotation        = initials.m_MyState.Rotation                                              + d.dRotation        * dt;	B.m_MyState.AngularVelocity = initials.m_MyState.AngularVelocity                                       + d.dAngularVelocity * dt;	B.m_MyState.Mass = initials.m_MyState.Mass;	B.m_MyState.MomentOfInertia = initials.m_MyState.MomentOfInertia;	B.ApplyForces();     	output.dPosition = B.m_MyState.Velocity;	output.dVelocity = Acceleration(B);	output.dRotation        = B.m_MyState.AngularVelocity;	output.dAngularVelocity = AngularAcceleration(B);	outputs.push_back(output);   }        return outputs;}void CSimulation::Integrate(vector<RigidBody> &bodies, double dt){	vector<Derivative> a = Evaluate(bodies);	vector<Derivative> b = Evaluate(bodies, dt*0.5, a);	vector<Derivative> c = Evaluate(bodies, dt*0.5, b);	vector<Derivative> d = Evaluate(bodies, dt,     c);    	//Euler	//Derivative b = a;	//Derivative c = a;	//Derivative d = a;    for(unsigned int i=0; i<bodies.size(); ++i)	{	Vector2D  dPosition = 1.0f/6.0f * (a.dPosition + 2.0f*(b.dPosition + c.dPosition) + d.dPosition);	Vector2D  dVelocity = 1.0f/6.0f * (a.dVelocity + 2.0f*(b.dVelocity + c.dVelocity) + d.dVelocity);		double    dRotation = 1.0f/6.0f        * (a.dRotation + 2.0f*                          (b.dRotation + c.dRotation) + d.dRotation);	double    dAngularVelocity = 1.0f/6.0f * (a.dAngularVelocity + 2.0f*     (b.dAngularVelocity + c.dAngularVelocity) + d.dAngularVelocity);			bodies.m_MyState.Position = bodies.m_MyState.Position + dPosition*dt;	bodies.m_MyState.Velocity = bodies.m_MyState.Velocity + dVelocity*dt;	bodies.m_MyState.Rotation        = bodies.m_MyState.Rotation + dRotation*dt;	bodies.m_MyState.AngularVelocity = bodies.m_MyState.AngularVelocity + dAngularVelocity*dt;	}}void CSimulation::Update(vector<RigidBody> &currents, double DeltaTime){  double accumulator = 0.0f;  vector<RigidBody> previous;  previous = currents;  if(DeltaTime > 0.25)	  DeltaTime = 0.25;  if(DeltaTime < dt)	  DeltaTime = dt;  accumulator += DeltaTime; while(accumulator>=dt)  {	  accumulator -= dt;	  previous = currents;	  Integrate(currents, dt);  } //RigidBody r; //r = Interpolate(previous, current, accumulator/dt);// current = r;}


It seem like it really slows thing down with all the extra loops...
Thanks.
"We've all heard that a million monkeys banging on a million typewriters will eventually reproduce the entire works of Shakespeare. Now, thanks to the internet, we know this is not true." -- Professor Robert Silensky
MrRowl
MrRowl
Quote:
Original post by daniel_i_l
Is this right?


I have no time to do more than glance at the code, but the best way is to test it yourself with a situation where you have a known analytical result (e.g. mass hanging from damped spring). Doesn't tell you if you have a _perfectly_ correct RK4 implemtation, but it will be a good indication that you have a correctly working integrator of some kind :) You can also compare against a Euler/RK2 integrator, and check that you get better accuracy.

Quote:

It seem like it really slows thing down with all the extra loops...


Of course, but you should find it's more stable/accurate than Euler scheme(s) for bigger time-steps, which compensates to some degree.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.

Sign in to reply to this topic.