Group Abstract Group Abstract

Message Boards Message Boards

How can I avoid complex components in a quaternion? (Seems like a bug)

GROUPS:
In the Quaternion package tutorial, it says:

A four-dimensional analog of de Moivre's theorem is used for calculating powers of quaternions
In[15]:=
Quaternion[1, 2, 0, 1]^2.5
Out[15]= Quaternion[-9.0604, 2.20741, 0., 1.1037]

That's fine, but if I take the 2.5 root again, I get

In[16]:= Quaternion[-9.0604, 2.20741, 0., 1.1037]^2.5
Out[16]:=Quaternion[
6.51112*10^-14 + 212.669 I, -4.56519*10^-14 -
149.11 I, -1.39769*10^-30 - 4.56519*10^-15 I, -2.2826*10^-14 -
74.5552 I]

The quaternion components are now complex (e.g. the first component above is 6.51112*10^-14 + 212.669 I). This seems to occur when the real part of the quaternion is negative. 

This seems like a bug. The output ought to be another well formed quaternion.

I thought I could simplify this by going:

Quaternion[ 6.51112*10^-14 + 212.669 I, -4.56519*10^-14 - 149.11 I, -1.39769*10^-30 - 4.56519*10^-15 I, -2.2826*10^-14 - 74.5552 I]
= (6.51112*10^-14 + 212.669 I) + (-4.56519*10^-14 - 149.11 I)*I + (-1.39769*10^-30 - 4.56519*10^-15 I)*J + (-2.2826*10^-14 - 74.5552 I)*K
=  6.51112*10^-14 + 212.669 I -4.56519*10^-14 I - 149.11 I^2  - 1.39769*10^-30 J - 4.56519*10^-15 IJ  - 2.2826*10^-14 K - 74.5552 IK
= 6.51112*10^-14 + 212.669 I -4.56519*10^-14 I + 149.11  - 1.39769*10^-30 J - 4.56519*10^-15 K - 2.2826*10^-14 K + 74.5552 J
= (6.51112*10^-14 + 149.11) + (212.669 - 4.56519*10^-14) I  + (74.5552 - 1.39769*10^-30) J - (4.56519*10^-15 + 2.2826*10^-14) K
= Quaternion[149.11, 212.669, 74.5552, -2.73911*10^-14]

But this is incorrect as Quaternion[149.11, 212.669, 74.5552, -2.73911*10^-14]^0.4 <> Quaternion[-9.0604, 2.20741, 0., 1.1037]

What is going on here? How can I either avoid getting complex components when I take quaternion roots or eliminate the complex components when they appear?
POSTED BY: Marc Widdowson
Answer
11 months ago
Hi Marc, 

this is a bug
In[1]:= Needs["Quaternions`"]

In[2]:= Quaternion[1, 2, 0, 1]^2.5
Out[2]= Quaternion[-9.0604, 2.20741, 0., 1.1037]

In[3]:= QuaternionQ[(Quaternion[1, 2, 0, 1]^(2.5))^(2.5)]
Out[3]= False
but, normally, if one speaks about the root of a number, one  has   to mean the principal root (cmp. http://mathworld.wolfram.com/PrincipalSquareRoot.html).
This seems to be implemented in the Quaternions package
In[4]:= (Quaternion[1, 2, 0, 1]^(1/2))^(5) // N
Out[4]= Quaternion[-9.0604, 2.20741, 0., 1.1037]

In[5]:= (Quaternion[1, 2, 0, 1]^5)^(1/2) // N
Out[5]= Quaternion[9.0604, -2.20741, 0., -1.1037]
try
Quaternion[1,2,0,1]^(1/n) // N
for integers n. Therefore one has simply to modify the rules for Power in the package
Quaternion[a,b,c,d]^(p/q) => (Quaternion[a,b,c,d]^p)^(1/q)
adding to the original ones
 (* *************** Rules for Power *************** *)
 (* Power is essentially de Moivre's theorem for quaternions. *)
 
 Quaternion /:
     Power[a:Quaternion[__?ScalarQ],0]:= 1
 
 Quaternion /:
     Power[a:Quaternion[__?ScalarQ],n_]:= Power[1/a, -n] /; n < 0 && n != -1
 
Quaternion /:
    Power[a:Quaternion[__?ScalarQ],-1]:= Conjugate[a] ** 1/Norm[a]

Quaternion /:
    Power[a:Quaternion[__?ScalarQ], n_]:=
    Module[{angle, pqradius},
    angle = toAngle[a];
    pqradius = angle[[2]]^n Sin[n angle[[3]]];
    Quaternion[
        angle[[2]]^n Cos[n angle[[3]]],
        Cos[angle[[4]]] pqradius,
        Sin[angle[[4]]] Cos[angle[[5]]] pqradius,
        Sin[angle[[4]]] Sin[angle[[5]]] pqradius
    ]  // QSimplify
    ] /; ScalarQ[n] && n > 0
one more rule and modifying another one
 (* *************** Rules for Power *************** *)
 (* Power is essentially de Moivre's theorem for quaternions. *)
 
 Quaternion /:
     Power[a:Quaternion[__?ScalarQ],0]:= 1
 
 Quaternion /:
     Power[a:Quaternion[__?ScalarQ],n_]:= Power[1/a, -n] /; n < 0 && n != -1
 
Quaternion /:
    Power[a:Quaternion[__?ScalarQ],-1]:= Conjugate[a] ** 1/Norm[a]

Quaternion /:
    Power[a:Quaternion[__?ScalarQ], n_]:=
    Module[{angle, pqradius, r},
    angle = toAngle[a];
    pqradius = angle[[2]]^n Sin[n angle[[3]]];
    Quaternion[
        angle[[2]]^n Cos[n angle[[3]]],
        Cos[angle[[4]]] pqradius,
        Sin[angle[[4]]] Cos[angle[[5]]] pqradius,
        Sin[angle[[4]]] Sin[angle[[5]]] pqradius
    ]  // QSimplify
    ] /; ScalarQ[n] && (Numerator[Rationalize[n]]==1 || Denominator[Rationalize[n]]==1)

Quaternion /:
    Power[a:Quaternion[__?ScalarQ], n_]:=Block[{r=Rationalize[n], y},
        y=If[Element[n,Rationals],1,1.,1.];
        Power[Power[a,Numerator[r]],y/Denominator[r]]
    ] /; ScalarQ[n] && n > 0
finally (after restarting Mma and reloading the package) arriving at
 In[1]:= Needs["Quaternions`"]
 
 In[2]:= Quaternion[1, 2, 0, 1]^2.5
 Out[2]= Quaternion[9.0604, -2.20741, 0., -1.1037]
 
 In[3]:= (Quaternion[1, 2, 0, 1]^2.5)^(2.5)
 Out[3]= Quaternion[212.669, -149.11, 4.56519*10^-15, -74.5552]
 
 In[6]:= QuaternionQ[(Quaternion[1, 2, 0, 1]^2.5)^(2.5)]
Out[6]= True
be aware that the arithmetics of quaternions is not too similar to the arithmetics of complex numbers. For example there is a continuum of square roots
of -1 in the skew field of quaternions showing another weakness of the Quaternions package
In[4]:= Sqrt[Quaternion[-1, 0, 0, 0]]
During evaluation of In[4]:= Power::infy: Infinite expression 1/0 encountered. >>
During evaluation of In[4]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>
Out[4]= Quaternion[Sqrt[\[Pi]/2], 0, Indeterminate, Indeterminate]
Compare this to the well-known 
In[5]:= Sqrt[-1]
Out[5]= I
again, note that it says I, not -I  ...

Regards
Udo.
POSTED BY: Udo Krause
Answer
11 months ago
Hi again, 

my proposal does not work in some cases. The strange thing with the original Quaternions package is that one can do things like that:
 In[1]:= Needs["Quaternions`"]
 
 In[2]:= (* Orig *) Quaternion[-1, 2, -7, 1]^(1/9.)
 Out[2]= Quaternion[1.15911 + 0.421881 I, -0.0507501 - 0.0184715 I,
  0.177625 + 0.0646503 I, -0.025375 - 0.00923576 I]
 
 In[3]:= % ** % ** % ** % ** % ** % ** % ** % ** %
 Out[3]= Quaternion[-1. + 2.24647*10^-16 I,
  2. - 8.1532*10^-16 I, -7. + 2.81719*10^-15 I, 1. - 4.0766*10^-16 I]

despite the fact, that Out[2] gives a Bi-Quaternion, non-commutative multiplication brings the original quaternion back. It is necessary to check the Quaternions package in more detail. Here is a reference

http://www.mathematica-journal.com/issue/v8i3/features/turner/contents/html/Links/index_lnk_3.html

"Real Quaternions and Rotations" by Robert Piziak and Danny W. Turner, The Mathematica Journal, vol. 8 issue 3.

Using the article's formulae one gets the ninth zeroes of Quaternion[-1, 2, -7, 1] as follows
 In[64]:= m = 9;
 q = Quaternion[-1, 2, -7, 1];
 w = PolarQtrig[q]
 w /. deg -> Degree // Simplify
 phi = angle[q]
 
 Out[66]= Sqrt[55] (Cos[
     97.7494 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[97.7494 deg])
Out[67]= (-1. + 2. I) - 7. J + 1. K
Out[68]= 97.7494 deg

In[69]:= r = Abs[q]^(1/m)
Out[69]= 55^(1/18)

In[72]:= theta = (phi/m + d 360. deg/m);
n = AdjustedSignIJK[q];
ninthrootsofq =
Table[Flatten[
   r (Cos[theta] (Quaternion[1, 0, 0, 0] // FromQuaternion) +
      Sin[theta] (n // FromQuaternion))], {d, 0, m - 1}]
numroots = ninthrootsofq /. deg -> Degree // N // Simplify

Out[74]= {55^(1/18) (Cos[0. + 10.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[0. + 10.861 deg]),
55^(1/18) (Cos[50.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[50.861 deg]),
55^(1/18) (Cos[90.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[90.861 deg]),
55^(1/18) (Cos[130.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[130.861 deg]),
55^(1/18) (Cos[170.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[170.861 deg]),
55^(1/18) (Cos[210.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[210.861 deg]),
55^(1/18) (Cos[250.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[250.861 deg]),
55^(1/18) (Cos[290.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[290.861 deg]),
55^(1/18) (Cos[330.861 deg] + (1/3 I Sqrt[2/3] - (7 J)/(3 Sqrt[6]) + K/(
       3 Sqrt[6])) Sin[330.861 deg])}

Out[75]= {(1.22698 + 0.0640715 I) - 0.22425 J + 0.0320357 K,
          (0.788599 + 0.263735 I) - 0.923072 J + 0.131867 K,
          (-0.0187746 + 0.339994 I) - 1.18998 J + 0.169997 K,
          (-0.817363 + 0.257166 I) - 0.90008 J + 0.128583 K,
          (-1.2335 + 0.0540071 I) - 0.189025 J + 0.0270036 K,
          (-1.07247 - 0.174422 I) + 0.610477 J - 0.087211 K,
          (-0.409615 - 0.321237 I) + 1.12433 J - 0.160619 K,
          (0.4449 - 0.317742 I) + 1.1121 J - 0.158871 K,
          (1.09124 - 0.165572 I) + 0.579501 J - 0.0827858 K}

(* only real coefficients, let's check the reproduction of q *)
Remove[q1, q2, q3, q4, q5, q6, q7, q8, q9]
q1 = ToQuaternion[%75[[1]]];
q2 = ToQuaternion[%75[[2]]];
q3 = ToQuaternion[%75[[3]]];
q4 = ToQuaternion[%75[[4]]];
q5 = ToQuaternion[%75[[5]]];
q6 = ToQuaternion[%75[[6]]];
q7 = ToQuaternion[%75[[7]]];
q8 = ToQuaternion[%75[[8]]];
q9 = ToQuaternion[%75[[9]]];

In[95]:= q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1 ** q1
Out[95]= Quaternion[-1., 2., -7., 1.]

In[96]:= q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2 ** q2
Out[96]= Quaternion[-1., 2., -7., 1.]

In[97]:= q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3 ** q3
Out[97]= Quaternion[-1., 2., -7., 1.]

In[98]:= q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4 ** q4
Out[98]= Quaternion[-1., 2., -7., 1.]

In[99]:= q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5 ** q5
Out[99]= Quaternion[-1., 2., -7., 1.]

In[100]:= q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6 ** q6
Out[100]= Quaternion[-1., 2., -7., 1.]

In[101]:= q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7 ** q7
Out[101]= Quaternion[-1., 2., -7., 1.]

In[102]:= q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8 ** q8
Out[102]= Quaternion[-1., 2., -7., 1.]

In[103]:= q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9 ** q9
Out[103]= Quaternion[-1., 2., -7., 1.]

Finally two solutions are possilbe
(a) you stay with the original Quaternions package but find a clever way to transform its Bi-Quaternions to Quaternions
(b) you use my proposal supplemented with the formulae of Piziak and Turner for the (1/q) part

Sincerely
Udo.
POSTED BY: Udo Krause
Answer
11 months ago