Integrators

Author

122

April 28, 2002 01:05 PM

Hey, I am not sure that my Runge-kutta4 and Midpoint integrators are correct! I have a box which "flies" with speed of 10 and I always apply to it a force vector (perpendicular to speed vector) which is m*v*v/r Newton. So it goes on a circle track. But while Euler integrator increases the radius (because of bad precision), the midpoint and rk4 integrators decreases it and it becomes 0 rather fastly. Is this the normal functioning of these integrators? Other: could you point me out how runge-kutta2 works? In google I just found foreign(italian IIRC) homepages which I can''t translate

That is my actual code: Is this the right way? ( t0-t1 is the time, RigidBody::dydt(time) is the dydt function from David Baraff''s document; it calculates the derivatives)

  
void OdeSolver::Solve(SCALAR t0, SCALAR t1, RigidBody *obj)
{
	SCALAR step_mul = OdeSolver::StepSize / (t1-t0);
	if (step_mul>1) step_mul=1;
	switch (Func) {
		case ODE_EULER: {
			for (SCALAR i=t0; i<t1; i += OdeSolver::StepSize) {
				*obj += obj->dydt(i) * step_mul;
			}
			break;
		}
		case ODE_MIDPOINT: {
			for (SCALAR i=t0; i<t1; i += OdeSolver::StepSize) {
				RigidBody m = obj->dydt(i) * step_mul/2;
				*obj += (*obj + m).dydt(i+OdeSolver::StepSize/2) * step_mul/2;
			}
			break;
		}
		case ODE_RUNGEKUTTA2: break;
		case ODE_RUNGEKUTTA4: {
			for (SCALAR i=t0; i<t1; i += OdeSolver::StepSize) {
				RigidBody k1,k2,k3,k4;
				k1 = obj->dydt(i) * step_mul;
				k2 = (*obj + k1).dydt(i+OdeSolver::StepSize/2) * step_mul;
				k3 = (*obj + k2).dydt(i+OdeSolver::StepSize/2) * step_mul;
				k4 = (*obj + k3).dydt(i+OdeSolver::StepSize) * step_mul;
				*obj += (k1 + k2*2 + k3*2 + k4)/6;
			}
			break;
		}
	}
}

Thanks in advance for your help! -- tSG --

-- tSG --

gumby

122

April 28, 2002 02:14 PM

Have you tried different step sizes? Maybe you
need to set it adaptively....

tSG

Author

122

April 28, 2002 03:05 PM

I have tried some versions, but haven''t tried adaptive step sizes... But 0.001 (OdeSolver::StepSize) should be enough I think. If I change it to 0.0001, it runs ~0.004 fps, so I can''t lower it with a great amount...

But maybe that is how it should work, do you have any experiences?

And thanks gumby for the suggestion, I forgot to tell you the stepsize...!

-- tSG --

-- tSG --

tSG

Author

122

April 28, 2002 04:55 PM

And I also forgot to ask you about an another interesting issue: why does Baraff calculate the Rotation quaternion on such a strange way in dydt? dydt should calculate the _change_ of rotation and he also uses the current Rotation. For me this would be more unambiguous:
s.Rotation = Quaternion(Vector3(omega), omega.Length()); <- this is a quaternion which rotates around the vector omega by angle |omega|
(instead of s.Rotation = (Quaternion ( omega.x, omega.y, omega.z, 0 )*Rotation) * 0.5; )

Because the ode adds together the current rotation and the one which dydt results... What am I thinking incorrectly? (the document can't be wrong, but for me only my code works, so I misunderstood something...)

Thx and cheers

-- tSG --

[edited by - tSG on April 28, 2002 5:56:29 PM]

[edited by - tSG on April 30, 2002 1:09:48 PM]

-- tSG --

jobless

122

May 01, 2002 01:43 AM

The derivative of an orientation quaternion is:
dQ/dq = 1/2 * omega * q;
where omega is the quaternion (omega.x, omega.y, omega.z, 0)

It''s just a formula; I can''t explain it. If it doesn''t work, make sure you normalize q after the euler step:

q += dq * (t1-t0);
Normalize(q);

tSG

Author

122

May 01, 2002 09:59 AM

quote:
Original post by jobless
The derivative of an orientation quaternion is:
dQ/dq = 1/2 * omega * q;
where omega is the quaternion (omega.x, omega.y, omega.z, 0)

It''s just a formula; I can''t explain it. If it doesn''t work, make sure you normalize q after the euler step:

q += dq * (t1-t0);
Normalize(q);

Thanks for the reply! I also don''t know how this formula comes out, but in theory |omega| says how much the object rotates in a second and that is the value what we are looking for in dydt, not?

Although I normalize q after the euler/rk4 step it doesn''t work

I am really perplexed why it doesn''t work! It just simply doesn''t rotate the object.

That is what I am doing at the moment:
The angular momentum is 1/3 so it should rotate one round for 3 seconds. I set up the inertia tensor correctly for the 10x10x10 box.

In dydt:
Matrix4 R = Rotation;
Iinv = R * IbodyInv * R.Transpose();
omega = Iinv * Momentum;
...
Return s.Rotation = Quaternion(omega.x,omega.y,omega.z,0)*0.5*Rotation;

After returning the derivatives I multiply the current Rotation with the returned s.Rotation with the quaternion multiplication.
Is this the correct way how I should do that? (the matrix and quaternion operations are good I suppose, because they worked in a lot of situations so far...)

Thanks!

-- tSG --

-- tSG --

jobless

122

May 01, 2002 12:34 PM

Not sure I can help you further... When you get the derivative of your orientation, you multiply it by the time step and ADD it to the old orientation. Also, angular VELOCITY determines how fast an object rotates, not angular momentum. You can try to crank up the torque and see what happens.

tSG

Author

122

May 01, 2002 02:04 PM

quote:
Original post by jobless
Not sure I can help you further... When you get the derivative of your orientation, you multiply it by the time step and ADD it to the old orientation. Also, angular VELOCITY determines how fast an object rotates, not angular momentum. You can try to crank up the torque and see what happens.

Yes, but if I want to sum the effect of two quaternions I have to use the multiplication function.
Indeed, angular velocity is which I wanted to say, so omega

What does "crank up" mean? (sorry, my mother language is not english, and the dictionary didn't help)

Angular velocity is similar to linear velocity. With euler integrator the following can be written:
Position += linvelocity*time;
Why isn't Rotation += angvelocity*time? (of course here += means the quaternion multiplication, which adds together the rotates)

-- tSG --

[edited by - tSG on May 1, 2002 3:05:37 PM]

[edited by - tSG on May 1, 2002 3:06:50 PM]

-- tSG --

jobless

122

May 01, 2002 04:24 PM

I believe in the Baraff paper, the derivative and previous orientation quaternions are added. That is, literally,
Qnew = Qold + dQ =
(Qold.x + dQ.x, Qold.y + dQ.y, Qold.z + dQ.z, Qold.w + dQ.w)

Assuming dQ has been multiplied by the time step. After the addition you normalize Qnew. You don''t multiply dQ because it''s not a rotation, it''s a derivative.

What I meant about the torque was, since you said the object wasn''t rotating, maybe you need to increase ("crank up") the torque in order to see the rotation.

tSG

Author

122

May 01, 2002 06:01 PM

Hmm, I see now a difference!
I use the following code for "summing up" the current rotation and my "derivative":

  Quaternion Quaternion::operator * (const Quaternion &q) const {	Quaternion n;	n.x = w * q.x + x * q.w + y * q.z - z * q.y;	n.y = w * q.y + y * q.w + z * q.x - x * q.z;	n.z = w * q.z + z * q.w + x * q.y - y * q.x;	n.w = w * q.w - x * q.x - y * q.y - z * q.z;	n.Normalize();	return n;}

That is completely other than what you have written.
I''ve tried the version which you noticed but it works really strange. It seems that it rotates the cube by k* Pi/2 radians very fastly, because I can see the cube staying in the original position, but the color of its faces are added together, so it is almost white... (except the top and bottom points which are positioned on the rotation axis; they have the appropriate color)
Are there quaternion gurus? What''s your opinion what is happening here?

Now I see the "crank up" expression. Actually the torque is 0 because I set the angular momentum at the beginning and it doesn''t change while running(because torque=0), so it should rotate with a constant angular velocity. (the "joke" is, that if I use my version, it rotates with the good speed...ie. 3sec is 360 degree)

Thanks for your help!

-- tSG --

-- tSG --

Integrators

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Integrators

This topic is closed to new replies.

Popular Topics

Recommended Tutorials

Reticulating splines