# Writing intelligent replacement rules for symbols

I am learning about quaternions right now, and I wanted to play around with some concepts I am learning.

e.g. take rotation of a point about an axis with specified angle.

th = 3.4;
v  = Sqrt[1/3.0]*(i + j + k)
q    = Cos[th] + Sin[th]*v;
qinv =  Cos[th] - Sin[th]*v;

p = 2*i + 3*j + 4*k;

r = q * p *qinv // Expand
r = r//.{ i*i->-1, j*j->-1,k*k->-1,
i*j->k , j*i->-k,
j*k->i,  k*j->-i,
k*i->j,  i*k->-j}


Here I want to rotate the point p in three dimensions by the angle th with the axis given by the unit vector v.

So for that I basically just need to use replacement rules like did above.

The result I get is

0.391808 + 2.26121 i - 0.0435342 i^3 + 3.1959 j - 0.0653013 j^3 +
4.1306 k - 0.0870683 k^3


How do I tell Mathematica to replace 'i^3' etc. by '-i'? I would have thought Mathematica would be smart enough to use my 'i^2=-1' rule to figure this out.

Generalizing this question, if I am manipulating symbols (not just quaternions) how do I specify replacements so that the final answer is as 'simple' as possible?

-
Have a look at this question –  Jens Sep 4 '13 at 6:25
also in Mathematica i*j and j*i are the same. So when you write a rule for i*j and a different rule for j*i, then Mathematica will use only one, depending on the full form. Look at FullForm[i*j] and FullForm[j*i], they are the same. But for i^3 to -i part, try this rule: p1 = Power[m_, n_] :> If[n > 2, -m^(n - 2), -m]; then do i^3//.p1 –  Nasser Sep 4 '13 at 6:29
@smilingbuddha you might want to check out the Quaternions package in Mathematica. –  chuy Sep 4 '13 at 14:40

Replacement rules, at least in a direct form, are frequently not the right choice for mathematical manipulation because they are structural tools.

In this case /. {i^3 -> -i} works as expected, because that matches the structure of the expression. You will have to write a rule that matches the structure and then performs your transformation if you want it to work with various powers of i. For example:

expr /. i^x_ /; x >= 2 :> (-1) i^(x - 2)


0.391808 + 2.30474 i + 3.1959 j - 0.0653013 j^3 + 4.1306 k - 0.0870683 k^3

-
This is part of the answer. However, the rules also require non-commuting operations, so even the starting point using Power or Times here isn't actually correct for the desired application. –  Jens Sep 4 '13 at 6:26
@Jens Okay, I await your posting the correct answer. :-) –  Mr.Wizard Sep 4 '13 at 6:28

For this problem, I wouldn't use a rule-based approach to realize the algebraic relations. Instead, if you don't want to load the Quaternions package, the most robust way to implement quaternions is using their matrix basis representation, see e.g. this Mathworld post.

Essentially, the basis vectors can be represented as multiples of the Pauli matrices, so I just modified my recent answer on the Pauli algebra to make it work in this case:

For simplicity, I use the symbols K[0] for 1, K[1] for i, K[2] for j and K[3] instead of k because the programming is easier when using indices for the four components. Here is a function that takes a quaternion expression in the component notation and simplifies it using the matrix representation. The multiplication between quaternion components is always required to be written as a Dot product because the underlying non-commutative algebra is realized with matrices:

Clear[quaternionReduce]
quaternionReduce[a_] := Module[{
x,
symbolicQIndices,
expression,
component = ({1, I, I, I} PauliMatrix[{0, 3, 2, 1}])
},
x = Array[\[FormalX], 4];
symbolicQIndices =
DeleteDuplicates[Cases[a, K[i_Symbol] :> i, Infinity]];
expression =
Array[K, 4,
0].(x /.
First[Solve[
x.component == a /.
K[i_] :>
Sum[KroneckerDelta[i, k] component[[k + 1]], {k, 0, 3}],
x]]);
FullSimplify[expression,
Assumptions ->
Map[# \[Element] Integers && 1 <= # <= 3 &, symbolicQIndices]]]


Here are some examples for doing algebra with this function:

quaternionReduce[K[1].K[2]]

(* ==> K[3] *)

quaternionReduce[K[2].K[1]]

(* ==> -K[3] *)

quaternionReduce[K[1].K[2].K[3]]

(* ==> -K[0] *)

quaternionReduce[MatrixPower[K[1], 3]]

(* ==> -K[1] *)

quaternionReduce[MatrixPower[(K[0] + K[1]), 2]]

(* ==> 2 K[1] *)


The K[0] above is the real unit 1, as mentioned above. In principle you could make the input and output fancier (i.e., more traditional symbols), but I wanted to keep the programming effort small.

The last two relations address the problem that led to the hang-up in your question. Here, powers are realized using MatrixPower and everything works as expected.

-