Evaluating standard functions of general quaternions symbolically

Question

Mathematica can easily compute exponentials and logarithms of concrete quaternions:

Needs["Quaternions`"]

Exp[Quaternion[1.3, 3.5, -0.7, 0.9]]
Log[Quaternion[1.3, 3.5, -0.7, 0.9]]

Quaternion[-3.14824, -1.79206, 0.358412, -0.460816]

Quaternion[1.36196, 1.17075, -0.23415, 0.301051]

Even for exact numbers:

Exp[Quaternion[1, 3, 7/2, -5]]

Quaternion[E Cos[Sqrt[185]/2], (6 E Sin[Sqrt[185]/2])/Sqrt[185],
           (7 E Sin[Sqrt[185]/2])/Sqrt[185], -2 Sqrt[5/37] E Sin[Sqrt[185]/2]]

But it seems to be unable to expand this for the general quaternion:

Exp[Quaternion[a, b, c, d]]

E^Quaternion[a, b, c, d]

I've tried using Simplify, FunctionExpand, ToQuaternion on this, but still it doesn't give a Quaternion object. But there exists a closed form for such functions.

I can, of course, define such functions myself, like:

exp[q_] = With[{v = q - Re[q]}, Exp[Re[q]] (Cos[Abs[v]] + v/Abs[v] Sin[Abs[v]]) ];
ln[q_] = Log[Abs[q]] + (q - Re[q])/AbsIJK[q] ArcCos[Re[q]/Abs[q]];

But is there a way to get Mathematica itself expand these functions of general quaternions, without me having to redefine them all?

Answer 1

The reason for the behavior you were seeing can be seen if you look through the definitions within the package. In particular, the extended elementary functions are defined with a format like func[a:Quaternion[__?ScalarQ]] := (* stuff *) where ScalarQ[] is a private package function that checks if the stuff within the Quaternion[] object are real numbers. Since the components of Quaternion[a, b, c, d] are manifestly not real (in the sense of not satisfying ScalarQ[]), the formulae within are not applied to your symbolic quaternion.

But, since this is a package function, you can modify the code yourself so that it can work with quaternions with symbolic components.

Answer 2

You can use the following code to manipulate quaternions:

Clear["Global`*"]
Unprotect[Dot];
Dot[x_?NumberQ, y_] := x y;
Protect[Dot];
Unprotect[Power];
Power[0, 0] = 1;
Protect[Power];
Matrix /: Matrix[x_?MatrixQ] := 
  First[First[x]] /; x == First[First[x]] IdentityMatrix[Length[x]];
Matrix /: NonCommutativeMultiply[Matrix[x_?MatrixQ], y_] := 
  Dot[Matrix[x], y];
Matrix /: NonCommutativeMultiply[y_, Matrix[x_?MatrixQ]] := 
  Dot[y, Matrix[x]];
Matrix /: Dot[Matrix[x_], Matrix[y_]] := Matrix[x . y];
Matrix /: Matrix[x_] + Matrix[y_] := Matrix[x + y];
Matrix /: x_?NumericQ + Matrix[y_] := 
  Matrix[x IdentityMatrix[Length[y]] + y];
Matrix /: x_?NumericQ Matrix[y_] := Matrix[x y];
Matrix /: Matrix[x_]*Matrix[y_] := Matrix[x . y] /; x . y == y . x;
Matrix /: Power[Matrix[x_?MatrixQ], y_?NumericQ] := 
  Matrix[MatrixPower[x, y]];
Matrix /: Power[Matrix[x_?MatrixQ], Matrix[y_?MatrixQ]] := 
  Exp[Matrix[y] . Log[Matrix[x]]];
Matrix /: Im[Matrix[x_?MatrixQ]] := Matrix[Im[x]]
Matrix /: Re[Matrix[x_?MatrixQ]] := Matrix[Re[x]]
Matrix /: Arg[Matrix[x_?MatrixQ]] := Matrix[Arg[x]]

$Post2 = FullSimplify[FullSimplify[# /. i -> Matrix[( {
               {I, 0},
               {0, -I}
              } )] /. j -> Matrix[( {
              {0, 1},
              {-1, 0}
             } )] /. k -> Matrix[( {
             {0, I},
             {I, 0}
            } ) ] /. 
        f_[args1___?NumericQ, Matrix[mat_], args2___?NumericQ] :> 
         Matrix[MatrixFunction[f[args1, #, args2] &, mat]]] /. 
      Matrix[{{a_, c_}, {d_, b_}}] :> 
       Re[a] + Im[a] i + Re[c] j + Im[c] k ] /. 
    Matrix[{{a_, c_}, {d_, b_}}] :> 
     Re[a] + Im[a] i + Re[c] j + Im[c] k &;
$Post = Nest[$Post2, #, 3] /. Dot -> NonCommutativeMultiply &;

After running it, you can use quaternionic $i$,$j$,$k$ inside any expressions alongside other numbers. The multiplication is evaluated only if the quaternions commute. For non-commutative multiplication one can use non-commutative multiplication operator (x ** y), it is evaluated always.

Test it:

In=Re[(j + k)^a]

Out=2^(a/2) Cos[(a Pi)/2]

In=Log[j+k]

Out=1/4 (Sqrt[2] (j+k) Pi+Log[4])

In=Sin[j + k]

Out=((j+k) Sinh[Sqrt[2]])/Sqrt[2]

In=Sqrt[j+k]

Out=(Sqrt[2]+j+k)/2^(3/4)

In=Log[j + k] ** Log[i + k]

Out=1/8 ((-1+i+j-k) Pi^2+Sqrt[2] (i+j+2 k) Pi Log[2]+2 Log[2]^2)

In=Sign[j + k]

Out=(j+k)/Sqrt[2]

In=1/(i + k)

Out=1/2 (-i-k)

In=k^j

Out=i