Re quaternions, second plea for help
- To: MATHGROUP at yoda.physics.unc.edu
- Subject: Re quaternions, second plea for help
- From: PMG0393 at v2.qub.ac.uk
- Date: Tue, 13 Oct 92 9:16 GMT
I wrote this little package as a demonstration for a course I teach and it is not intended to be complete in any way, but it may give you some ideas. Notes: (1) It was written under Mathematica 1.2 but seems to work fine under Mathematica 2.0. (2) Error handling needs massive improvement---try entering i**j--j**i and you will see that the resulting message is far from informative. (3) Output formatting could be improved by, for example, replacing expressions like `1 j' by `j' but this would require far too many rules the way that I have done it. (4) There must be a far better way of doing this, perhaps such an inelegant solution as this will prompt someone to come up with something better. A.W. Wickstead Pure Mathematics Dept. The Queen's University of Belfast PMG0393 at VAX2.QUEENS-BELFAST.AC.UK (*------------------------------Cut Here----------------------------------*) BeginPackage["Quaternions`"] QuatAbs::usage="QuatAbs[a+b i+c j+d k] gives Sqrt[a2+b2+c2+d2]" QuatConjugate::usage="QuatConjugate[a+b i+c j+d k] gives a-b i-c j-d k" Quaternion::usage="Quaternion[a,b,c,d] is used internally to represent a quaternion a+b i+c j+d k. Quaternion multiplication is denoted by **." QuatInverse::usage="QuatInverse[quaternion] gives the inverse of quaternion." Unprotect[Plus,Times,Power,NonCommutativeMultiply,Divide,i,j,k,Equal,Unequal, QuatAbs,QuatConjugate,Quaternion,QuatInverse] Clear[Plus,Times,Power,NonCommutativeMultiply,Divide,i,j,k,Equal,Unequal, QuatAbs,QuatConjugate,Quaternion,QuatInverse] (*Begin["`Private`"]*) (*All quaternions will be kept in a structure Quaternion[a,b,c,d], which will correspond to a+b i +c j+d k. We need to know when something is a scalar, rather than a quaternion. The function Head[] gives the outermost function in an expression.*) Scalar[x_]:=Not[TrueQ[Head[x]==Quaternion]] (*in order to get some Quaternion[] objects to work with, we replace i, j and k by the appropriate object. If we tried to replace the number 1 by Quaternion[1,0,0,0] then we would have to replace the 1 here and so on*) i:=Quaternion[0,1,0,0] j:=Quaternion[0,0,1,0] k:=Quaternion[0,0,0,1] (*Quaternion addition, addition of scalars to quaternions and scalar multiplication of quaternions deals with the rest of input formatting*) Quaternion[a_,b_,c_,d_]+Quaternion[e_,f_,g_,h_]:=Simplify[Quaternion[a+e,b+f,c+g ,d+h]] x_ + Quaternion[a_,b_,c_,d_]:=Simplify[Quaternion[x+a,b,c,d]]/;Scalar[x] x_ Quaternion[a_,b_,c_,d_]:=Simplify[Quaternion[x a,x b,x c,x d]]/;Scalar[x] (*Here is the rest of quaternion arithmetic*) (* The quaternion multiplication is denoted by **, use of * between quaternions is a mistake to be flagged*) Quaternion[a_,b_,c_,d_]**Quaternion[e_,f_,g_,h_]:= Simplify[Quaternion[a e-b f-c g-d h,a f+b e+c h-d g, a g-b h+c e+d f,a h+b g-c f+d e]] Quaternion[a_,b_,c_,d_]*Quaternion[x_,y_,z_,t_]:="Error in quaternion expression" QuatConjugate[Quaternion[a_,b_,c_,d_]]:=Quaternion[a,-b,-c,-d] QuatModsquare[Quaternion[a_,b_,c_,d_]]:=Simplify[a*a+b*b+c*c+d*d] QuatAbs[Quaternion[a_,b_,c_,d_]]:=Sqrt[Simplify[a*a+b*b+c*c+d*d]] QuatInverse[Quaternion[x_,y_,z_,t_]]:=Simplify[QuatConjugate[Quaternion[x,y,z,t] ]/QuatModsquare[Quaternion[x,y,z,t]]] (*Quaternion powers are defined inductively in the obvious manner*) Quaternion[x_,y_,z_,t_]0:=Quaternion[1,0,0,0] Quaternion[x_,y_,z_,t_]1:=Quaternion[x,y,z,t] Quaternion[x_,y_,z_,t_]n_Integer?Positive:=Simplify[Quaternion[x,y,z,t]**(Quate rnion[x,y,z,t](n-1))] Quaternion[x_,y_,z_,t_]n_Integer?Negative:=Simplify[QuatInverse[Quaternion[x,y, z,t](-n)]] (*In general quaternion division is ambiguous, so we leave it undefined. The only exception is real/quaternion, which is OK*) a_/Quaternion[x_,y_,z_,t_]:=Simplify[a QuatInverse[Quaternion[x,y,z,t]]] (*Test for equality and inequality must compare real with quaternions as well as two quaternions*) Equal[Quaternion[a_,b_,c_,d_],Quaternion[x_,y_,z_,t_]]:=a==x&&b==y&&c==z&&d==t Unequal[Quaternion[a_,b_,c_,d_],Quaternion[x_,y_,z_,t_]]:=!(a==x&&b==y&&c==z&&d= =t) Equal[Quaternion[a_,b_,c_,d_],x_]:=a==x&&b==0&&c==0&&d==0/;Scalar[x] Unequal[Quaternion[a_,b_,c_,d_],x_]:=!(a==x&&b==0&&c==0&&d==0)/;Scalar[x] Equal[x_,Quaternion[a_,b_,c_,d_]]:=a==x&&b==0&&c==0&&d==0/;Scalar[x] Unequal[x_,Quaternion[a_,b_,c_,d_]]:=!(a==x&&b==0&&c==0&&d==0)/;Scalar[x] (*It is printing answers that is tedious to enter. The many cases are necessary to prevent results like 1+0 i+0 j+2 k, but results would be intelligible even if we were not too careful about this detail. The HoldForm[] prints its argument without evaluating it, otherwise we would take Quaternion[0,0,c,0] and try to print it as "c j" which is in turn evaluated to Quaternion[0,0,c,0] etc. giving an infinite loop*) Format[Quaternion[0,0,0,0]]:=HoldForm[0] Format[Quaternion[0,0,0,d_]]:=HoldForm[d k] Format[Quaternion[0,0,c_,0]]:=HoldForm[c j] Format[Quaternion[0,b_,0,0]]:=HoldForm[b i] Format[Quaternion[a_,0,0,0]]:=HoldForm[a] Format[Quaternion[0,0,c_,d_]]:=HoldForm[c j+d k] Format[Quaternion[a_,0,0,d_]]:=HoldForm[a+d k] Format[Quaternion[a_,b_,0,0]]:=HoldForm[a+b i] Format[Quaternion[0,b_,c_,0]]:=HoldForm[b i +c j] Format[Quaternion[a_,0,c_,0]]:=HoldForm[a+c j] Format[Quaternion[0,b_,0,d_]]:=HoldForm[b i +d k] Format[Quaternion[a_,b_,c_,0]]:=HoldForm[a+b i +c j] Format[Quaternion[a_,b_,0,d_]]:=HoldForm[a+b i +d k] Format[Quaternion[a_,0,c_,d_]]:=HoldForm[a +c j+d k] Format[Quaternion[0,b_,c_,d_]]:=HoldForm[b i +c j+d k] Format[Quaternion[a_,b_,c_,d_]]:=HoldForm[a+b i +c j+d k] (*End[]*) Protect[Plus,Times,Power,NonCommutativeMultiply,Divide,i,j,k,Equal,Unequal, QuatAbs,QuatConjugate,Quaternion,QuatInverse] EndPackage[]