mjbWorld

Update:
*Coriolis and centripetal forces are present in the object and these cause the change in the direction of W
*Modulus of W is (97%) constant
*Euler's equations hold for the local frame of the object (when I used the small numerical change in W for the current frame during the time step)
*A numerical simulation would be just as fast as a theoretical calculation for obtaining the new orientation (after a very small time step). Remember this is done in games (I am just curious about a theoretic solution :)) After several minutes of running the numerical simulation, the object begins to rotate about a principal axes, but this is not about the stable axes as it is caused by floating point errors in the numerical simulation that accumulate over time.

This is how the numerical simulaton works.

L is constant = IW {where W is the initial value and all terms are in the global frame}

O {this is the initial orientation}
O = O + dO = O + Wt {this is the orientation after the first time step with W}
W = (O * Io^-1 * O^-1 * L) {W is updated from the new O}
O = O + dO = O + Wt
W = (O * Io^-1 * O^-1 * L) {and so on}

This is my approach at solving Euler's equation:
Given,
d(w1)/dt = ((I2 - I3)/I1)w2w3
d(w2)/dt = ((I3 - I1)/I2)w3w1
d(w3)/dt = ((I1 - I2)/I3)w1w2

Let w1 = Ae^(Bix) where B is a polynomial that could be a possible solution with the differential equation. !It really is quite difficult!