4 | Rotationen - Basics
4.6 Rotation um eine Achse(4) - Quaternionen
Quaternionen lassen sich als Verallgemeinerung vonC als auch vonR3 auffassen:
q :=q0+iq1+jq2+kq3 = (q0, q1, q2, q3) = q0
|{z}
Skalarteil
+ ~q
|{z}
V ektorteil
= (q0, ~q) Die Addition ist isomorph zu der inR4:
p+q= (p0+q0) + (~p+~q) = (p0+q0, ~p+~q)
Für die Multiplikation standen die kartesischen Einheitsvektoren (mit dem Kreuzprodukt im Hinterkopf) undCPate:
i2 =j2 =k2 =ijk=−1
ijk=−1
· k ⇒ ij=k ijk=−1
i· ⇒ −jk=−i
j· ⇒ k=−ji daraus folgt die Nichkommutativität
ij=k=−ji und analoge Ausdrücke
Berechnet man das Produkt zweier Quaternionen ist dies ein längerer Ausdruck, der allerdings mit Skalar- und Vektorprodukt etwas kompakter geschrieben werden kann:
pq= (p0+ip1+jp2+kp3) (q0+iq1+jq2+kq3) = (p0q0−p~·~q, p0~q+q0~p+~p×~q) (4.11) Wieder sieht man, dass durch das Vektorprodukt die Kommutativität verletzt wird. Untersucht man die beiden Rechenoperationen genauer ergibt sich als Struktur ein Nicht-Abelscher Ring.
Mit der komplexen Konjugationq∗ :=q0−~q lässt sich eine Norm |q|festlegen:
q∗q = (q0,−~q)(q0, ~q) =q02+~q·~q=q20+q21+q22+q23 =qq∗ =:|q|2 (4.12) Man kann mit obiger Multiplikation leicht zeigen:
(pq)∗ =q∗p∗ ⇒ |pq|2= (pq) (pq)∗ = (pq) (q∗p∗) =p|q|2p∗ =|p|2|q|2 Mit der Definition des Inversen ergibt sich
q−1:= q∗
|q|2 ⇒qq−1 =q−1q= 1 Für Einheitsquaternionen gilt also:
q−1 =q∗ und mit der Beziehung 1 =|q|2=q q∗ 4.12
=↓ q20+~q·~q =q02+k~qk2
4.6 Rotation um eine Achse(4) - Quaternionen
gibt es daher für jedesq mit |q|= 1einθ∈[0, π]und einen Einheitsvektornˆ ∈R3 sodass gilt:
q = (s, ~q) = (cosθ,sinθˆn) Damit können wir jetzt den Rotationsoperator festlegen:
Rq(~v) =q ~v q−1 =q ~v q∗= (s, ~q) (0, ~v) (s,−~q) (s, ~q) = (cosθ,sinθˆn)
Wir zeigen, dass Rq(~v) eine Rotation des Vektors~v um die Achsenˆ mit Betrag2θdarstellt:
Rq(~v) = (s, ~q) (0, ~v) (s,−~q) 4.11
= (s, ~↓ q) (~v~q, s~v+~q×~v) 4.11
=↓
= (s~v~q−~q·(s~v+~q×~v), s(s~v+~q×~v) + (~v~q)~q+~q×(s~v+~q×~v)) =
= (0, s2~v+ 2s~q×~v+ (~v~q)~q+~q×~q×~v) = für den letzten Ausdruck verwenden wir 4.13
= (0, ~v(s2−~q·~q) + 2(~v·~q)~q+ 2s~q×~v) mit (s, ~q) = (cosθ,sinθˆn) ergibt sich
= (0, ~v(cos2θ−sin2θ) + (~v·n)ˆˆ n2 sin2θ+ 2 cosθsinθ(ˆn×~v)
= (0, ~v(cos 2θ) + (~v·n)ˆˆ n(1−cos(2θ)) + sin 2θ(ˆn×~v) Die letzte Zeile entspricht der Rodrigues-Formel 4.8 für 2θ.
Im folgendenwxMaxima- Programm zeigen wir auch die Darstellung eines Quaternions als Ma-trix. Einige Teile gehen auf das Macro “quaternion.mac” von Gosei Furuya (war im Verzeichnis
“share/contrib/quaternion.mac” meiner Maxima-Version zu finden). Wobei einige grundlegen-de Ängrundlegen-derungen von mir gemacht wurgrundlegen-den - auf die auch hingewiesen wird. Sehen wir uns dieses
“Ungetüm” an!
4.6.1 Quaternionen in wxMaxima
Wir benutzen für das Produkt der Quaternionen das “Matrixprodukt” (Punkt) und setzen es nichtassoziativ- das stellt sicher, dass beim Arithmetik-Baum beim Punkt immer nur 2 Kindknoten vorhanden sind (wichtig!) - damit kann applyb2 den Baum “von unten” (bottom) nach den Regeln ersetzen. Würden mehr als 2 Kindknoten möglich sein, gingen die Regeln für i ·j und Co ins “Leere”. Außerdem sollen Produkte wie i · i nicht durch i2 ersetzt werden (= dotexptsimp:false )!
Abb.42 : Einfluss vondotassoc
4. Rotationen - Basics
Punktregeln festlegen
(% i1) (dotdistrib:true, dotexptsimp:false, dotscrules:true, dotassoc:false)$
i,j,k sind besondere Dinger
(% i2) declare([i,j,k],nonscalar)$
Jetzt die Produktregeln - beachte die Nichtkommutavität
(% i3) (defrule(f1,j.j,-1),defrule(f2,i.i,-1),defrule(f3,k.k,-1),defrule(f4,i.j,k),defrule(f5,j.i,-k))$
(% i4) (defrule(f6,j.k,i), defrule(f7,k.j,-i), defrule(f8,k.i,j), defrule(f9,i.k,-j))$
Das Arbeitspferd für die Multiplikation - rekursiv werden die Regeln - von den “Blättern” des Baums ausgehend(applyb2 - “b” wie bottom) - angewandt und anschl.i,j,k herausgehoben (% i6) expand4(_ex):=block([_eex],_eex:applyb2(expand(_ex),f1,f2,f3,f4,f5,f6,f7,f8,f9),
collectterms(_eex,i,j,k))$
Jetzt folgen einige “utility”-Routinen wie z.B.: is Quaternion?
(% i7) quaternionp(q):= if (abs(coeff(q,i)) + abs(coeff(q,j)) + abs(coeff(q,k)) # 0) then true else false$
q → q∗
(% i8) conj4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
if (listp(_ex1)) then [_ex1[1], -_ex1[2], -_ex1[3], -_ex1[4]]
else scalarpart4(_eex1) - vectorpart4(_eex1) )$
Wir legen die Norm fest - gilt auch bei Übergabe als Liste!
(% i9) norm4(_ex1):=if (listp(_ex1)) then sqrt(lreduce("+",_ex1 * _ex1)) else sqrt(expand4(_ex1 . conj4(_ex1)))$
q → q−1
(% i10) inv4(_ex):= block( if (norm4(_ex)6=0) then conj4(_ex)/norm4(_ex)ˆ2 else false)$
Skalaranteil zurückgeben
(% i11) scalarpart4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
_ex1:coeff(_eex1,i)*i+coeff(_eex1,j)*j+coeff(_eex1,k)*k,(_eex1-_ex1))$
Vektoranteil zurückgeben
(% i12) vectorpart4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
(_eex1,i)*i+coeff(_eex1,j)*j+coeff(_eex1,k)*k)$
4.6 Rotation um eine Achse(4) - Quaternionen
Konvertierung in Zeilenvektor (Liste), Vorsicht: keine Matrix!
(% i13) convert2vec(_ex1):=block([_eex1],_eex1:expand4(_ex1),
[scalarpart4(_eex1),coeff(_eex1,i),coeff(_eex1,j),coeff(_eex1,k)])$
Definition von Skalaren: inpart(a·b,0)→ ∗, dot wird “normale” Multiplikation (% i14) declare([a,b,c,d,e,f,g,h],scalar)$
4 Dezimalen-Ausgabe, keine Konvertierungswarnungen, brauchen Levi-Civita Symbolεijk (% i15) (fpprintprec:4,ratprint:false, load(itensor))$
i,j und k werden fett gedruckt beitex - Ausgabe
(% i16) (texput(i,"\\textbf{i}"),texput(j,"\\textbf{j}"),texput(k,"\\textbf{k}"))$
Wir konstruieren 2 “Probe”-Quaternionen und berechnen das Produkt (% i17) (q1:a+b*i+c*j+d*k, q2:e+f*i+g*j+h*k)$
(% i18) p_q1q2:expand4(q1 . q2);
(a h+b g−c f +d e) k+(−b h+a g+d f +c e)j+(c h−d g+a f +b e)i−d h−c g−b f+a e Ausgabe wurde oben verwendet
(% i19) tex(p_q1q2,false);
Nur für Testzwecke
(% i20) vectorpart4(p_q1q2)+scalarpart4(p_q1q2);
(a h+b g−c f +d e) k+(−b h+a g+d f +c e)j+(c h−d g+a f +b e)i−d h−c g−b f+a e Ausgabe siehe oben
(% i21) tex(vectorpart4(p_q1q2)+scalarpart4(p_q1q2),false);
Konvertierung in einen Zeilen- bzw. Spaltenvektor (% i22) v1:convert2vec( p_q1q2 );
[−dh−cg−bf+ae, ch−dg+af+be,−bh+ag+df+ce, ah+bg−cf +de]
(% i23) v2:transpose(v1);
−dh−cg−bf+ae ch−dg+af+be
−bh+ag+df +ce ah+bg−cf+de
Diesen Vektor kann man offensicht-lich als Linearkombination mit den Linearfaktoren h, g, f, e darstellen - und somit nach 4.1.2 als Pro-dukt einer Matrix mit dem Vektor (h, g, f, e)t
4. Rotationen - Basics
(% i24) q2AsVec:convert2vec(q2); [e, f, g, h]
Die Zerlegung von(% i23) gelingt nur auf “Umwegen”:
inpart(v1[2],3)→ holt vom 2.-ten Listeneintrag den 3.-ten Summanden also → −dg partition(inpart(v1[2],3),d)→ erstes Listenelement ist Faktor von dalso → −g Wir holen die Ausdrücke woe, f, g und h vorkommen als Vektoren:
partition(inpart(v1[1],2),q2AsVec[2]) → [−b, f] Wir verallgemeinern:
partition(inpart(v1[i],j),q2AsVec[j]) →mit i∈ {1. . .4}, j∈ {1. . .4}
nehmen nur das erste Element der Liste (den Faktor) und stellen damit eine Liste zusammen (makelist), bauen einen Spaltenvektor(transpose) und stellen eine Liste der Spaltenvektoren zusammen (makelist) - nicht gerade elegant, aber funktioniert.
(% i25) vecList:makelist(transpose(makelist(first(partition(inpart(v1[i],j),q2AsVec[j])),i,1,4)),j,1,4);
Wir schalten den “Simplifier” ab, sonst multipliziert er wieder aus!
(% i26) simp:false$
Damit können wir das Produkt q1·q2 wie unten schreiben - nach 4.1.2 ist das gleich dem Produkt einer Matrix mit dem Spaltenvektor [e,f,g,h]
(% i27) ex2:vecList[1]*e+vecList[2]*f+vecList[3]*g+vecList[4]*h;
(% i29) is (equal(ex2 ,v2)); true Bestätigung der Gleichheit Wir bauen aus den Spaltenvektoren eine Matrix
(% i30) buildQMatrix():=block([A:matrix([]) ],for i: 1 thru 4 do A: addcol(A,vecList[i]),A)$
(% i31) B:buildQMatrix();
(% i32) is (equal(B . q2AsVec, v2)); true Passt ja!
4.6 Rotation um eine Achse(4) - Quaternionen
Konvertierung für später
(% i35) Vec2Complex(v):=v[1]+i*v[2]+j*v[3]+k*v[4]$
Wie ist diese MatrixMaufgebaut: 1) antisymmetrisch:M =L−Lt List linke Dreiecksmatrix 2) Hauptdiagonale ist der Skalarteil:M+D (wobei DDiagonalmatrix)
3) 1. Spalte die Vektordarstellung:L[i,1]←q[i], damit sind bis auf 3 alle Elemente bestimmt:
Genauere Untersuchung ergibt: B[i, j] :=
4
X
k=2
−εkij∗q[k] i∈ {3,4}, j ∈ {2. . .(i−1)}
dabei ergeben sich die richtigen Vorzeichen. (Man könnte auch die 3 Werte einfach in die Matrix schreiben, aber . . . ). Damit haben wir die Minor-Matrix erledigt.
(% i36) qAsMat(q):=block([q_row, m:zeromatrix(4,4), d ], /* convert input in row-vector (list)!
and leave inner block, when finished */
block(
if (listp(q)) then (q_row:q, return(0)),
if (scalarp(q)) then (q_row:cons(q,[0,0,0]), return(0)), if (matrixp(q)) then (q_row:flatten(args(q)), return(0)), if (quaternionp(q)) then (q_row:convert2vec(q))
),
/* end of inner block - q_row is a list now! */
d:diagmatrix(4,q_row[1]),
/* first column of matrix is set to q */
for i:2 thru 4 do m[i,1]: q_row[i], /* Now formula above is applied */
for i:3 thru 4 do
for j:2 thru (i-1) do m[i,j]: sum(-levi_civita([k,i,j])*q_row[k],k,2,4), /* skew matrix is built with q[1] as diagonal-elements */
m-transpose(m)+d )$
Jetzt der Test:
(% i37) q:[a,b,c,d]$
4. Rotationen - Basics Jetzt die Matrix des inversen Quaternions
(% i39) q1_inv_Mat:qAsMat(1/norm4(q)ˆ2*conj4(q))$
(% i40) q1_inv_Mat . q1Mat;
Jetzt lösen wir unsere Testaufgabe(Rotation) mit Quaternionen - zuerst die Angaben (% i41) (axis:[2,-2,1], M:[0.3,0.2,0.2], X:[1,0.5,0.5] , rotAngle:float(%pi/3) )$
Einheitsvektornˆ in Richtung Drehachse
(% i42) n_0:float(1/norm4(cons(0,axis)) * axis)$
Wir erstellen das Rotations-Quaternion - Achtung halber Drehwinkel!
(% i44) q_rot_vec:cons(cos(rotAngle/2), sin(rotAngle/2)*n_0)$
(% i45) q_rot:Vec2Complex(q_rot_vec)$
Der rotierte Punkt als Quaternion - Realteil (Rundungsfehler) verschwindet!
(% i51) expand4(q_rot . Vec2Complex(cons(0,X-M)) . conj4(q_rot)) + Vec2Complex(cons(0,M));
0.9885∗k+ 0.2566∗j+ 0.5124∗i+ 3.469∗10−18 Jetzt mit Quaternion-Matrizen - transpose nur aus Platzgründen
(% i53) transpose(qAsMat( qAsMat(q_rot) . cons(0,X-M) ) . conj4(q_rot_vec) + cons(0,M));
0.0 0.5124 0.2566 0.9885
4.6.2 Quaternionen in GNU-Octave
In GNU-Octave braucht man das Rad nicht neu erfinden - es existiert ein Paket, welches die wesentlichen Operationen zur Verfügung stellt - insbesondere die Funktion rot2q ist bereits vorhanden - Eingabe: Einheitsvektor der Rotationsachse und Rotationswinkel.
Hier das Listing von SkriptrotByQuat.m:
4.6 Rotation um eine Achse(4) - Quaternionen
1 ## make sure that the package i s loaded pkg load quaternion ;
3
## convert quaternion to column−v e c t o r
19 f u n c t i o n v=q2vec ( q )
41 X_r=q2vec ( q_r∗ vec2q ( point−axis_point ) ∗conj( q_r ) )+axis_point endfunction;
43## output
X_r=rotPoint (X, rotAngle ,axis,M)
Output:
Auf weitergehende Betrachtungen sei hier verzichtet (Differentiation, Integration, Exponential-und Logarithmusfunktion, algebraische Eigenschaften, usw.) .
Fazit: Sind Translationen im Spiel haben die Matrizen leichte Vorteile, da sie ebenfalls über homogene Koordinaten als Matrizen darstellbar sind, ansonsten sind Quaternionen wegen ihrer reichhaltigen algebraischen Struktur beinahe unschlagbar.
4. Rotationen - Basics
4.7 ANHANG
Definition 4.3 Levi-Civita Symbol
Seien~e1,~e2 und~e3 die Standardbasisvektoren des 3 dim Vektorraums, dann defi-nieren wir
εijk ist die Determinante der Matrix, die aus obigen Zeilenvektoren besteht!
Nach den Eigenschaften der Determinante giltεijk= (−1)·εjik. 2 Zeilenvertauschungen führen wieder zur Identität, alsoεjk`=ε`kj
Dies gilt natürlich für alle(!) Vertauschungen.
Wir benutzen jetzt folgende Eigenschaften der Determinante:
det(A·B) =det(A)·det(B) und det(AT) =det(A)
Theorem 4.4 Multiplikation von Levi-Civita Symbolen
εijkε`mn =
Theorem 4.5 Korollar aus 4.4
Mit dem Laplace’schen Entwicklungssatz(Entwicklung nach erster Zeile) undi=` folgt mit Theorem 4.4
εijkεimn:=
3
X
i=1
εijkεimn =δjmδkn−δkmδjn
4.7 ANHANG
Theorem 4.6 Graßmann-Identität
~a×(~b×~c)
| {z } d~
= (~a·~c)~b−(~a·~b)~c (4.13)
Beweis: Diej-te Komponente des Kreuzproduktes (mit Summenkonvention) ist (~a×d)~j =εjk`akd`
Wir berechnen diej-te Komponente von 4.13:
εjk`
=|{z}ε`jk
ak ε`mnbmcn
| {z } d`
= (δjmδkn−δjnδkm)akbmcn=akbjck−akbkcj
wobei wir Theorem 4.5 benutzt haben.
Theorem 4.7 Orthogonalität des Kreuzproduktes
(~a×~b)⊥~a ∧ (~a×~b)⊥~b
Beweis:
(~a×~b)·~a=εik`akb`ai
Da gilt εik` =−εki` kommen in obiger Summe die Summandenakb`aisowohl mit positiven als auch negativen Vorzeichen vor - die Summe verschwindet also!
Die analoge Situation ergibt sich für~b.
4. Rotationen - Basics
Theorem 4.8 Wiederholtes Kreuzprodukt mit einem Einheitsvektor Sei ˆnein Einheitsvektor des R3, dann gilt
[ˆn]3×=−[ˆn]× ⇔ ~v×××=−~v× ⇔ nˆ×(ˆn×(ˆn×~v)) =−(ˆn×~v) aus der ersten Identität folgt unmittelbar
[ˆn]k+2× =−[ˆn]k× k∈ {1,2,3. . .}
Beweis: Wir setzen ~a:= ˆn×~v wobei wir aus Theorem 4.7 wissen:~a⊥nˆ Wir müssen also zeigen: nˆ×(ˆn×~a) =−~a
Die i-te Komponente von nˆ×~a ist:
(ˆn×~a)i =εijknjak =:xi
Die α-Komponente von nˆ×~x ist:
(ˆn×~x)α =εαβinβxi =εαβinβεijknjak
Beachte: Der einzige unabhängige Index ist α - über alle anderen wird summiert. Es handelt sich also um einen Vektor!
Mit Theorem 4.5 lösen wir das Produkt der Levi-Civita Symbole auf zu ˆ
n×(ˆn×~a)α = [δαjδβk−δαkδβj]nβnjak =nkak
| {z } 0
nα−njnj
| {z } 1
aα =−aα
Der erste Term verschwindet wegen der Orthogonalität von nˆ und~aund njnj = 1, danˆ ein Einheitsvektor ist!