Crossposting this from Reddit.

Basically, I have a custom 2d physics engine, written in AS3, for entirely academic purposes.

Previously, it consisted solely of jacobson advanced character physics style verlet integrated particles and constraints, but recently I added rigid body rotating disks.

I use verlet both for the integration of particles and the integration of the rotation of disks, so it looks like this:

For the particles' positions:

var temp:Vector2 = new Vector2(x, y); //Store old position

var tempv:Vector2 = new Vector2(v.x, v.y); //and old velocity



x = x * 2 - old.x + f.x; //integrate the velocity + accumulated force

y = y * 2 - old.y + f.y;



old.x = temp.x;

old.y = temp.y;



v.x = ((x - old.x) + tempv.x)/2; //averaging the velocity like this

v.y = ((y - old.y) + tempv.y)/2; //seems to help with stability



f.x = 0; //zero the force accumulators

f.y = 0;

For the rotation of uniform objects:

acc_force += -velocity * friction

velocity = velocity * 2 - old_velocity + (acc_force*radius)/inertia;

rotation += velocity

old_velocity = velocity;

acc_force = 0;

As you can see, I integrate the velocity of rotating bodies and use that to derive the position. This is because if I do things the other way, it was incredibly unstable. Averaging the position a la the velocity of fixed points also had a negative effect.

So basically, my issue is that rotation is a lot less stable than position. It tends to 'shudder' at low speed and it can explode quite easily.

Has anyone had similar experiences? Is there a better way to integrate rotation? Or is this more likely a numerical instability issue?

(acc_force*radius)/inertia

This makes little sense - the rasius value would affect the calculation of the inertia, but then you should not use it again during integration.

I wonder why you use inertia at all. You use unit mass for all particles, but then why individual inertia?

But this has probably nothing to do with your problem.

Actually for a object rotating at constant speed it's angle would never stop increasing until floating point precission issues creep in and small changes can't be handled anymore. (Although it's unlikely you get quickly to this point)

You could avoid that by keeping your rotation value always in the range -PI, +PI.

Of course using this you need to handle special cases because vel = new - old no longer works that simple, e.g:

actual = 1.5

prev = 1.4

vel = actual - prev // 0.1

next = actual + vel // 1.6 oops larger than PI/2

while (next > PI) next -= PI*2 // rotate 360 degrees back so we get the same orientation, but a clamped number with good precission

while (next < -PI) next += PI*2

Additionally a common trick is to set a max angular velocity per integration step as well so velocity can't go over 360 and you can replace while with if.

Edit: Oops, forgot to continue.

in the next step we have different values:

actual = -1.54

prev = 1.5

vel = actual - prev // -3.04 that's wrong, we want 0.1

// add your extarnal torque but make sure the added velocity is not too large, e.g. (-PI/2, PI/2)

extV = acc_force

if (extV > PI/2) extV = PI/2;

if (extV < -PI/2) extV = -PI/2;

vel += extV

if (vel > PI) vel -= PI*2 // do the same for velocity

if (vel < -PI) vel += PI*2 // so this line produces the correct value of 0.1 assuming external is zero

next = actual + vel

if (next > PI) next -= PI*2

if (next < -PI) next += PI*2

So assuming i did nothing wrong this should work.

Edit: Fixed terrible bugs

velocity = velocity * 2 - old_velocity + (acc_force*radius)/inertia;

That doesn't look like verlet integration to me. You are supposed to do that with the angle, not the angular velocity:

angle = 2 * angle - old_angle + angular_acceleration * delta_t * delta_t;

You may run into trouble with this formulation if you don't keep track of how many whole rotations you have performed. If you want to subtract some multiple of 2*pi from the angle because it's gotten too large, you have to do the same thing to old_angle, so the difference still makes sense.

You could also reformulate verlet integration so you'll remember (angle, angle-old_angle) instead of (angle, old_angle). This will make the code look more like Euler integration, but if you are careful you should get the same results:

angle_increment += angular_acceleration * dt * dt;
angle += angle_increment;

Then you can feel free to add or subtract multiples of 2*pi to angle.


NOTE: I have advocated many times in these forums the use of unit-length complex numbers instead of angles. If you were to do that, the code would look like this:

angle_increment += angular_acceleration * dt * dt;
attitude *= std::exp(std::complex<float>(0, angle_increment));
attitude /= std::abs(attitude); // make sure the length stays at 1

NOTE: I have advocated many times in these forums the use of unit-length complex numbers instead of angles. If you were to do that, the code would look like this: angle_increment += angular_acceleration * dt * dt; attitude *= exp(std::complex(0, angle_increment)); attitude.normalize();

Iteresting, this could help me to understand the use of complex numbers for quats or even better the usecases of complex numbers in general.

Could you link to an older thread or repeat some explantation for the 101th time? :)

You can have a look in Box2D. Erin doesn't use complex numbers, but more like a 2D stripped quaternion iirc.

Here goes a brief explanation.

If you think of what you are going to do with an angle, chances are you are going to use it to rotate something. In that case, you'll have to compute its cosine and its sine. But it turns out we can do all the operations we want to do with angles using a representation that consists of two numbers: The cosine and the sine of the angle.

So how do you rotate a point (x,y) by a rotation (C,S)? (x',y') = (C*x - S*y, S*x + C*y)

It turns out, that's exactly the formula for multiplying the complex number x+i*y times the complex number C+i*S. C+i*S is not an arbitrary complex number, because C^2+S^2=1. This means that we are representing rotations using complex numbers of length 1.

Common operations:

 typedef std::complex<float> C;

C normalize(C z) {
  return z / std::abs(z);
}

// The next two functions are so trivial that you don't even need them. Just multiply!
C apply_rotation_to_point(C rotation, C point) {
  return rotation * point; 
}

C compose_rotations(C r1, C r2) {
  return r1 * r2;
}

C inverse_rotation(C rotation) {
  return std::conj(rotation);
}

C rotation_towards(C original, C target) {
  return normalize(target / original);
}

C limited_rotation_towards(C original, C target, C max_rotation) {
  C rotation = rotation_towards(target, original);
  if (rotation.real() < max_rotation.real())
    return rotation.imag() > 0.0f ? max_rotation : conj(max_rotation);
  return rotation;
}

C integrate_angular_velocity_over_time(float angular_velocity, float delta_t) {
  return std::exp(C(0.0f, angular_velocity * delta_t)); // Check out Euler's formula
}

Notice how this code has no special cases to deal with angles going over 2*pi, which plague code that does similar things using angles. Also importantly, notice how you can do many operations without trigonometric functions.

Try to implement limited_rotation_towards using angles if you want to feel the power of this technique. :)

As an added bonus, this is a good entry into using quaternions for 3D rotations. Complex numbers for 2D rotations is much easier to understand, but there are many parallels.

That's really cool. I've added some more useful functions but i don't know what i'm doing here yet :)

I need to implement complex myself to get a better understanding, for now the exp is the biggest mystery to me.

Thanks, Alvaro!

C rotation_from_angle(float angle) {
    return C (cos(angle), sin(angle));
}

float angle_from_rotation(C rotation) {
    float angle = acos(rotation.real());
    return (rotation.imag() < 0.0f ? -angle : angle);
}

C lerp_rotation(C original, C target, float t) {
    C r = target*t + original*(1.0f-t);
    return normalize(r); // in comparision with quaternions this result is pretty bad for large angles
}

C slerp_rotation(C original, C target, float t) {
    C d = rotation_towards(original, target);
    float angle = angle_from_rotation(d) * t;
#if 0
    return original * rotation_from_angle(angle);
#else
    return original * std::exp(C(0.0f, angle)); // todo: how does this work and which method is faster?
#endif
}

It's a good idea to have that rotation_from_angle. You can implement it as you did, or with the exponential of an imaginary number, or using std::polar. Your method is probably fastest. Using exp is natural only if you are a mathematician and you are used to working with Lie groups, or if you are a big fan of Euler's formula. So just ignore it and go with your implementation.

I would prefer this:

float angle_from_rotation(C rotation) {
  return std::atan2(rotation.imag(), rotation.real());
}

// EDIT: An alternative, using Euler's formula
float angle_from_rotation(C rotation) {
  return std::log(rotation).imag();
}

I don't understand the comment "in comparision with quaternions this result is pretty bad for large angles". The results of lerp are just as bad with quaternions.

Here's a simpler implementation of slerp:

C slerp_rotation(C original, C target, float t) {
  return original * std::pow(target / original, t);
}

Ha ok - i'll stick at the trig which is familiar, although it boils down to the same math we would use when using 2x2 matrices or just cos and sin per angle.

(So rotations alone seem not the best example to illustrate the usefullness of complex numbers to a math noob like me)

Those log, exp, pow, polar things you do look very interesting, but maybe i should learn more basic things first like basic calculus. Life is too short to catch up anyways :)

I don't understand the comment "in comparision with quaternions this result is pretty bad for large angles". The results of lerp are just as bad with quaternions.

There is a big difference, check it yourself.

With the complex numbers method the result is is as expected: The same as interpolating two unit vectors and normalizing the result (same math). For angles close to 180 degrees, the error can go up to almost 90 degrees.

But with quaternions the error of lerping 2D rotations (both has the same axis) is never larger than maybe 5 degrees. The question is: Why does it work so well with quats?

I think the reason is related to different properties: Complex numbers store sin(angle) and cos(angle), but quats store axis*sin(angle / 2) and cos(angle / 2).

Probably this is also the reason why a quat can store one extre revolution (IIRC), but this complex number method can not.

Guys I can't thank you all enough for this. I have some reading to do, it seems.