DE102023100648B3 - Verfahren und Messeinrichtung zur Bestimmung der Lage eines Objekts - Google Patents

Abstract

Es wird ein Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehwinkel des Objektes beschrieben, wobei relative Drehwinkelinkremente Δθkan aufeinanderfolgenden Messzeitpunkten tkmit k = 1, ..., n gemessen und für die jeweiligen Messzeitpunkte tkzu kumulierten Drehwinkelinkrementen θkzusammengefasst werden. Es erfolgt ein Berechnen der kumulierten Drehwinkelinkremente θkfür die Messzeitpunkte tkiterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθk, die an aufeinanderfolgenden Messzeitpunkten tkmit k = 1, ..., n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θk-1und der transzendenten FunktionB(ζ)=1ζ2(1−ζ cot ζ) mit ζ als Funktion des kumulierten Drehwinkelinkrementes θk-1zum jeweils vorhergehenden Messzeitpunkt tk-1, wobei ein Approximieren der transzendenten FunktionB(ζ)=1ζ2(1−ζ cot ζ)durch Padé-Approximationen. Durch diese Padé-Approximationen der transzendenten Funktion kann auch sehr ressourceneffizient und genau eine automatische rechnergestützte Berechnung von trigonometrischen Funktionen von Winkeln θkoder Drehungen aus gegebenen Winkeln θkvorgenommen werden.

Description

• Das Projekt, das zu dieser Anmeldung geführt hat, wurde durch das Forschungs- und Innovationsprogramm Horizont 2020 der Europäischen Union unter der Fördervereinbarung Nr. 101004205 gefördert.
• Die Erfindung betrifft ein Verfahren zum Bestimmen der Lage eines Objektes aus den Drehraten des Objektes, wobei relative Drehwinkelinkremente Δθ^k an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen und für die jeweiligen Messzeitpunkte t_k zu kumulierten Drehwinkelinkrementen θ^k zusammengefasst werden.
• Die Erfindung betrifft weiterhin eine Messeinrichtung zum Bestimmen der Lage eines Objektes aus den Drehraten des Objektes, wobei die Messeinrichtung eine Sensoreinheit zum Messen von relativen Drehwinkelinkrementen Δθ^k für aufeinanderfolgende Messzeitpunkte t_k mit k = 1, ..., n und eine Datenverarbeitungseinheit hat, die eingerichtet ist, die an aufeinanderfolgenden Messzeitpunkten t_k gemessenen relativen Drehwinkelinkremente Δθ^k für die jeweiligen Messzeitpunkte t_k zu kumulierten Drehwinkelinkrementen θ^k zusammenzufassen.
• Gyroskope dienen der Messung der Drehrate (Winkelgeschwindigkeit) um eine vordefinierte Achse. Bei Verwendung mindestens dreier Gyroskope, die so orientiert sind, dass die Richtungsvektoren entlang der Achsen ihrer Messung linear unabhängig sind, lässt sich der Winkelgeschwindigkeitsvektor ω einer Drehung um eine beliebig im Raum orientierte Drehachse bestimmen.
• Wenn die Gyroskope fest mit einem Objekt, wie insb. einem Fahrzeug (z.B. einer Rakete) verbunden sind, lässt sich aus der Drehrate durch Integration des Signals die (relative) Lageänderung des Fahrzeugs (Orientierung im Raum) berechnen. Bei bekannter Anfangsorientierung ist dann die absolute Orientierung des Fahrzeugs im Raum bekannt. Die Kenntnis der Orientierung ist die Eingangsgröße für die Lageregelung und für die Bestimmung der Position, wenn orientierungsabhängige Kräfte auf das Fahrzeug wirken, z.B. verursacht durch fahrzeugfeste Antriebs- oder Bremssysteme.
• Üblicherweise werden die Gyroskope mit Beschleunigungssensoren zur Messung der Linearbeschleunigung in einem Sensorpaket (einer Inertialmesseinheit bzw. englisch: Inertial Measurement Unit (IMU)) gebündelt, welches als Modul verfügbar ist und zentraler Baustein für die Bestimmung der Orientierung und Position des Fahrzeugs ist. Die IMU enthält neben den Sensoren auch noch Stufen der Signalverarbeitung, mit dem Ziel an der Benutzerschnittstelle ein im Boardcomputer weiterzuverwendendes Signal zu erzeugen. Nach der Digitalisierung und evtl. Filterung der Signale erfolgt eine Reduktion der Frequenz, so dass an der Schnittstelle ein die Spezifikationen des Schnittstellenprotokolls erfüllendes hinreichend niederfrequentes Signal zur Weiterverarbeitung zur Verfügung steht.
• Im Fall der Gyroskope werden oft nicht direkt die Drehraten ω gemessen, sondern ihre Integrale, d.h. die über einen Abtastzeitschritt Δt_i integrierte (d.h. kumulierte) Drehrate, die z.B. im Zeitschritt von t_i-1 nach t_i, i ∈ N durch $$\Delta {\phi }^i={\displaystyle {\int }_{t_{i-1}}^{t_{i-1}+\Delta t_i}\omega \left(t\text{'}\right)\text{d}t\text{'}}$$

  

  gegeben ist. Im Fall einer zeitlich konstanten Drehachse entspricht dieses Integral dem (relativen) Drehwinkel (Winkelinkrement) Δθⁱ, um den sich die IMU während dieses Abtastzeitschrittes Δt_i gedreht hat. Soll nun an der Benutzerschnittstelle ein kumulierter Drehwinkel $\Delta {\theta }_{out}^I$

  

  ausgegeben werden, der sich auf einen längeren Zeitbereich $\Delta t_{I,out}={\displaystyle {\sum }_{k=1}^n\Delta t_k,n\in \text{N}}$

  

  von t_I-1 = t_i-1 nach $t_I=t_{i-1}+{\displaystyle {\sum }_{k=1}^n\Delta t_k}$

  

  (also eine gegenüber der Eingangsfrequenz f niedrigere Frequenz f_out < f) bezieht, so sind die n individuellen Drehwinkel Δθ^i-1+k, k = 1,..., n, die innerhalb des Zeitintervalls [t_I-1, t_I] gemessen werden, aufzuaddieren.
• Zum Verständnis der Notation wird darauf hingewiesen, dass das Zeitintervall [t_I-1, t_I], das im Folgenden auch als Zeitbereich I bezeichnet wird, durch die Folge von n Messzeitpunkten t_meas,k, k = 1,..,n, wobei t_n = t_I ist, in n Zeitintervalle [t_meas,k-1, t_meas,k] aufgeteilt wird, wobei die individuellen Längen Δt_k = t_meas,k-t_meas,k-1 unterschiedlich sein können. Der I-te Zeitbereich kann sich an einen vorausgehenden, (I-1)-ten Zeitbereich [t_I-2, t_I-1]) anschließen.
• Beträgt beispielsweise die Abtastfrequenz der Gyroskope f = 2 kHz, so dass die gemessenen Winkelinkremente auf eine konstante Abtastrate von Δt = f^-1 = 5 * 10^-4 s bezogen sind und soll die Frequenz an der Benutzerschnittstelle f_out = 100 Hz betragen, so sind jeweils $n=\frac{f}{f_{out}}=20$

  

  der mit Zeitabständen Δt = 5 * 10^-4 s gemessenen Drehwinkelinkremente Δθ^k zu summieren, um das Winkelinkrement bezogen auf einen Zeitschritt Δt_out = 1 / f_out = 0,01 s am Ausgang zu erhalten, also den relativen Drehwinkel, um den sich die IMU während eines Zeitschritts Δt_out gedreht hat.
• Diese Überlegungen sind allerdings nur korrekt, wenn die Rotation um eine zeitlich konstante Drehachse erfolgt. Diese Annahme ist jedoch in der Regel nicht erfüllt. Rotiert beispielsweise ein starrer Körper ohne Einwirkung äußerer Drehmomente um eine Achse, die keine der drei Hauptträgheitsachsen ist, so führt der Winkelgeschwindigkeitsvektor im körperfesten Bezugssystem eine Präzessionsbewegung aus und seine Richtung ist nicht zeitlich konstant. Im Kontext der Navigation wird dieses Phänomen als „Coning“ bezeichnet. Im Fall einer zeitlich veränderlichen Drehachse ergeben sich nun folgende Änderungen:
• 1 . Das von den Sensoren gemessene Integral Δφ^k der Drehrate ω ergibt nicht direkt das Winkelinkrement Δθ^k, das direkt Drehachse und den Drehwinkel der im Abtastintervall erfolgten Drehung angibt. Stattdessen ist das relative Drehwinkelinkrement Δθ^k aus den Integralen der aktuellen und vergangenen Drehraten zu berechnen.
• 2. Die direkte Summation zeitlich aufeinanderfolgender mit den Zeitintervallen [t_i-2+k, t_i-1+k], [t_i-1+k, t_i+k] assoziierter Winkelinkremente Δθ^k, Δθ^k+1 ergibt nicht das zur Drehung im Zeitintervall [t_i-2+k, t_i+k] gehörende Winkelinkrement. Anstelle der einfachen Summation der Winkelinkremente ist eine kompliziertere mathematische Operation auszuführen.
• John E. Bortz. „A New Mathematical Formulation for Strapdown Inertial Navigation." In: IEEE Transactions on Aerospace and Electronic Systems AES-7.1 (1971), pp. 61-66. DOI: 10.1109/TAES.1971.310252, schlägt eine (nichtlineare) Differentialgleichung für den Drehwinkel θ(t) als Funktion der Zeit vor, die von der Drehrate ω(t) abhängt und wie folgt lautet: $$\begin{array}{l}\dot{\theta }=\frac{\mathrm{\theta }}{2}\text{cot}\frac{\mathrm{\theta }}{2}\mathrm{\omega }+\left(1-\frac{\mathrm{\theta }}{2}cot\frac{\mathrm{\theta }}{2}\right)\frac{\mathrm{\theta }\cdot \mathrm{\omega }}{{\mathrm{\theta }}^2}\mathrm{\theta }+\frac{1}{2}\mathrm{\theta }\times \mathrm{\omega },\\ \text{mit }\mathrm{\theta }\cdot \dot{\theta }=\mathrm{\theta }\cdot \mathrm{\omega =\theta }\dot{\theta },\mathrm{\theta =}{\Vert \theta \Vert }_2.\end{array}$$

  

  In der obigen Gleichung ist der mehrkomponentige Drehwinkel θ(t) in Fettdruck dargestellt, um ihn von seiner Norm θ = ||θ||₂ zu unterscheiden.
• Diese Differentialgleichung kann für kleine Drehwinkel θ vereinfacht werden, indem sie um Drehwinkel θ = 0 in eine Taylorreihe entwickelt wird, wobei nur die führende(n) Ordnung(en) beibehalten werden. Sie lässt sich dann beispielsweise in einem Zeitintervall [t_I-1, t_I], also für das Winkelinkrement $\Delta {\theta }_{out}^I$

  

  lösen. Hierzu wird ω(t) in diesem Intervall durch ein Polynom (n - 1)-ten Grades approximiert durch $$\omega \left(t\right)={\displaystyle {\sum }_{k=1}^n{\omega }_0^{\left(k-1\right)}\frac{{\left(t-t_{i-1}\right)}^{k-1}}{\left(k-1\right)!}}.$$

  
• Dessen n Koeffizienten können aus n gemessenen integrierten Drehraten Δφ^k, k = 1, ..., n bestimmt werden. Hierzu müssen die Koeffizienten in die Polynomapproximation eingesetzt, das Integral berechnet und dann die k = 1,..., n Gleichungen nach den Koeffizienten aufgelöst werden. Dieses Verfahren benötigt bereits n gemessene Winkelinkremente Δφ^k von Zeitintervallen, die zusammen den Zeitbereich des kumulierten Drehwinkels Δθ^k bilden, und liefert den kumulierten Drehwinkel $\Delta {\theta }_{out}^I$

  

  des Intervalls [t_I-1, t_I].
• M. B. Ignagni. „Efficient dass of optimized coning compensation algorithms." In: Journal of Guidance, Control, and Dynamics 19.2 (1996), pp. 424-429. DOI: 10.2514/3.21635. eprint: https://doi.org/10.2514/3.21635 beschreibt einen Algorithmus höherer Ordnung. Allerdings steigt für höhere Ordnungen n der Rechenaufwand stark an und die Approximation der Messdaten durch Polynome hohen Grades n - 1 führt zu starken Oszillationen an den Rändern des Interpolationsintervalls, die als Runges Phänomen bekannt sind.
• R.A. McKern. „A study of transformation algorithms for use in a digital computer." MA thesis. Massachusetts Institute of Technology, Dept. of Aeronautics and Astronautics, 1968. URL: http://hdl.handle.net/1721.1/14164 offenbart einen Algorithmus zur Kombination der relativen Rotationen aus aufeinanderfolgenden Zeitintervallen. Dieser Algorithmus berechnet über den Umweg der Einheitsquaternionen aus den Winkelinkrementen zweier zeitlich aufeinanderfolgender Zeitintervalle das Einheitsquaternion der relativen Lageänderung während des gesamten Zeitintervall. Dabei wird der zuvor beschriebene Algorithmus zur Berechnung des kumulierten Drehwinkels $\Delta {\theta }_{out}^I$

  

  für das Intervall [t_I-1, t_I] genutzt, allerdings nur für n = 2 aufeinanderfolgende Winkelinkemente. Dann wird mit dem erhaltenen Drehwinkel $\Delta {\theta }_{out}^I$

  

  das Einheitsquaternion der relativen Lageänderung in dem Zeitintervall [t_I-1, t_I] berechnet zu: $$\begin{array}{l}{}_{I-1}^{I}q\\ \begin{array}{cc}& =e^{\frac{\overrightarrow{\Delta }{\theta }_{out}^I}{2}}=cos\frac{\Delta {\theta }_{out}^I}{2}+\frac{\overrightarrow{\Delta }{\theta }_{out}^I}{\Delta {\theta }_{out}^I}sin\frac{\Delta {\theta }_{out}^I}{2}\cong 1-\frac{{\left(\Delta {\theta }_{out}^I\right)}^2}{8}+\frac{\overrightarrow{\Delta }{\theta }_{out}^I}{2}\end{array}\\ \begin{array}{cc}& \left(1-\frac{{\left(\Delta {\theta }_{out}^I\right)}^2}{24}\right),\end{array}\\ mit\Delta {\theta }_{out}^I={\Vert \Delta {\theta }_{out}^I\Vert }_2.\end{array}$$

  
• Der Drehwinkel-Vektor im Argument der Exponentialfunktion und auf der rechten Seite der Gleichung ist als sogenanntes reines Quaternion zu interpretieren, das nur einen mit den Vektorkomponenten des Winkelinkrements $\mathrm{\Delta }{\theta }_{out}^I$

  

  gebildeten Imaginärteil besitzt. Darüber hinaus werden die trigonometrischen Funktionen durch ihre Taylorpolynome bis zur 3. Ordnung approximiert.
• Das Ergebnis wird dann verwendet, um ausgehend von einem Startzeitpunkt t₀ und dem Quaternion ${}_0^{0}q=1$

  

  das Quaternion der relativen Lageänderung rekursiv zu berechnen gemäß: $${}_0^{I}q{=}_0^{I-1}{q}_{I-1}^{I}q.$$

  

  Hierbei sind die Quaternionen mittels des Quaternionproduktes zu multiplizieren.
• In Analogie zu der Formel (2.3) werden die zu den Winkelinkrementen Δθ^k-1, Δθ^k gehörenden Quaternionen ${}_{k-2}^{k-1}q,{}_{k-1}^{k}q$

  

  berechnet und dann in Analogie zur Formel (2.4) multipliziert, um das Quaternion $${}_{k-2}^{k}q{=}_{k-2}^{k-1}{q}_{k-1}^{k}q.$$

  

  zu erhalten. Anschließend ist dann das Winkelinkrement der kombinierten Drehungen durch Inversion der Funktion (2.3) zurückzugewinnen. Die entsprechende Relation lautet für das aus den Vektorkomponenten des Winkelinkrements $\mathrm{\Delta }{\theta }_{out}^I$

  

  gebildete Quaternion wie folgt: $$\begin{array}{l}\overrightarrow{\mathrm{\Delta }}{\theta }_{out}^I=\frac{arcsinq}{q}{}_{k-2}^k\overrightarrow{q}\cong \left(1+\frac{q^2}{6}+\frac{3q^4}{40}+\frac{5q^6}{112}+\frac{35q^8}{1152}+\frac{63q^{10}}{2816}\right){}_{k-2}^k\overrightarrow{q},\hfill \\ mitq=\Vert {}_{k-2}^k\overrightarrow{q}\Vert .\hfill \end{array}$$

  
• Diese Lösung für den Spezialfall von n=2 zeitlich aufeinanderfolgenden Winkelinkrementen Δθ^k-1, Δθ^k kann auf beliebige Anzahl n von Winkelinkrementen für einen Zeitbereich und einen zugehörigen Drehwinkel $\mathrm{\Delta }{\theta }_{out}^I$

  

  verallgemeinert werden.
• Die iterative Berechnung des Drehwinkels $\mathrm{\Delta }{\theta }_{out}^I=\mathrm{\Delta }{\theta }_{out}^{I,n}$

  

  des $\left[t_{i-1},t_{i-1}+{\displaystyle {\sum }_{k=1}^n\mathrm{\Delta }{t}_k}\right]$

  

  kann dann durch die Prozedur beschrieben werden:
• Prozedur $\mathrm{\Delta }{\mathrm{\theta }}_{out}^I\left(\left(\mathrm{\Delta }{\mathrm{\theta }}^1,\dots ,\mathrm{\Delta }{\mathrm{\theta }}^n\right)\right)$
  
• Setze ${}_0^{0}q=1$
  
• Für k = 1,...,n
  • Berechne ${}_{k-1}^{k}q$
    
    aus Δθ^k mit der Gleichung (2.3)
  • Berechne ${}_0^{k}q$
    
    aus ${}_0^{k-1}q$
    
    und ${}_{k-1}^{k}q$
    
    mit der Gleichung (2.5)
• Berechne $\mathrm{\Delta }{\theta }_{out}^I=\mathrm{\Delta }{\theta }_{out}^{I,n}$
  
  aus ${}_0^{n}q$
  
  mit der Gleichung (2.6)
• Die Berechnung des Drehwinkelinkrementes $\mathrm{\Delta }{\theta }_{out}^{I,n}$

  

  erfolgt indirekt über Quaternionen und beruht daher auf den drei Gleichungen (2.3), (2.4) und (2.5). Dies erfordert die Berechnung verschiedener transzendenter Funktionen (cos, sin, arcsin) z.B. nach Approximation durch ihre Taylorpolynome, wodurch der Rechenaufwand erhöht wird. Nach jeder n-ten rekursiven Aktualisierung des Quaternions wird die Gleichung (2.6) ausgewertet, um aus dem Quaternion das iterative Drehwinkelinkrement $\mathrm{\Delta }{\theta }_{out}^{I,n}$

  

  zurückzugewinnen. Auch dies ist mit erhöhtem Rechenaufwand verbunden, der jeden n-ten Rekursionsschritt anfällt.
• Bei Approximation der trigonometrischen Funktionen in der Gleichung (2.3) ist das Quaternion kein Einheitsquaternion mehr. Die Relation der Gleichung (2.6) gilt jedoch nur für Einheitsquaternionen. Ihre Anwendung auf nicht normierte Quaternionen führt zu einem zusätzlichen (numerischen) Fehler.
• Die Berechnung der Gleichung (2.6) bei hinreichend großen Winkelinkrementen $\mathrm{\Delta }{\theta }_{out}^I,$

  

  wie sie z.B. entstehen, wenn eine große Anzahl n von Inkrementen Δθ^k summiert wird, erfordert bei Approximation durch ihr Taylorpolynom eine Berücksichtigung recht hoher Ordnungen, um die gewünschte numerische Genauigkeit zu erreichen. Dieses erhöht den Rechenaufwand.
• Die bekannten Verfahren zur Bestimmung der Lage aus einem Drehwinkel $\mathrm{\Delta }{\theta }_{out}^I$

  

  eines Zeitbereichs, bspw. des Intervalls [t_I-1, t_I] durch Zusammenfassung einer Anzahl von n Winkelinkrementen Δφ^k an Zeitpunkten t_meas,k, bspw. k = 1,...,n, die zusammen den Zeitbereich bilden, sind sehr rechenaufwändig und benötigen daher relativ viel Rechenzeit, Rechenkapazität und erfordern eine relativ leistungsfähige und damit aufwändige Elektronik mit Hardware-Rechenleistung.
• DE 10 2022 121 662 B3 offenbart ein Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehwinkel des Objektes, wobei Drehwinkelinkremente für aufeinanderfolgende Messzeitpunkte gemessen und zu einem kumulierten Drehwinkel für einen die Zeitpunkte beinhaltenden Kumulations-Zeitbereich zusammengefasst werden. Es erfolgt ein Berechnen des kumulierten Drehwinkels für einen Kumulations-Zeitbereich als Funktion von Exponenten einer Approximation mindestens zweiter Ordnung der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3), wobei iterativ kumulierte Drehwinkelinkremente für die Zeitpunkte jeweils aus der Funktion der Exponenten mit jeweils einem gemessenen Drehwinkelinkrement für das jeweilige Zeitintervall und dem für den vorhergehenden Zeitpunkt bestimmten kumulierten Drehwinkelinkrement oder einem vorgegebenen Ausgangs- Drehwinkelinkrement als Eingangsgrößen bestimmt werden und der kumulierte Drehwinkel für den Zeitbereich durch das am letzten Messzeitpunkt bestimmte letzte kumulierte Drehwinkelinkrement bestimmt ist. Die Berechnung der iterativen Drehwinkelinkremente erfolgt direkt ohne Umwege über Quaternionen.
• US 2018/0231385 A1 offenbart ein Verfahren zum Bestimmen der Posenparameter eines Inertialmesseinheitssensors (IMU) mit den Schritten des Sammelns von durch IMU-Sensoren erzeugten Messdaten, der Verwendung eines Prozessors zum zeitlichen Integrieren der Messdaten, einschließlich etwaiger Fehler, des Erzeugens eines zeitlich kontinuierlichen Fehlerausbreitungsmodells und des zeitlichen Integrierens. Mit dem Fehlerausbreitungsmodell lassen sich Kompensationsgradienten für die Posenparameter generieren.
• EP 1 637 840 A1 offenbart ein Verfahren zur Kompensation der Konizität in einem Strapdown-Inertialnavigationssystem, das Gruppen von fünf aufeinanderfolgenden inkrementellen Drehwinkeln eines körperfesten Koordinatensystems verwendet, die von orthogonal montierten Kreiseln in regelmäßigen Messintervallen gemessen werden. Jede Gruppe von fünf Messungen wird während eines Gruppenintervalls erhalten, das fünf Messintervallen entspricht. Die kegelkompensierte Winkelverschiebung des körperfesten Koordinatensystems um eine feste Achse im Raum während eines p-ten Gruppenintervalls wird durch Summieren der fünf gemessenen Inkrementalwinkel und eines Kegelkompensationsterms erhalten. Der Konuskompensationsterm besteht aus der Summe von der Hälfte des Kreuzprodukts einer ersten und einer zweiten Vektorsumme und der gewichteten Summe von drei Vektorkreuzprodukten. Die zweite Vektorsumme ist die Summe der fünf inkrementellen Drehwinkel in einer Gruppe. Die erste Vektorsumme ist die Summe der zweiten Vektorsumme über p-Gruppen. Der Multiplikator und der Multiplikand jedes Vektorkreuzprodukts ist eine gewichtete Summe von fünf gemessenen Inkrementalwinkeln. Die kegelkompensierte Winkelverschiebung kann über p-Gruppen summiert werden, um eine genaue Schätzung des Vektordrehwinkels über mehrere Gruppenintervalle zu erhalten.
• Ausgehend hiervon ist es Aufgabe der vorliegenden Erfindung, ein verbessertes Verfahren zu schaffen.
• Die Aufgabe wird mit dem Verfahren mit den Merkmalen des Anspruchs 1 und die Messeinrichtung mit den Merkmalen des Anspruchs 11 und 12 gelöst. Vorteilhafte Ausführungsformen sind in den Unteransprüchen beschrieben.
• Es wird vorgeschlagen, dass zum Bestimmen der Lage eines Objektes kumulierte Drehwinkelinkremente θ^k für die Messzeitpunkte t_k iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθ^k, die an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ^k-1 und der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

  

  mit $\zeta =\sqrt{\left({\theta }^{k-1}\text{/}2\right)\cdot \left({\theta }^{k-1}\text{/}2\right)}$

  

  als Funktion des kumulierten Drehwinkelinkrementes θ^k-1 zum jeweils vorhergehenden Messzeitpunkt t_k-1 berechnet werden. Hierzu erfolgt ein Approximieren der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

  

  durch Padé-Approximationen.
• Die Drehwinkelinkremente und Drehwinkel enthalten Winkelinformationen für die erforderlichen Raumrichtungen und können Vektoren mit den Vektorkomponenten der drei Einheiten i, j, k des oben beschriebenen Einheitsquaternions bzw. der drei Raumrichtungen x, y, z haben. Auch wenn die Drehwinkelinkremente und Drehwinkel in der vorliegenden Anmeldung und den Ansprüchen zur Vereinfachung nicht explizit als Vektoren bezeichnet sind, ist die Ausführungsform als Vektor oder Darstellungen durch schiefsymmetrische Matrizen davon umfasst.
• Das Verfahren hat den Vorteil, dass die Berechnung der iterativen Drehwinkelinkremente direkt ohne Umwege über Quaternionen erfolgt. Dabei ist nur eine transzendente Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

  

  mit ζ als Funktion des kumulierten Drehwinkelinkrementes θ^k-1 zum jeweils vorhergehenden Messzeitpunkt t_k-1 zu berechnen. Dies erfolgt durch Pade-Approximationen der transzendenten Funktion. Auf diese Weise kann die einzige zu berechnende transzendente Funktion mit hoher Präzision durch eine rationale Funktion approximiert und somit effizient und mit vorab bekannter Anzahl an Rechenoperationen berechnet werden.
• Pade-Approximationen selbst sind bekannte Werkzeuge für numerische Approximationen. Die Approximationen der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

  

  $$

  

  mit ζ als Funktion des kumulierten Drehwinkelinkrementes θ^k-1 durch Pade-Approximationen nutzt die Tatsache aus, dass die Funktion B(ζ) spezielle Eigenschaften besitzt, durch die ihre Pade-Approximationen erst so genau und damit effizient werden. Damit lässt sich auch der maximale Approximationsfehler abschätzen.
• Zudem werden nur gerade Potenzen des Argumentes $\zeta =\sqrt{\left({\theta }^{k-1}\text{/}2\right)\cdot \left({\theta }^{k-1}\text{/}2\right)}$

  

  benötigt, so dass die Berechnung der Quadratwurzel entfällt. Die Approximationen der transzendenten Funktion sind nur für den Bereich $0\le \zeta \le \frac{\pi }{2}$

  

  notwendig, da sich die Funktion für den doppelten Winkel daraus konstruieren lässt. Damit lässt sich auch direkt die Singularität der Funktion bei ζ = π vermeiden, wo die Approximation zusammenbricht.
• Zusätzlich zu relativen kumulierten Winkelinkrementen θ^k, deren Norm zur Gewährleistung einer vorgegebenen Präzision begrenzt sein muss, können auch absolute Winkelinkremente und somit insbesondere auch der die Lage des Objektes (Körpers) parametrisierende Drehwinkel mit vorgegebener Präzision berechnet werden. Die einzige dabei zu berechnende transzendente Funktion S(C;) wird dennoch mit hoher Präzision durch eine rationale Funktion approximiert und kann somit effizient und mit vorab bekannter Anzahl an Rechenoperationen berechnet werden. Zudem bietet sich die Möglichkeit, die trigonometrischen Funktionen (Sinus und Cosinus) der kumulierten Drehwinkelinkremente θ^k und daher auch die zugehörige Rotationsmatrix ohne weitere Approximation oder Berechnung von transzendenten Funktionen aus dem Ergebnis für die transzendente Funktion S(C;) zu berechnen. Damit lassen sich auch mit reduzierter Prozessorleistung sehr präzise und schnell trigonometrische Funktionen und Drehungen automatisiert rechnergestützt allgemein aus Winkeln θ^k aus dem Ergebnis der transzendenten Funktion B(ζ) bestimmen.
• Damit werden folgende technische Probleme gelöst:
• 1) Es können statt relativer Winkelinkremente auch absolute Winkelinkremente an einen Bordcomputer weitergegeben werden. Dadurch kann bei einem temporären Ausfall der Kommunikation zwischen IMU und Bordcomputer und somit dem Verlust auch nur einzelner Winkelinkremente gewährleistet werden, dass sich dieser Ausfall / Fehler auch ohne andere Informationsquellen (Sensoren) kompensieren lässt. Damit ist eine Lagebestimmung aus den Messdaten der Gyroskope direkt innerhalb der IMU, unabhängig vom Bordomputer möglich.
• 2) Die Notwendigkeit, Wertetabellen (eng. lookup tables) und/oder zeitverzögernde Iterationsalgorithmen (mit ggf. einer Genauigkeitsüberprüfung als Abbruchbedingung und damit einer nicht festgelegen Anzahl an Rechenoperationen), für die Berechnung der transzendenten Funktionen zu verwenden, entfällt. Stattdessen kann eine direkte Berechnung der erforderlichen transzendenten Funktion mit vorab festgelegter Anzahl von Rechenoperationen bei garantierter Genauigkeit erfolgen.
• 3) Mit Hilfe der Rotationsmatrix können aus den durch die Beschleunigungssensoren der IMU gemessenen Beschleunigungen / relativen Geschwindigkeitsinkremente die translatorische Geschwindigkeit und die Position berechnet werden. In Kombination mit der Lage und mit einem Modell für die Gravitation erscheint somit eine vollständige Trägheitsnavigation (eng. inertial navigation) innerhalb der IMU möglich, unabhängig vom Bordcomputer.
• 4) trigonometrische Funktionen und Rotationsmatrizen können auch in anderen Anwendungsgebieten wie beispielsweise Computergrafik / Computer Vision, CAD bei gegebenen Drehwinkeln effizient und präzise berechnet werden.
• Die Berechnung der kumulierten Drehwinkelinkremente θ^k für einen jeweiligen Messzeitpunkt t_k kann als Funktion des in mindestens zweiter Ordnung im zweiten Argument approximierten Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) erfolgen. Dabei werden iterativ kumulierte Drehwinkelinkremente θ^k für die Messzeitpunkte t_k mit k = 1,..., n jeweils aus der Funktion des approximierten Exponenten der Baker-Hausdorff-Formel mit jeweils einem zum Messzeitpunkt t_k gemessenen relativen Drehwinkelinkrement Δθ^k und dem für den vorhergehenden Zeitpunkt t_k-1 bestimmten kumulierten Drehwinkelinkrement θ^k-1 oder einem vorgegebenen Start-Drehwinkelinkrement θ⁰ als Eingangsgrößen bestimmt.
• Damit ist eine präzise Summation / Kumulation der Winkelinkremente auch im Fall einer zeitlich nicht konstanten Drehachse auf sehr effiziente Weise mit reduziertem Rechenaufwand möglich. Liegt zum Abtastzeitpunkt t_k ein neues (gemessenes) relatives Winkelinkrement vor, welches abkürzend mit Δθ^k bezeichnet wird, so wird dieses durch $${\theta }^k=2\gamma \left(\frac{{\theta }^{k-1}}{2},\frac{\Delta {\theta }^k}{2}\right)$$

  

  zu dem bestehenden kumulierten Winkelinkrement θ^k-1 hinzugefügt und ergibt dann das neue kumulierte Winkelinkrement θ^k. Dieser Ablauf wird nachfolgend in abkürzender Schreibweise als Algorithmus 1 dargestellt, wobei das Zwischenergebnis mit ${\theta }_{out}^{I,k}$

  

  bezeichnet ist, welches vom Zeitpunkt t_I-1 ausgehend, eine ganze Anzahl von k relativen Winkelinkrementen θ^k berücksichtigt. Das gespeicherte Winkelinkrement ${\theta }_{out}^{I,k}$

  

  wird nach einer Anzahl von n Updates wieder auf Null gesetzt.
• Damit wird ein kumuliertes, aber immer noch relatives (also auf den jeweiligen Startzeitpunkt t_I-1 der Kumulation bezogenes) Winkelinkrement ${\theta }_{out}^{I,k}$

  

  erzeugt, welches die aus der nacheinander ausgeführten Anzahl von n relativen Drehungen mit Winkeln Δθ^k, k = 1, ...,n resultierende Gesamtdrehung parametrisiert.

  

  
• Die Berechnung der Funktion γ in Zeile 4 des Algorithmus 1 erfordert hierbei insbesondere die lediglich approximative Berechnung nur einer einzigen transzendentalen Funktion $$B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$$

  

  die mit dem Wert der halben quadrierten Norm des letzten kumulierten Winkelinkrementes, also $\zeta =\frac{{\theta }^{k-1}}{2}$

  

  mit ${\theta }^{k-1}=\sqrt{{\theta }^{k-1}\cdot {\theta }^{k-1}}$

  

  auszuwerten ist (vgl. die expliziten Ausdrücke in Gleichung (3.5a) weiter unten). Hierbei ist die Funktion B(ζ) gerade, so dass nur die quadrierte Norm (θ^k-1)² benötigt wird und damit keine Berechnung der Quadratwurzel erforderlich ist. Darüber hinaus basiert der Algorithmus auf einer einfacheren Iteration und sollte die Analyse der Fortpflanzung von Messunsicherheiten und Messrauschen gegenüber dem zum Stand der Technik einleitend beschriebenen quaternionenbasierten Algorithmus vereinfachen. Diese Eigenschaften erlauben eine direkte Implementierung auf der Hardware (zum Beispiel einem Field Programmable Gate Array (FPGA)) der IMU. Der Algorithmus kann damit als Teil eines Sensorpaketes und der Signalverarbeitung der Sensordaten für die Ausgabe an der Benutzerschnittstelle implementiert werden. Er berechnet Winkelinkremente mit reduzierter Frequenz, die im Bordcomputer dann für die Berechnung der Lage verwendet werden können.
• Die Pade-Approximation erlaubt eine besonders effiziente Berechnung der transzendenten Funktion im Zusammenhang mit dem iterativen Algorithmus 1.
• Dabei können die Drehwinkelinkremente Δθ^k an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen (bzw. aus gemessenen Δφ^k berechnet) und zu einem kumulierten Drehwinkel $\Delta {\theta }_{out}^I$

  

  für einen die Messzeitpunkte t_k umfassenden Kumulations-Zeitbereich I zusammengefasst werden, indem jeweils ein kumulierter Drehwinkel $\Delta {\theta }_{out}^I$

  

  für einen Kumulations-Zeitbereich I als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) berechnet wird. Dabei werden iterativ kumulierte Drehwinkelinkremente $\Delta {\theta }_{out}^{I,k}$

  

  für die Zeitpunkte t_k mit k = 1,..., n jeweils aus der Funktion der Exponenten mit jeweils einem Drehwinkelinkrement Δθ^k des Messzeitpunktes t_k und dem für das vorhergehende Zeitintervall [t_i-1, t_i-1+k] bestimmten kumulierten Drehwinkelinkrement $\Delta {\theta }_{out}^{I,k-1}$

  

  oder einem vorgegebenen Ausgangs- Drehwinkelinkrement $\Delta {\theta }_{out}^{I\mathrm{,0}}$

  

  als Eingangsgrößen bestimmt. Der kumulierte Drehwinkel $\Delta {\theta }_{out}^I$

  

  für den Zeitbereich I ergibt sich als Funktion des letzten kumulierten Drehwinkelinkrements $\Delta {\theta }_{out}^{I,n},$

  

  das für den letzten Zeitpunkt t_n des Intervalls des Zeitbereichs I bestimmt wurde. Es kann diesem kumulierten Drehwinkelinkrement $\Delta {\theta }_{out}^{I,n}$

  

  entsprechen.
• Für die Iteration kann das Ausgangs-Drehwinkelinkrement $\Delta {\theta }_{out}^{I,k=0}$

  

  für k = 0 bei einer Initialisierung mit dem Wert Null vorgegeben werden. Zur Bestimmung des kumulierten Drehwinkels $\Delta {\theta }_{out}^I$

  

  für den Zeitbereich [t_I-1, t_I] mit den gemessenen Drehwinkelinkrementen Δθ^k für die jeweiligen Zeitpunkte t_k mit k = 1,..., n als Eingangsgrößen wird vor dem ersten Iterationsschritt die Initialisierung durchgeführt.
• Anschließend werden die kumulierten Drehwinkelinkremente $\Delta {\theta }_{out}^{I,k}$

  

  als Zwischenergebnisse iterativ mit k = 1 bis n aus dem jeweils für den vorhergehenden Iterationsschritt k-1 bestimmten oder vorgegebenen kumulierten Drehwinkelinkrement $\Delta {\theta }_{out}^{I,k-1}$

  

  und dem für den Zeitpunkt t_k, der zum jeweiligen Iterationsschritt k gehörig ist, gemessenen Drehwinkelinkrement $\Delta {\theta }_{out}^{I,k}$

  

  als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) berechnet. Dieser Exponent muss hierzu nur in seinem zweiten Argument durch sein Taylorpolynom approximiert werden, da der erste Exponent exakt ist.
• In jedem Iterationsschritt ist die Auswertung nur einer Gleichung, die später als Gleichung (3.5) näher erläutert wird, erforderlich. Dies hat die folgenden Konsequenzen:
• - Es ist nur eine einzige transzendente Funktion, $B\left(\zeta \right)=B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)$
  
  zu berechnen, die effizient durch die Padé-Approximation approximiert wird.
• - Die Iteration ist vereinfacht, da (bis auf die in jedem Fall notwendige Initialisierung) keine Ausführung eines zusätzlichen Rechenschrittes bei jedem n-ten Zeitschritt erforderlich ist. Dies ist insbesondere für die Implementierung auf stark limitierter Hardware (wie beispielsweise einem Field Programmable Gate Array (FPGA)) ein wichtiger Vorteil.
• - Die Analyse der Fortpflanzung von Messunsicherheiten und Messrauschen wird erleichtert.
• - Die Padé-Approximation der einzigen transzendente Funktion, $B\left(\frac{{\theta }^{k-1}}{2}\right)$
  
  ist genauer als die Approximation der Gleichung (1.3) durch ihr Taylorpolynom gleicher Ordnung. Für große kumulierte Drehwinkel $\mathrm{\Delta }{\theta }_{out}^I$
  
  und damit insbesondere für eine größere Anzahl n wird eine größere Genauigkeit erreicht.
• - Durch eine getrennte Wahl der beiden Ordnungen für die Approximation lässt sich das Verfahren recht einfach an die gegebenen Anforderungen (bspw. die maximale Drehrate, Sensorfrequenz, Teiler zur Frequenz, Anzahl der zu summierenden Winkelinkremente) für den Benutzer anpassen.
• Im Fall einer zeitlich veränderlichen Drehachse ergeben sich nun folgende Änderungen (im vorgeschlagenen Verfahren berücksichtigt):
1. 1. Das von den Sensoren gemessene Integral Δφ^k der Drehrate ω ergibt nicht direkt das Winkelinkrement Δθ^k, das direkt Drehachse und den Drehwinkel der im Abtastintervall erfolgten Drehung angibt. Stattdessen ist Δθ^k aus den Integralen der aktuellen und vergangenen Drehraten zu berechnen.
2. 2. Die direkte Summation zeitlich aufeinanderfolgender mit den Zeitintervallen [t_meas,k-2, t_meas,k-1] [t_meas,k-1, t_meas,k] assoziierter Winkelinkremente Δθ^k-1 , Δθ^k ergibt nicht das zur Drehung im Zeitintervall [t_meas,k-2, t_meas,k] gehörende Winkelinkrement. Anstelle der einfachen Summation der Winkelinkremente ist eine kompliziertere mathematische Operation auszuführen.
• Das Verfahren beruht auf den bei zeitlich veränderlichen Drehachsen erforderlichen komplexen mathematischen Operationen anstelle einer einfachen Summation der Winkelinkremente. Die Ursache hierfür ist, dass Drehungen um nichtparallele Drehachsen miteinander nicht vertauschen. Vielmehr bilden die Drehungen im (dreidimensionalen) Raum eine Gruppe, die Drehgruppe SO(3), welche eine sogenannte Lie-Gruppe ist. Das mathematische Gebiet der Lie-Gruppen bildet den theoretischen Rahmen, aus welchem sich die Darstellungen der Gruppenelemente in Form von Rotationsmatrizen oder in Form der Einheitsquaternionen ergeben.
• Die Rotationsmatrizen sind spezielle (ihre Determinante ist 1) orthogonale (ihr Transponiertes ist ihr eigenes Inverses) Matrizen in drei Dimensionen. Der Name der Gruppe, SO(3), leitet sich aus den Begriffen Spezielle-Orthogonale Rotationsmatrizen in 3 Dimensionen ab.
• Einheitsquaternionen sind die Erweiterung komplexer Zahlen auf drei unabhängige imaginäre Einheiten, deren aus Real- und Imaginärteil berechnete Norm 1 ist und mittels eines nichtkommutativen, normerhaltenden Produktes multipliziert werden. Eine Rotation kann durch ein Einheitsquaternion dargestellt werden, wobei die Darstellung nicht eindeutig ist, sondern sowohl q als auch -q zur Darstellung einer konkreten Drehung verwendet werden können.
• In Abhängigkeit der gewählten Darstellung (in Form von Vektoren, Einheitsquaternionen, ...) folgen dann auch die oben beschriebenen Operationen für die Berechnung der Winkelinkremente und deren Summation.
• Der kumulierte Drehwinkel $\mathrm{\Delta }{\mathrm{\theta }}_{out}^I$

  

  für einen Zeitbereich I, der eine Anzahl von n aufeinanderfolgenden Zeitpunkten t_meas,k mit k = 1, ..., n beinhaltet, kann hierzu mit der Funktion des Exponenten $\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  rekursiv (d.h. iterativ) aus einer Anzahl n von Drehwinkelinkrementen Δθ^k, die an aufeinanderfolgenden Messzeitpunkten t_meas,k mit k = 1, ..., n gemessen wurden, und dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ^k-1 berechnet werden. Der Exponent nullter Ordnung kann dabei ${\mathrm{\gamma }}^{\left(0\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{{\mathrm{\theta }}^{k-1}}{2}$

  

  $$

  

  sein. Der Exponent erster Ordnung kann dabei ${\mathrm{\gamma }}^{\left(1\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}+\frac{{\mathrm{\theta }}^{k-1}}{2}\times \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}+\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^2\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  $$

  

  und der Exponent zweiter Ordnung ${\mathrm{\gamma }}^{\left(2\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{1}{2}{\left(\frac{2}{{\mathrm{\theta }}^{k-1}}\right)}^2$

  

  $\left(1-3B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(1-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^{2}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\right)\right)\left({\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^2{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\right)\frac{{\mathrm{\theta }}^{k-1}}{2}+\text{B}$

  

  $\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(1-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^{2}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\right)\left(\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}-{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\frac{{\mathrm{\theta }}^{k-1}}{2}\right)+\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot$

  

  $\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\frac{{\mathrm{\theta }}^{k-1}}{2}\times \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}$

  

  sein. Hierzu wird die transzendente Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)\cong {\tilde{B}}_{Padé,N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\displaystyle {\sum }_{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m{\zeta }^{2m}}}{1+{\displaystyle {\sum }_{n=1}^{\left[\frac{N_0+1}{2}\right]}{b}_n{\zeta }^{2n}}},$

  

  $$

  

  wobei N₀ ein vorgegebener ganzzahliger Wert und die Parameter a_m und b_m von dem Wert N₀ abhängig vorgegebene Koeffizienten der Pade-Approximation sind, und sich das für den Zeitpunkt t_k kumulierte Drehwinkelinkrement $\frac{{\mathrm{\theta }}^k}{2}=\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  jeweils aus der Summe der Entwicklungskoeffizienten ${\mathrm{\gamma }}^{\left({\text{m}}_{\alpha }\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  des Exponenten der Baker-Campbell-Hausdorff-Formel in der nullten, ersten und mindestens zweiten Ordnung m_α = 0, 1 und 2 ergibt.
• Es wird für den Messzeitpunkt t_meas,k ein kumuliertes Drehwinkelinkrement $\frac{\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k}}{2}=\gamma \left(\frac{\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  $$

  

  jeweils aus der Summe der Entwicklungskoeffizienten γ^(mα) $\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  des Exponenten der Baker-Campbell-Hausdorff-Formel in der nullten, ersten und mindestens zweiten Ordnung m_α = 0, 1 und 2 bestimmt.
• Die Berechnung beruht auf nur einer einzigen transzendenten Funktion, was den Zeit- und Rechenaufwand reduziert.
• Die Funktion des Exponenten $\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

  

  kanndie Summe der Argumente der Approximationen in der nullten, ersten, zweiten und dritten Ordnung m_α = 0, 1, 2 und 3 umfassen. Durch die Berücksichtigung auch eines Entwicklungskoeffizienten der höheren dritten Ordnung kann die Genauigkeit weiter gesteigert werden. Auch hierbei ist nur die eine einzige transzendente Funktion enthalten, die durch die Pade-Approximation effizient und genau berechnet werden kann.
• Es kann weiterhin die (automatisierte, rechnergestützte) Berechnung einer Rotationsmatrix der aktiven Rotation zwischen einem ersten Koordinatensystem F des Objektes und einem nach Anwendung einer durch den kumulierten Drehwinkel θ^k beschriebenen Drehung erhaltenen zweiten Koordinatensystem G mit den Approximation von trigonometrischen Funktionen enthaltenden Ergebnissen der Pade-Approximationen und Berechnung der translatorischen Geschwindigkeit und/oder Position des Objektes aus den ebenfalls gemessenen Beschleunigungen oder relativen Geschwindigkeitsinkrementen unter Verwendung eben dieser Rotationsmatrix erfolgen.
• Über die Verarbeitung von Sensordaten einer IMU und der effizienten Berechnung der Lage und Position von Objekten hinaus kann eine effiziente Berechnung trigonometrischer Funktionen von gegebenen Winkeln mit reduzierter Rechenleistung und hoher Genauigkeit automatisiert rechnergestützt realisiert werden. Dies kann in vielfältige Anwendungsmöglichkeiten in verschiedenen Bereichen eingesetzt werden, z.B. in Computer Vision.
• Aus einem gegebenen Winkel θ^k können durch die Padé-Approximation der Funktion B(ζ) effizient und insbesondere ohne Approximation irgendwelcher weiterer transzendenter Funktionen trigonometrische Funktionen, wie beispielsweise die Funktionen $\frac{sin2\zeta }{2\zeta }$

  

  und cos2ζ, erhalten werden.
• Die Erfindung wird nachfolgend anhand der beigefügten Zeichnung mit einem Ausführungsbeispiel näher erläutert. Es zeigt:
• 1 - Zeitverlauf der k = 1,...,n gemessenen bzw. des aus den Messungen zu den Messzeitpunkten t_meas,k ermittelten Inputs Δθ^k, dem k-ten Iterationsschritt zur Bestimmung des k-ten kumulierten Winkelinkrements $\Delta {\theta }_{out}^{I,k}$
  
  und dem Output $\Delta {\theta }_{out}^{I,n}$
  
  zum Zeitpunkt t_I;
• 2 - Skizze einer Messeinrichtung zur Bestimmung der Lage eines Objektes;
• 3 - Flussdiagramm des Verfahrens zur Bestimmung eines kumulierten Drehwinkelinkrementes aus gemessenen Drehwinkelinkrementen (der Einfachheit halber für konstante Abtastzeitschritte Δt_k = Δt);
• 4 - Diagramme der relativen Fehler bei der Verwendung der Pade-Approximation zur Berechnung der transzendenten Funktion B(ζ);
• 5 - Diagramme der relativen Fehler bei direkter Verwendung der Pade-Approximation von B(ζ) und eines regularisierten Ausdruckes zur Berechnung der transzendenten Funktion B⁽²⁾(ζ);
• 6 - Relativer Fehler der Approximation von sin_N0 (θ) und cos_N0 (θ), die auf einer Padé-Approximation von B(ζ) basieren;
• 7 - Relativer Fehler der Approximation durch die Taylor-Polynome mit den Graden 2n+1, 2n von sin(θ) und cos(θ);
• 8 - Relativer Fehler der Approximationen von sin θ und cos θ durch die Taylor-Polynome mit den Graden 2n+1, 2n von $sin\left(\frac{\theta }{2}\right)$
  
  und $cos\left(\frac{\theta }{2}\right)$
  
  verwenden in den Ausdrücken für sin θ und cos θ mittels der Funktionen der halben Winkel;
• 9 - Weiteres Flussdiagram des Verfahrens.
• 1 zeigt einen Zeitverlauf der k = 1,...,n gemessenen bzw. des aus den Messungen zu den Messzeitpunkten t_k ermittelten Messwerte für die Drehwinkelinkremente, d.h. den Inputs Δθ^k, dem k-ten Iterationsschritt zur Bestimmung des k-ten kumulierten Winkelinkrements $\Delta {\theta }_{out}^{I,k}$

  

  und dem kumulierten Drehwinkel, d.h. dem Messergebnis als Output $\Delta {\theta }_{out}^{I,n}$

  

  zum Zeitpunkt t_I. Damit wird die Zeitabfolge mit den Messzeitpunkten t_meas,k, der Aufsummierung für unterschiedliche Zeitintervalle k = 1, ...., n und der Zusammenfassung zum niederfrequenten (Low Frequency) Messergebnis deutlich. Zur Vereinfachung der Darstellung wurden zeitlich konstante Abtastzeitschritte Δt_k = Δt angenommen.
• Die Fettschreibung der Größen bezeichnet Vektoren und die Normalschrift der Größen Skalare, wie z.B. die Norm des Drehwinkelvektors.
• 2 zeigt eine Skizze einer Messeinrichtung 1 zur Bestimmung der Lage eines Objektes 2. Die Messeinrichtung hat mindestens ein Inertiale Messeinheit IMU, die zur Messung der Drehwinkelinkremente Δθ^k aufeinanderfolgender Zeitpunkte t_k für die Erfassung rotierender (kreisender) Bewegungen der Inertialen Messeinheit IMU und des damit verbundenen Objektes 2 in den drei zueinander orthogonal stehenden Raumachsen (X- bzw. Y- bzw. Z-Achse) des kartesischen Koordinatensystems eingerichtet ist. Diese gemessenen Drehwinkelinkremente Δθ^k werden als Eingangswerte einer Datenverarbeitungseinheit zugeführt, die zur Berechnung eines kumulierten Drehwinkels $\Delta {\theta }_{out}^I=\mathrm{\Delta }{\theta }_{out}^{I,n}$

  

  für einen Kumulations-Zeitbereich I, der aufeinanderfolgende Zeitpunkte t_k mit i = 1,..., n umfasst, eingerichtet. Außerdem ist es optional denkbar, für die Berechnung der Geschwindigkeit und ggf. Position die Beschleunigungen oder relativen Geschwindigkeitsinkremente als Komponenten entlang der drei Raumachsen zu messen und der Datenverarbeitungseinheit zuzuführen.
• Die Datenverarbeitungseinheit 3 kann beispielsweise ein Mikroprozessor oder Mikrocontroller sein, der ein Computerprogramm mit Programmbefehlen ausführt, die den Mikroprozessor oder Mikrocontroller veranlassen, die Verfahrensschritte des anspruchsgemäßen Verfahrens auszuführen. Die Datenverarbeitungseinheit 3 kann aber auch als festverdrahtete Hardware-Logik, wie beispielsweise ein Field-Programmable-Gate-Array o.ä. ausgeführt sein.
• Mit dem Verfahren wird eine direkte Berechnung des kumulierten Drehwinkels $\Delta {\theta }_{out}^I$

  

  iterativ aus jeweils dem im letzten Schritt gespeicherten und dem aktuell gemessenen Drehwinkelinkrement Δθ^k durchgeführt. Die dabei ebenfalls in diesem Verfahren erhaltene Rotationsmatrix kann zur Berechnung der Geschwindigkeit und ggf. Position aus den ebenfalls gemessenen Beschleunigungen oder Geschwindigkeitsinkrementen verwendet werden. Grundlage für die Summation der gemessenen Drehwinkelinkremente bildet die aus der Theorie der Lie-Gruppen stammende Baker-Campbell-Hausdorff (BCH) Formel, die ein Vertauschungsgesetz für bestimmte lineare Operatoren angibt. Sie liefert für die Multiplikation zweier durch die Exponentialabbildung erhaltenen Elemente e^α, e^β einer Lie-Gruppe den Exponenten γ(α, β) des resultierenden Elementes e^{γ(α, β)}, welches also die Relation $${\text{e}}^{\gamma \left(\mathrm{\alpha },\mathrm{\beta }\right)}={\text{e}}^{\mathrm{\alpha }}{\text{e}}^{\mathrm{\beta }}$$

  

  erfüllt. Das Ergebnis γ(α, β) wird durch die Baker-Campbell-Hausdorff Formel als formale Reihenentwicklung in der Gesamtordnung m von α und β angegeben und hat daher schematisch folgende Form: $$\gamma \left(\alpha ,\beta \right)={\displaystyle {\sum }_{m=1}^{\infty }{\displaystyle \sum m_{\alpha }+m_{\beta }}=mC\left(m_{\alpha },m_{\beta }\right){\alpha }^{m_{\alpha }}}{\beta }^{m_{\beta }}$$

  

  mit dem Koeffizienten c(m_α, m_β), wobei die tatsächliche Form wegen der allgemeinen Nichtvertauschbarkeit der Gruppenelemente α β ≠ β α wesentlich komplizierter ist.
• Die allgemeine Baker-Campbell-Hausdorff (BCH) Formel ist als formale Reihe noch nicht in einer Form, die für die Summation zeitlich aufeinanderfolgender mit den Zeitintervallen [t_i-3+k, t_i-2+k], [t_i-2+k, t_i-1+k] assoziierter Winkelinkremente Δθ^k-1, Δθ^k geeignet ist, um das tatsächliche zur Drehung im Zeitintervall [t_i-3+k, t_i-1+k] gehörende Drehwinkelinkrement, d.h. den Drehwinkel $\Delta {\theta }_{out}^I$

  

  für den Zeitbereich I (in diesem Fall bestehend aus zwei Unterintervallen), zu erhalten. Die BCH-Formel muss noch auf den Spezialfall der Lie-Gruppe SO(3) spezialisiert werden. Zudem muss eine geeignete Auswertung / Berechnung der Reihenentwicklung erfolgen. Nach der Spezialisierung auf SO(3) wird eine geeignete Auswertung des dann immer noch als Reihe vorliegenden Ausdrucks durchgeführt.
• Die Spezialisierung auf die Lie-Gruppe SO(3) liefert zunächst nur die benötigte Relation für die Summation zweier Winkelelemente, für den zuvor beschriebenen Fall $${\theta }^k=2\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$$

  

  wobei y nach wie vor als Reihenentwicklung in der Gesamtordnung m wie in Gleichung (3.2) gegeben ist und θ^k das neu kumulierte Drehwinkelinkrement ist, das aus dem neu gemessenen relativen Drehwinkelinkrement Δθ^k und dem vorhergehend berechneten kumulierten Drehwinkelinkrement θ^k-1 berechnet wird.
• Die Auswertung der Reihenentwicklung gemäß Gleichung (3.2) der Funktion y in der Gleichung (3.3) erfolgt unter der Annahme, dass nur ihr zweites Argument klein ist und daher nur in diesem Argument eine Approximation durch das Taylorpolynom vorgenommen werden kann. Dieses ist erforderlich, denn bei iterativer Anwendung der Relation (3.3) auf insgesamt n Elemente gemäß der Vorschrift: $$\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,\text{ 0}}=\mathrm{0,}\frac{\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k}}{2}=\mathrm{\gamma }\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right),\mathrm{\Delta }{\mathrm{\theta }}_{out}^I=\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,n}$$

  

  beinhaltet das erste Argument von γ das (iterative) kumulative Drehwinkelinkrement $\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k-1}$

  

  aus den vorangegangenen k-1 Schritten. Dieses ist also für große k ≤ n nicht notweniger Weise klein. Daher kann eine Approximation der Taylorreihe im ersten Argument durch ihr Taylorpolynom zu einem inakzeptabel großen numerischen Fehler führen.
• Die BCH-Formel lässt sich nun für die Lie-Gruppe SO(3) wie folgt approximieren, um die Exponenten γ^(mθ) der Ordnungen m_β = 0, 1, 2 zu erhalten: $$\begin{array}{l}\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\cong {\displaystyle {\sum }_{m_{\beta }=0}^2{\gamma }^{\left(m_{\beta }\right)}}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right),\hfill \\ {\mathrm{\gamma }}^{\left(0\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{{\mathrm{\theta }}^{k-1}}{2},\hfill \\ {\mathrm{\gamma }}^{\left(1\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{\Delta {\mathrm{\theta }}^k}{2}+\frac{{\theta }^{k-1}}{2}x\frac{\Delta {\mathrm{\theta }}^k}{2}+B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\hfill \\ \left(\left(\frac{\Delta {\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\Delta {\mathrm{\theta }}^k}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}-{\left(\frac{\Delta {\mathrm{\theta }}^{k-1}}{2}\right)}^2\frac{\Delta {\mathrm{\theta }}^k}{2}\right),\hfill \\ {\gamma }^{\left(2\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{1}{2}{\left(\frac{2}{{\mathrm{\theta }}^{k-1}}\right)}^2\left(1-3B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(1-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^{2}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\right)\right)\hfill \\ \left({\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^2{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\right)\frac{{\mathrm{\theta }}^{k-1}}{2}+B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(1-{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^{k-1}}{2}\right)}^{1}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\right)\hfill \\ \left(\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}-{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\frac{{\mathrm{\theta }}^{k-1}}{2}\right)+B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\frac{{\mathrm{\theta }}^{k-1}}{2}\times \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2},\hfill \end{array}$$

  

  mit der transzendenten Funktion $$B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)\cong {\tilde{B}}_{Pade,N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\displaystyle {\sum }_{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m{\zeta }^{2m}}}{1+{\displaystyle {\sum }_{n=1}^{\left[\frac{N_0+1}{2}\right]}{b}_n{\zeta }^{2n}}}.$$

  
• Die einzige, in dem Gleichungssystem enthaltene und durch Approximation zu berechnende Funktion ist die transzendenten Funktion B(ζ).
• Die Herleitung der Approximation mindestens zweiter Ordnung m_β der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) wird später noch im Detail beschrieben.
• Es reicht aus, wenn lediglich im zweiten Argument ${\theta }_2=\frac{\mathrm{\Delta }{\mathrm{\theta }}^2}{2}$

  

  eine Approximation mindestens in den Ordnungen m_β = 0, 1, 2 vorgenommen wird. Die Ausdrücke sind hingegen im ersten Argument ${\theta }_1=\frac{{\mathrm{\theta }}^{k-1}}{2}$

  

  exakt.
• Es tritt nur eine einzige transzendente Funktion $\text{B}\left(\frac{{\mathrm{\theta }}_1}{2}\right)=\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)$

  

  auf, die mit Hilfe der Padé-Approximation effizient und genau approximiert werden kann.
• Die maximal erforderlichen Ordnungen der Taylorpolynome für $\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)$

  

  und den Ausdruck (2.6) hängen von der gewünschten numerischen Genauigkeit bei gegebenen maximal zu erwartenden Werten für die Winkelinkremente und deren Anzahl n ab, die in einem einzigen Winkelinkrement $\mathrm{\Delta }{\mathrm{\theta }}_{out}^I$

  

  zusammengefasst werden sollen.
• Bei Verwendung der Gleichungen (3.5) kann die erforderliche Ordnung bei gleicher numerischer Genauigkeit deutlich niedriger gewählt werden als bei Verwendung der Gleichungen (2.3), (2.5) und (2.6). Damit kann die Komplexität reduziert und der Rechenaufwand, d.h. Berechnungszeit und Hardwareaufwand, reduziert werden.
• Die iterative Berechnung des kumulierten Drehwinkels $\mathrm{\Delta }{\theta }_{out}^I$

  

  des Intervalls [t_I-1, t_I] mit Hilfe der Approximation mindestens zweiter Ordnung der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) kann auf Grundlage der Gleichungen (3.4) und (3.5a), (3.5) durch die Prozedur 1 wie folgt durchgeführt werden:
• Prozedur 1: $\mathrm{\Delta }{\mathrm{\theta }}_{out}^I\left(\mathrm{\Delta }{\mathrm{\theta }}^1,\dots ,\mathrm{\Delta }{\mathrm{\theta }}^n\right):$
  
• Setze $\mathrm{\Delta }{\theta }_{out}^{I,0}=0$
  
• Für k = 1,...,n iterativ:
  • Berechne $\mathrm{\Delta }{\theta }_{out}^{I,k}$
    
    aus dem (iterativen) kumulierten Drehwinkelinkrement $\Delta {\theta }_{out}^{I,k-1}$
    
    und dem gemessenen Drehwinkelinkrement Δθ^k gemäß Gleichung (3.4) mit der expliziten Form (3.5a), (3.5)
• Setze den resultierenden kumulierten Drehwinkel $\Delta {\theta }_{out}^I$
  
  gleich dem letzten (iterativen) kumulierten Drehwinkelinkrement $\Delta {\theta }_{out}^{I,n}.$
  
• 2 zeigt ein zugehöriges Flussdiagramm dieses Verfahrens bzw. der obigen Prozedur zur Bestimmung eines kumulierten Drehwinkelinkrementes aus mit einem Drehwinkelsensor für aufeinanderfolgende Zeitintervalle gemessenen Drehwinkelinkrementen.
• Berechnung des expliziten Ausdrucks für die BCH Formel für die Lie-Gruppe SO(3) Im folgenden wird ein Ausdruck für die BCH Formel für die Lie-Gruppe SO(3) hergeleitet, die Rotationen im dreidimensionalen Euklidischen Raum beschreibt. Strategie ist hierbei, die Entwicklung der Formel in Ordnungen nur ihres zweiten Argumentes $\frac{{\mathrm{\theta }}_2}{2}=\frac{\Delta {\mathrm{\theta }}_{out}^k}{2}$

  

  durchzuführen und die Abhängigkeit vom ersten Argument $\frac{{\mathrm{\theta }}_1}{2}=\frac{\Delta {\mathrm{\theta }}_{out}^{I,k-1}}{2}$

  

  $$

  

  in allen Ordnungen zu resummieren, so dass die Formel bezüglich ihres ersten Argumentes $\frac{{\mathrm{\theta }}_1}{2}$

  

  exakt ist. Das Ergebnis stellt also eine Entwicklung der Formel in Ordnungen ihres zweiten Argumentes $\frac{{\mathrm{\theta }}_2}{2}$

  

  bei gleichzeitiger Resummation aller Abhängigkeiten von ihrem ersten Argument $\frac{{\mathrm{\theta }}_1}{2}$

  

  dar.
• M. Müger. „Note on the theorem of Baker-Campbell-Hausdorff-Dynkin.", https://www.math.ru.nl/~mueger/PDF/BCHD.pdf gibt einen Überblick über einige Grundlagen, Fakten und Beweise zur BCH-Formel. Dort wird im Abschnitt 8 und insbesondere in Bemerkung 8.2 auf Seite 24 darauf hingewiesen, dass auf Grundlage einer speziellen Stufung der Algebra, basierend auf der Wortlänge im zweiten Argument anstelle der Gesamtwortlänge in beiden Argumenten, eine Entwicklung der BCH Formel im zweiten Argument erfolgen kann. Der entsprechende Ausdruck ist bis zu ersten Ordnung im zweiten Argument angegeben, ohne allerdings eine Spezialisierung auf die Lie-Gruppe SO(3) vorzunehmen und die die Bernoulli-Zahlen beinhaltende Serie explizit zu resummieren.
• 1. Berechnung der nullten und ersten Ordnung, m_β = 0,1
• Wird die in erster Ordnung einfachen Schritte der Spezialisierung und der Resummation unter Verwendung der Erzeugendenfunktion der Bernoullizahlen aus, welche durch $$\frac{z}{e^z-1}={\displaystyle \sum {}_{n=0}^{\infty }}\frac{B_n}{n!}{z}^n$$

  

  gegeben ist, ergibt sich für den in M. Müger „Note on the theorem of Baker-Campbell-Hausdorff-Dynkin.“ Als Reihe gegebenen Ausdruck der BCH-Formel die folgende explizite Darstellung in resummierter Form: $$\begin{array}{l}\gamma \left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\\ \frac{{\theta }_1}{2}+\frac{{\theta }_2}{2}+\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}+{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}{2}^{2n}}{\left(1n\right)!}}{\left(-{\left(\frac{{\theta }_1}{2}\right)}^2\right)}^{n-1}\frac{{\theta }_1}{2}\times \left(\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\right)+\mathcal{O}\left({\left(\frac{{\theta }_2}{2}\right)}^2\right)\\ \frac{{\theta }_1}{2}+\frac{{\theta }_2}{2}+\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}+B\left(\frac{{\theta }_1}{2}\right)\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right)+\mathcal{O}\left({\left(\frac{{\theta }_2}{2}\right)}^2\right)\end{array}$$

  
• Die Taylorreihe der Funktion $$\begin{array}{l}B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)=4{\displaystyle \sum {}_{k=1}^{\infty }\frac{B_{2k}}{\left(2k\right)!}{\left(2i\zeta \right)}^{2\left(k-1\right)}}=\frac{1}{3}+\frac{{\zeta }^2}{45}+\frac{2{\zeta }^4}{945}+\frac{{\zeta }^6}{4725}+\\ \frac{2\left({\zeta }^8\right)}{93555}+\frac{1382{\zeta }^{10}}{638512875}+O\left({\zeta }^{12}\right)\end{array}$$

  

  enthält nur gerade Ordnungen von ζ, so dass die Berechnung der Gleichung (A.2) keine Berechnung der in der Norm θ₁ des Vektors θ₁ enthaltenen Quadratwurzel erfordert. Dieses ist bei numerischen Berechnungen ein Vorteil.
• 2. Grundlagen zum Formalismus und Zusammenhang zwischen der Differentialgleichung der Winkelinkremente und der BCH Formel
• Die zur Herleitung der zweiten Ordnung erforderlichen Grundlagen und Schritte werden im Folgenden dargestellt.
• Ausgangspunkt bildet die Exponentialabbildung, die auch für Argumente α definiert ist, die Elemente einer Lie-Algebra, also des Tangentialraums einer Lie-Gruppe am Element der Identität sind, und die durch Matrizen dargestellt werden können. Die Exponentialabbildung ist über ihre Taylorreihe definiert $$e^{\alpha }={\displaystyle \sum {}_{n=0}^{\infty }}\frac{{\alpha }^n}{n!},$$

  

  wobei α ein Element einer Darstellung einer Lie-Algebra ist. Dieses Element α kann z.B durch die reinen Quaternionen im Fall der Rotationen gegeben sein. Potenzen in der obigen Gleichung sind so zu verstehen, dass sie mit dem entsprechenden Produkt durchzuführen sind, was für die Darstellung definiert ist, also z.B. dem Matrixprodukt im Fall einer Matrixdarstellung oder dem Quaternionenprodukt im Fall einer Darstellung durch Quaternionen.
• Da verschiedene Elemente α, β einer Lie-Algebra in der Regel nicht miteinander vertauschen, also α β ≠ β α gilt, ist die Ableitung der Exponentialabbildung nicht durch den bekannten Ausdruck gegeben. Vielmehr liefert die Anwendung der Ableitung auf die Reihendarstellung gemäß Gleichung (A.4) den folgenden Ausdruck:

  

  wobei der n-fache Kommutator folgendermaßen definiert ist: $$\begin{array}{lllll}{\left[\alpha ,\beta \right]}_{\left(n\right)}=\left[\alpha ,{\left[\alpha ,\beta \right]}_{\left(n-1\right)}\right],\hfill & \hfill & {\left[\alpha ,\beta \right]}_{\left(0\right)}=\beta ,\hfill & \hfill & {\left[\alpha ,\beta \right]}_{\left(1\right)}=\left[\alpha ,\beta \right]=\alpha \beta -\beta \alpha \hfill \end{array}$$

  

  und die adjungierten Wirkung des Elements α auf ein anderes Element β repräsentiert.
• Aus der Ableitung gemäß Gleichung (A.5) folgt nun die Relation

  

  
• Diese lässt sich invertieren. Dazu wird die Erzeugendenfunktion der Bernoulli-Zahlen verwendet, welche in Gleichung (A.1) gegeben ist. Das Inverse von der Relation (A.7) ist dann gegeben durch

  
• Es wird nun eine für numerische Berechnungen in der Praxis bedeutsame Eigenschaft erkennbar, nämlich, dass die Reihe in der Inverse gemäß (A.8) nur gerade Potenzen von α enthält. Deshalb hängt der dann auf die Lie-Gruppe SO(3) spezialisierte Ausdruck nur von geraden Potenzen der Norm von α ab und beinhaltet somit keine Quadratwurzeln. Zunächst wird die Relation für die Erzeugendenfunktion auf ihre geraden und ungeraden Potenzen projiziert, so dass die Gleichungen erhalten werden: $$\begin{array}{l}\frac{z}{e^z-1}-\frac{-z}{e^{-z}-1}=-z={\displaystyle \sum _{n=0}^{\infty }\frac{B_n}{n!}}\left(z^n-{\left(-z\right)}^n\right)=2B_{1}z+2{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n+1}}{\left(2n+1\right)!}{z}^{2n+1}},\hfill \\ \frac{z}{e^z-1}+\frac{-z}{e^{-z}-1}=z\frac{e^{\frac{z}{2}}+e^{-\frac{z}{2}}}{e^{\frac{z}{2}}-e^{-\frac{z}{2}}}=z\text{ coth }\frac{z}{2}={\displaystyle \sum _{n=0}^{\infty }\frac{B_n}{n!}}\left(z^n+{\left(-z\right)}^n\right)=2{\displaystyle \sum _{n=0}^{\infty }\frac{B_{2n}}{\left(2n\right)!}{z}^{2n}.}\hfill \end{array}$$

  
• Daraus folgt unmittelbar, dass $B_1=-\frac{1}{2},B_{2n+1}=0$

  

  für n ≥ 1 ist. Es gilt also: $$\frac{z}{e^z-1}=\frac{1}{2}z+{\displaystyle \sum _{n=0}^{\infty }\frac{B_{2n}}{\left(2n\right)!}{z}^{2n}}.$$

  
• Dieses Ergebnis wird verwendet, um die Inverse gemäß Gleichung (A.8) in folgende Form zu bringen:

  
• Dieser Ausdruck (A.11) wird nun spezialisiert auf die Lie-Gruppe SO(3), bzw ihrer doppelten Überlagerung, die Lie-Gruppe SU(2), dargestellt durch die Einheitsquaternionen. Für zwei Elemente θ, ω aus ihrer Lie-Algebra

  

  (2) reduzieren sich die multiplen Kommutatoren wie folgt: $$\begin{array}{cc}{\left[\frac{\theta }{2},\frac{\omega }{2}\right]}_{\left(n+2\right)}=-{\theta }^2{\left[\frac{\theta }{2},\frac{\omega }{2}\right]}_{\left(n\right)},& f\ddot{u}r\end{array}n\ge 1.$$

  
• Durch wiederholte Anwendung dieser Relation ergibt sich für die Gleichung (A.11) dann der Ausdruck

  
• Nach Identifikation des einfachen Kommutators mit dem Kreuzprodukt und Reduktion des doppelten Kommutators / Kreuzproduktes $$\begin{array}{cc}\left[\frac{\theta }{2},\frac{\omega }{2}\right]=2\frac{\theta }{2}\times \frac{\omega }{2},& \left[\frac{\theta }{2},\left[\frac{\theta }{2},\frac{\omega }{2}\right]\right]=4\frac{\theta }{2}\times \left(\frac{\theta }{2}\times \frac{\omega }{2}\right)=4\left(\left(\frac{\theta }{2}\cdot \frac{\omega }{2}\right)\frac{\theta }{2}-\frac{{\theta }^2}{2}\frac{\omega }{2}\right)\end{array}$$

  

  ist das Ergebnis (A.13) in der Tat die Differentialgleichung für das Winkelinkrement θ in Abhängigkeit von der Winkelgeschwindigkeit ω, die in der Gleichung (2.1) gegeben ist.
• Das Ergebnis der Gleichung (A.8) und im Spezialfall der Lie-Gruppe SO(3) dessen expliziter Ausdruck durch (A.13) sind eng verknüpft mit der Baker-Campbell-Hausdorff Formel, wie im Folgenden gezeigt wird. Die Baker-Campbell-Hausdorff Formel liefert das Lie-algebra-wertige Element γ als Funktion zweier Lie-algebrawertiger Elemente α, β, als Lösung folgender Relation $$e^{\gamma \left(\alpha ,\beta \right)}=e^{\alpha }{e}^{\beta }.$$

  
• Um diesen Zusammenhang zu zeigen, wird der Exponent y in Ordnungen von β entwickelt, wie dieses z.B. in M. Müger „Note on the theorem of Baker-Campbell-Hausdorff-Dynkin.“ Dargestellt ist. Hierzu wird einen Parameter t ∈ R eingeführt und die Entwicklung geschrieben als $$\begin{array}{cc}\gamma \left(\alpha ,\beta \right)={\displaystyle \sum _{m_{\beta }=0}^{\infty }{\gamma }^{\left(n\right)}\left(\alpha ,\beta \right),}& {\gamma }^{\left(m_{\beta }\right)}=\frac{1}{m_{\beta }!}{\frac{{\text{d}}^{m_{\beta }}\gamma \left(\alpha ,t\beta \right)}{\text{d}{t}^{m_{\beta }}}|}_{t=0}\end{array},$$

  

  wobei klarer Weise γ⁽⁰⁾= γ(α, 0) = α ist.
• Das Ergebnis für γ⁽¹⁾ in der Entwicklung der Gleichung (A.16) kann direkt aus der Gleichung (A.11) erhalten werden, wenn α durch γ(α, t β) ersetzt wird. Dieses führt dann auf die folgende Differentialgleichung für y:

  
• Durch Einsetzen von t = 0 in die obige Gleichung wird wegen der Gleichung (A.16) ein Ergebnis für γ⁽¹⁾ erhalten, das wie folgt lautet:

  
• Ein Vergleich dieses Ergebnisses mit der Gleichung (A.11) zeigt, dass die Differentialgleichung für das Lie-algebra-Element α bei gegebenem Lie-algebra-Element β = e^-α ∂_t e^α identisch mit der Gleichung für den in β linearen Beitrag in der BCH Reihe gemäß Gleichung (A.16) ist. Hinzu kommt, dass durch Einsetzen von t = Δt die Approximation von γ durch $\gamma \cong \alpha +{\gamma }^{\left(1\right)}\left(\alpha ,\mathrm{\Delta }t\beta \right)+\mathcal{O}{\left(\mathrm{\Delta }t\beta \right)}^2$

  

  zu einer Lösung der Differentialgleichung (A.11) zur linearen Ordnung in Δt β wird.
• 3 Berechnung der zweiten Ordnung, m_β = 2
• Im Folgenden wird der Term γ⁽²⁾ der BCH Reihe (A.16) hergeleitet. Dieser beinhaltet die Terme 2. Ordnung in β.
• Zunächst müssen hierfür die Ableitung der Relation (A.11) berechnet werden. Dies führt zu:

  
• Eine Gleichung für γ⁽²⁾ ist analog zum Vorgehen im Fall für γ⁽¹⁾ zu erhalten. Hierzu kann α durch γ(α, t β) ersetzt und außerdem β^̇ = 0,t = 0 eingesetzt werden. Damit ergibt sich

  
• Um ein explizites Ergebnis für den Spezialfall der Lie-Gruppe SO(3) zu erhalten, wird die Reihe in (A.20) zu einem geschlossenen Ausdruck summiert. Es kann analog wie zur Herleitung von (A.13) vorgegangen werden, wobei die Gleichung (A.12) verwendet wird, um die multiplen Kommutatoren zu reduzieren. Außerdem wird direkt $\mathrm{\alpha }=\frac{{\mathrm{\theta }}_1}{2},\mathrm{\beta }=\frac{{\mathrm{\theta }}_2}{2}$

  

  gesetzt, um das Ergebnis in den im Haupttext verwendeten Variablen zu erhalten. Um die Reduktionsformel (A.12) nutzen zu können, kann zunächst festgestellt werden, dass sich eine allgemeine Folge mit Argument f_k wie folgt in Folgen der nur geraden und nur ungeraden Anteile aufteilen lässt: $${\displaystyle \sum _{k=1}^n{f}_k={\displaystyle \sum _{k=1}^{\left[\frac{n}{2}\right]}{f}_{2k}}}+{\displaystyle \sum _{k=0}^{\left[\frac{n-1}{2}\right]}{f}_{2k+1}}.$$

  
• Mit dieser Identität lässt sich die innere Summe in (A.20) wie folgt aufteilen: $$\begin{array}{l}{\displaystyle \sum _{k=1}^{2n}{\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(k-1\right)}\right]\right]}_{\left(2n-k\right)}}\\ ={([\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]]}_{\left(2n-1\right)}+[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2n-1\right)}])\vartheta \left(n-1\right)\\ +{{\displaystyle \sum _{k=1}^{n-1}([\frac{{\theta }_1}{2},[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2k-1\right)}]]}}_{\left(2n-k\right)}+{[\frac{{\theta }_1}{2},[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2k\right)}]]}_{\left(2n-2k-1\right)}\\ ={\left(i{\theta }_1\right)}^{2\left(n-1\right)}(\left({\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]}_{\left(2n-1\right)}+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right)\vartheta \left(n-1\right)\\ +\left(n-1\right){\left(i{\theta }_1\right)}^{2\left(n-2\right)}{([\frac{{\theta }_1}{2},[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]]]}_{\left(2\right)}+[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right]}_{\left(2\right)}]]),\end{array}$$

  

  wobei ϑ die Heavyside Stufenfunktion ist. Einsetzen dieses Ergebnisses in die Reihe in (A.20) liefert die Taylorreihe: $$\begin{array}{l}{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}}{{\displaystyle \sum _{k=1}^{2n}\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(k-1\right)}\right]\right]}}_{\left(2n-k\right)}\\ ={\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}{\left(i{\theta }_1\right)}^{2\left(n-1\right)}\left(\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right)}\\ {\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}\left(n-1\right){\left(i{\theta }_1\right)}^{2\left(n-1\right)}\left({\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right]}_{\left(2\right)}+\left[\frac{{\theta }_1}{2},[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2\right)}\right]]\right)}.\end{array}$$

  
• Ein geschlossener Ausdruck für die erste der beiden Reihen ist bereits aus (A.9) bekannt. Ein geschlossener Ausdruck für die zweite Reihe kann durch Anwenden eines geeigneten Ableitungsoperators auf das bekannte Ergebnis erhalten werden: $$\begin{array}{l}{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}}{z}^{\left(n-1\right)}=\frac{\frac{z}{2}\text{coth}\frac{z}{2}-1}{z^2},\\ {\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}}\left(n-1\right)z^{2\left(n-2\right)}=\frac{1}{2z}\frac{\text{d}}{\text{d}z}{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}}{z}^{2\left(n-1\right)}\\ =\frac{2{\text{sinh}}^2\frac{z}{2}-\frac{z}{2}\text{sinh}\frac{z}{2}\text{cosh}\frac{z}{2}-{\left(\frac{z}{2}\right)}^2}{2z^4{\text{sinh}}^2\frac{z}{2}}\\ =\frac{1}{8z^2}\left({\left(\frac{2}{z}\right)}^2{\left(\frac{z}{2}\text{coth}\frac{z}{2}-1\right)}^2+3{\left(\frac{2}{z}\right)}^2\left(\frac{z}{2}\text{coth}\frac{z}{2}-1\right)-1\right).\end{array}$$

  
• Analytische Fortsetzung der Ausdrücke durch Identifikation z = i θ₁ führt dann für den Term der BCH Reihe, der quadratisch in θ₂ ist, auf das Ergebnis $$\begin{array}{l}{\gamma }^{\left(2\right)}=\frac{1}{4}\left[{\gamma }^{\left(1\right)},\beta \right]+\frac{1}{8}B\left(\frac{{\theta }_1}{2}\right)\left(\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right)\\ -\frac{1}{16{\theta }_1^2}\left(1-3B\left(\frac{{\theta }_1}{2}\right)+{\left(\frac{{\theta }_1}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)\right)\left({\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right]}_{\left(2\right)}+\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2\right)}\right]\right]\right),\end{array}$$

  

  das nur von einer einzigen transzendenten Funktion, B(ζ), abhängt, die bereits im Ausdruck für γ⁽¹⁾ erscheint und in (A.3) gegeben ist. Dieser Umstand ist für numerische Berechnungen von Vorteil, da nur eine einzige transzendente Funktion approximiert und berechnet werden muss.
• 4. Expliziter Ausdruck für die BCH Formel bis zur zweiten Ordnung, m_β = 0, 1, 2
• Das explizite Ergebnis der BCH Formel für SO(3) in der Entwicklung (A.16), d.h. bis einschließlich der Terme quadratischer Ordnung im zweiten Argument $\frac{{\theta }_2}{2}$

  

  wird im Folgenden zusammengefasst und vereinfacht. Der Term γ⁽⁰⁾ ist einfach durch das erste Argument $\frac{{\theta }_1}{2}$

  

  der Funktion gegeben. Der Term γ⁽¹⁾ ist durch (A.13) gegeben, wenn darin θ _→ θ₁ und ω _→ θ₂ ersetzt oder wahlweise direkt der Ausdruck (A.2) verwendet wird. Der Term γ⁽²⁾ ist in (A.25) dargestellt. Damit ergibt sich für die BCH Formel $$\begin{array}{l}\frac{\theta }{2}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\gamma \left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)={\displaystyle \sum _{m_{\beta }=0}^2{\gamma }^{\left(m_{\beta }\right)}}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)+\mathcal{O}\left({\left(\frac{{\theta }_2}{2}\right)}^3\right),\hfill \\ {\gamma }^{\left(0\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{{\theta }_1}{2},\hfill \\ {\gamma }^{\left(1\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{{\theta }_2}{2}+\frac{1}{2}\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]+\frac{1}{4}B\left(\frac{{\theta }_1}{2}\right)\left[\frac{{\theta }_1}{2},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right],\hfill \\ {\gamma }^{\left(2\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{1}{4}\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]+\frac{1}{8}B\left(\frac{{\theta }_1}{2}\right)\left(\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right)\hfill \\ -\frac{1}{16{\theta }_1^2}\left(1-3B\left(\frac{{\theta }_1}{2}\right)+{\left(\frac{{\theta }_1}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)\right)\left({\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right]}_{\left(2\right)}+\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2\right)}\right]\right]\right).\hfill \end{array}$$

  
• Nun werden die folgenden Identitäten verwendet: $$\begin{array}{c}\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]=2\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2},\\ \left[\frac{{\theta }_1}{2},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]=4\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right),\\ \left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]=2{\gamma }^{\left(1\right)}\times {\theta }_2,\\ \left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]=4\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right){\gamma }^{\left(1\right)}+\left({\gamma }^{\left(1\right)}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-2\left({\gamma }^{\left(1\right)}\cdot \frac{{\theta }_1}{2}\right)\frac{{\theta }_2}{2}\right),\\ {\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right]}_{\left(2\right)}+\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2\right)}\right]\right]=-8\left({\gamma }^{\left(1\right)},\frac{{\theta }_1}{2}\right)\left[\frac{{\theta }_1}{2},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\end{array}$$

  

  $$=32\left({\gamma }^{\left(1\right)}\cdot \frac{{\theta }_1}{2}\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right),$$

  

  und außerdem $$\begin{array}{c}{\gamma }^{\left(1\right)}\times {\theta }_2=\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}-{\left(\frac{{\theta }_2}{2}\right)}^2\frac{{\theta }_1}{2}+B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2},\\ {\gamma }^{\left(1\right)}\cdot \frac{{\theta }_1}{2}=\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2},\\ {\gamma }^{\left(1\right)}\cdot \frac{{\theta }_2}{2}={\left(\frac{{\theta }_2}{2}\right)}^2-B\left(\frac{{\theta }_1}{2}\right){\left(\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\right)}^2={\left(\frac{{\theta }_2}{2}\right)}^2+B\left(\frac{{\theta }_1}{2}\right)\left({\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2-{\left(\frac{{\theta }_1}{2}\right)}^2{\left(\frac{{\theta }_2}{2}\right)}^2\right).\end{array}$$

  

  um zunächst die in der Gleichung (A.26) auftretenden multiplen Kommutatoren zu berechnen und zu vereinfachen $$\begin{array}{l}\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\frac{{\theta }_2}{2}\right]\right]+\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]=4\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right){\gamma }^{\left(1\right)}+\left({\gamma }^{\left(1\right)}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-2\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}\right)\hfill \\ =4\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(\frac{{\theta }_2}{2}+\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}+B\left(\frac{{\theta }_1}{2}\right)(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right)\right)\hfill \\ +\left({\left(\frac{{\theta }_2}{2}\right)}^2+B\left(\frac{{\theta }_1}{2}\right)\left({\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2-{\left(\frac{{\theta }_1}{2}\right)}^2{\left(\frac{{\theta }_2}{2}\right)}^2\right)\right)\frac{{\theta }_1}{2}\hfill \\ -2\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2})\hfill \\ =4\left(2{\left({\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)+\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)(\frac{{\theta }_2}{2}\right)}^2\right)\frac{{\theta }_1}{2}\hfill \\ -\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(1+{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\frac{{\theta }_2}{2}\hfill \\ +\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}),\hfill \\ {\left[\frac{{\theta }_1}{2},\left[{\gamma }^{\left(1\right)},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\right]}_{\left(2\right)}+\left[\frac{{\theta }_1}{2},[{\gamma }^{\left(1\right)},{\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]}_{\left(2\right)}\right]]=-8\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)\left[\frac{{\theta }_1}{2},\left[\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right]\right]\hfill \\ =-32\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right).\hfill \end{array}$$

  
• Damit kann die Formel (A.26) schrittweise auf folgende Form gebracht werden: $$\begin{array}{l}\gamma \left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)={\displaystyle \sum _{m_{\beta }=0}^2{\gamma }^{\left(m_{\beta }\right)}}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)+\mathcal{O}\left({\left(\frac{{\theta }_2}{2}\right)}^3\right),\hfill \\ {\gamma }^{\left(0\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{{\theta }_1}{2},\hfill \\ {\gamma }^{\left(1\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{{\theta }_2}{2}+\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}+B\left(\frac{{\theta }_1}{2}\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right),\hfill \\ {\gamma }^{\left(2\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{1}{2}\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}-{\left(\frac{{\theta }_2}{2}\right)}^2\frac{{\theta }_1}{2}+B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\right)\hfill \end{array}$$

  

  $$\begin{array}{l}+\frac{1}{2}B\left(\frac{{\theta }_1}{2}\right)\left(2{(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)+\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right){\left(\frac{{\theta }_2}{2}\right)}^2\right)\frac{{\theta }_1}{2}-\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(1+{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\frac{{\theta }_2}{2}\\ +\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2})\\ +\frac{1}{2}{\left(\frac{2}{{\theta }_1}\right)}^2\left(1-3B\left(\frac{{\theta }_1}{2}\right)+{\left(\frac{{\theta }_1}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)\right)\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right)\\ =\frac{1}{2}{(-{\left(\frac{{\theta }_2}{2}\right)}^2+2{\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)+B\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\left(\frac{{\theta }_2}{2}\right)}^2\\ +{\left(\frac{2}{{\theta }_1}\right)}^2\left(1-3B\left(\frac{{\theta }_1}{2}\right)+{\left(\frac{{\theta }_1}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)\right){\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2)\frac{{\theta }_1}{2}\\ +\frac{1}{2}\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)-B\left(\frac{{\theta }_1}{2}\right)\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\left(1+{\left(\frac{{\theta }_1}{2}\right)}^{2}B(\frac{{\theta }_1}{2}\right)\right)-\left(1-3B\left(\frac{{\theta }_1}{2}\right)+{\left(\frac{{\theta }_1}{2}\right)}^2{B}^2\left(\frac{{\theta }_1}{2}\right)\right)\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right))\frac{{\theta }_2}{2}\\ +B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\\ =\frac{1}{2}\left({\left(\frac{2}{{\theta }_1}\right)}^2\left(1-3B\left(\frac{{\theta }_1}{2}\right)(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\right){\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2-\left(1-B\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\right){\left(\frac{{\theta }_2}{2}\right)}^2)\frac{{\theta }_1}{2}\\ +B\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}+B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\\ =\frac{1}{2}{\left(\frac{2}{{\theta }_1}\right)}^2\left(1-3B\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\right)\left({\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2-{\left(\frac{{\theta }_1}{2}\right)}^2{\left(\frac{{\theta }_2}{2}\right)}^2\right)\frac{{\theta }_1}{2}\\ +B\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^{2}B\left(\frac{{\theta }_1}{2}\right)\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}-{\left(\frac{{\theta }_2}{2}\right)}^2\frac{{\theta }_1}{2}\right)+B\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}.\end{array}$$

  
• 5. Berechnung höherer Ordnungen, m_β ≥ 3
• Das beschriebene Verfahren lässt sich direkt auch für die Berechnung höherer Ordnungen anwenden. Insbesondere die Ordnung m_β = 3 in θ₂, also der Term γ⁽³⁾ in der BCH Reihe (A.16) kann genutzt werden, um eine noch größere Genauigkeit zu erhalten.
• Die Berechnung der Exponenten γ^(mα) in der nullten, ersten und mindestens zweiten Ordnung m_α = 0, 1 und 2 erfordert eine approximative Berechnung lediglich einer einzigen transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right).$

  

  Diese ist mit dem Wert der halben quadrierten Norm des letzten kumulierten Winkelinkrementes, also $$

  

  $\zeta =\frac{{\theta }^{k-1}}{2}$

  

  mit ${\mathrm{\theta }}^{k-1}=\sqrt{{\mathrm{\theta }}^{k-1}\cdot {\mathrm{\theta }}^{k-1}}$

  

  auszuwerten. Hierbei ist die Funktion B(ζ) gerade, so dass nur die quadrierte Norm (θ^k-1)² benötigt wird und damit keine Berechnung der Quadratwurzel erforderlich ist.
• Für hinreichend kleine 0 ≤ ζ << 1, also kleine Norm (θ^k)² << 1 des kumulierten Winkelinkrementes θ^k kann die Funktion B(ζ) in der Gleichung (1.3) durch ihr Taylorpolynom in niedriger Ordnung approximiert werden. Dieser Fall relativ kleiner kumulierter Winkel tritt beispielsweise auf, wenn lediglich eine endliche Anzahl n relativer Winkelinkremente in Δθ^k kumuliert werden soll, deren Norm ebenfalls jeweils klein ist, also Δθ^k << 1 erfüllt. Da nach n Zeitschritten der kumulierte Winkel ausgegeben und dann dessen intern gespeicherter Wert wieder auf 0 gesetzt wird, handelt es ich bei den kumulierten Winkelinkrementen θ^k ebenfalls um relative Winkelinkremente, die auf ihren jeweiligen Startzeitpunkt der Kumulation bezogen sind. Die Genauigkeit der Approximation der Funktion B(e;) durch ihr Taylorpolynom wird somit dadurch bestimmt, wie groß n und die zu summierenden relativen Winkelinkremente Δθ^k, mit k = 1, ...,n sind. Beispielsweise ist für n relative Winkelinkremente Δθ^k konstanter Norm Δθ¹ = · · · = Δθⁿ = Δθ die Größe des Produktes nΔθ ausschlaggebend für die Genauigkeit der für die Auswertung von B (ζ) aus Gleichung (1.3) verwendeten Approximation.
• Im Folgenden wird die Renormierung der Winkelinkremente erläutert.
• Um den absoluten Drehwinkel und somit die Lage des Körpers im Raum zu erhalten, sind die ab einem Zeitpunkt, an dem die Lage bekannt war, alle relativen Winkelinkremente Δθ^k auzuaddieren. Dieses stellt also eine Kumulation beliebig vieler relativer Winkelinkremente Δθ^k ohne periodisches Zurücksetzen des Ergebnisses dar. Die Zahl n und daher auch das Argument ζ, mit welchem die Funktion B(ζ) in (1.3) auszuwerten ist, sind dann a priori nach oben nicht beschränk. Die Approximation von (1.3) durch das Taylorpolynom führt dann bei größer werdenden ζ zu einem immer größer werdenden und letztendlich unbeschränkt wachsenden und somit inakzeptablen Fehler.
• Dieses Problem der nach oben nicht beschränkten Norm der Winkelinkremente und somit eines unbeschränkten Argumentes e; kann wie im Folgenden dargestellt gelöst werden. Dieses ist eine der Voraussetzungen dafür, die periodische, transzendente Funktion (1.3) mit nach oben begrenztem Genauigkeitsverlust durch z.B. ihr Taylorpolynom zu approximieren.
• Es wird darauf hingewiesen, dass eine Rotation in zwei Dimensionen durch einen einzigen Drehwinkel ϕ bestimmt ist, der auf das Intervall [-π, π] oder bei Angabe des Vorzeichens ± der Drehrichtung auf das Intervall [0, π] eingeschränkt werden kann. Das bedeutet, dass der Drehwinkel ϕ und alle Drehwinkel ϕ_k = ϕ ± 2πk,k ∈ N dieselbe Rotation beschreiben. Sollte also ein Drehwinkel ϕ > π angegeben sein, so ist es immer möglich, diesen durch Subtraktion geeigneter ganzzahliger Vielfache von 2π auf seinen Repräsentanten ϕ' = ϕ -2πκ₀, der innerhalb des Intervalls [-π, π] liegt, zu reduzieren.
• Wird bei gegebenem ϕ > 0 die Größe ${\mathrm{\kappa }}_0=\text{sgn}\left(\mathrm{\varphi }\right)\left[\frac{\left|\mathrm{\varphi }\right|}{2\mathrm{\pi }}\right]$

  

  gewählt, wobei [ ] die Floor-Funktion ist, dann wird mit ϕ' = ϕ -2πκ₀ in der Tat der Repräsentant erhalten, der die Bedingung |ϕ'| ≤ π erfüllt.
• Bei einer Rotation in drei Dimensionen ist ein ähnliches Vorgehen möglich. Eine solche Rotation wird durch einen vektorwertigen Drehwinkel θ parametrisiert, dessen Norm (Länge) $\theta =\sqrt{\theta \cdot \theta }$

  

  den totalen (skalaren) Drehwinkel und dessen Richtung (Einheitsvektor) $e_{\theta }=\frac{\theta }{\theta }$

  

  die Drehachse festlegen. Da θ > 0 ist, legt e_θ auch gleich die Drehrichtung fest. Hierbei kann die Rechte-Faust- oder Korkenzieher-Regel angewendet werden, bei der der Daumen der rechten Hand in Richtung des Vektors zeigt und die Fingerspitzen der zur Faust geballten rechten Hand die Drehrichtung für Drehungen mit positivem, Winkel anzeigen, d.h. die Norm des vektorwertigen Drehwinkels.
• Analog zu den Rotationen in zwei Dimensionen kann daher in drei Dimensionen jede Rotation durch einen vektorwertigen Drehwinkel θ' = θ'e_θ, parametrisiert werden, dessen Norm 0 ≤ θ' < π erfüllt. Wenn ein allgemeiner Drehwinkel θ gegeben ist, dessen Norm θ nicht im Intervall [0, π] liegt, kann dieser durch einen Repräsentanten θ' mit minimaler Norm 0 ≤ θ' < π ersetzt werden, der aus θ wie folgt berechnet werden kann $$\mathrm{\theta }\text{'}=\left(\mathrm{\theta }-2\mathrm{\pi }{\mathrm{\kappa }}_0\right)e_{\theta }$$

  
• Hierbei ist k₀ ∈ ℕ. Norm und Richtungsvektor des Repräsentanten mit minimaler Norm sind dann durch $$\begin{array}{cc}\mathrm{\theta }\text{'}=\left|\mathrm{\theta }-2\mathrm{\pi }{\mathrm{\kappa }}_0\right|,& e_{\mathrm{\theta }\text{'}}=\text{sgn}\left(\mathrm{\theta }-2\mathrm{\pi }{\mathrm{\kappa }}_0\right)\end{array}{e}_{\theta }$$

  

  gegeben, d.h. die Bedingung 0 ≤ θ' ≤ π ist erfüllt.
• Im nächsten Schritt wird der Wert bestimmt, der für κ₀ zu wählen ist, um den Repräsentanten mit minimaler Norm zu erhalten. Ein allgemeiner Drehwinkel θ besitzt eine Norm 0, die entweder im Intervall 2kπ ≤ θ < (2k + 1)π oder im Intervall (2k + 1)π ≤ θ < 2(k + 1)kπ, k ∈ ℕ, liegt. In ersterem und letzterem Fall führt die Wahl der Konstanten κ₀ = k bzw. κ₀ = k + 1 auf folgende Ausdrücke: $$\begin{array}{c}\theta \text{'}=\{\begin{array}{ll}\left(\theta -2\pi k\right){\text{e}}_{\theta }\hfill & 2k\pi \le \theta <\left(2k+1\right)\pi \hfill \\ \left(\theta -2\pi \left(k+1\right)\right){\text{e}}_{\theta }\hfill & \left(2k+1\right)\pi \le \theta <2\left(k+1\right)\pi \hfill \end{array}=\\ \{\begin{array}{ll}\left|\theta -2\pi k\right|{\text{e}}_{\theta }\hfill & 2k\pi \le \theta <\left(2k+1\right)\pi \hfill \\ \left|\theta -2\pi \left(k+1\right)\right|\left(-{\text{e}}_{\theta }\right)\hfill & \left(2k+1\right)\pi \le \theta <2\left(k+1\right)\pi \hfill \end{array}\end{array}$$

  
• Diese können auch wie folgt dargestellt werden $$\theta \text{'}=\{\begin{array}{ll}\left(1-\frac{2k\pi }{\theta }\right)\theta \hfill & 2k\pi \le \theta <\left(2k+1\right)\pi \hfill \\ \left(1-\frac{2\left(k+1\right)\pi }{\theta }\right)\theta \hfill & \left(2k+1\right)\pi \le \theta <2\left(k+1\right)\pi \hfill \end{array}$$

  
• In der Tat erfüllt die Norm θ' nun die Bedingung 0 ≤ θ' < π. Die beiden Fälle lassen sich mit Hilfe der Floor-Funktion [ ] kompakt zu $$\begin{array}{cc}\theta \text{'}=\left(1-\frac{2\left[\frac{\kappa +1}{2}\right]\pi }{\theta }\right)\theta ,& \kappa \pi \le \theta <\left(\kappa +1\right)\pi \end{array}$$

  

  zusammenfassen.
• Ein Vergleich der Ausdrücke aus den Gleichungen (3.5) und (3.1) ergibt, dass für einen Drehwinkel θ, dessen Norm im Intervall κπ ≤ θ < (κ + 1)π liegt, die Konstante ${\mathrm{\kappa }}_0=\left[\frac{\mathrm{\kappa }+1}{2}\right]$

  

  gewählt werden muss, um dann mittels der Gleichung (3.1) den Repräsentanten θ' mit minimaler Norm zu erhalten. Das lässt sich unmittelbar verifizieren, wenn die Norm des Repräsentanten θ' unter Zuhilfename der Relation $k-2\left[\frac{k}{2}\right]=\frac{1-{\left(-1\right)}^k}{2},$

  

  welche für beliebige k ∈ ℕ gültig ist, berechnet wird. Damit wird die Bedingung 0 ≤ θ' ≤ π erfüllt.
• Bei dem oben detailliert beschriebenen Algorithmus für die Kumulation relativer Winkelinkremente liegt zu jedem Abtastzeitpunkt t_k das Quadrat der Norm des gespeicherten kumulierten Winkels oder Winkelinkrementes θ^k-1 vor, welches hier wieder abkürzend mit (θ^k-1)² bezeichnet wird. Auch, wenn die Norm des Winkelinkrements θ^k-1 minimal ist, also die Bedingung 0 ≤ θ^k-1 ≤ π erfüllt ist, ist dieses nicht notwendiger Weise mehr für das kumulierte Drehwinkelinkrement θ^k der Fall, das mit dem auf der Gleichung (1.2) beruhenden Algorithmus und dem aktuellen relativen Winkelinkrement Δθ^k berechnet wurde.
• Unter der Annahme, dass θ^k-1 bereits bei Bedarf durch seinen Repräsentanten mit minimaler Norm ersetzte Winkel ist, also die Bedingung 0 ≤ θ^k-1 ≤ π erfüllt ist, und dass das aktuelle relative Winkelinkrement Δθ^k die Norm κ_kπ ≤ Δθ^k < (κ_k + 1)π besitzt, also beispielsweise für 0 ≤ θ^k < π die Konstante κ_k = 0 ist, für π ≤ Δθ^k < 2π dann κ_k = 1 ist usw, ergibt sich für die quadrierte Norm bis einschließlich zur 2.Ordnung in Δθ^k. $$\begin{array}{l}{\left({\theta }^k\right)}^2={\left|2\gamma \left(\frac{{\theta }^{k-1}}{2},\frac{\Delta {\theta }^k}{2}\right)\right|}^2=\\ {\left({\theta }^{k-1}+\mathrm{\Delta }{\theta }^k\right)}^2-B\left(\frac{{\theta }^{k-1}}{2}\right){\left({\theta }^{k-1}\times \mathrm{\Delta }{\theta }^k\right)}^2\le {\left({\theta }^{k-1}+\mathrm{\Delta }{\theta }^k\right)}^2<{\left({\kappa }_k+2\right)}^2{\pi }^2\end{array}$$

  

  wobei eine Obergrenze mittels der Dreiecksungleichung für Vektoren ermittelt werden kann. Ein Vergleich dieser Abschätzung mit der Gleichung (3.4) zeigt, dass die Norm des aktuellen relativen Winkelinkrementes den Wert der Konstanten κ bestimmt, der für die Berechnung des Repräsentanten minimaler Norm erforderlich ist.
• Im Fall, dass alle relativen Winkelinkremente eine hinreichend kleine Norm haben, also 0 ≤ Δθ^k < π für alle k gilt, sind alle κ_k = 0 und damit ergeben sich in der Gleichung (3.5) die Fälle κ = 0, in denen das kumulierte Winkelinkrement bereits minimale Norm hat oder κ = 1. In der Praxis sollte die Einschränkung 0 ≤ Δθ^k < π die Regel sein, denn üblicherweise wird die Abtastzeit der Gyroskope angepasst an die zu erwartende Systemdynamik so gewählt, dass nur eine geringe Drehbewegung während eines einzigen Abtastzeitschrittes erfolgt, also Δθ^k << π ist.
• Im generischen Fall, wenn also individuelle κ_k so zu wählen sind, dass κ_kπ ≤ Δθ^k < (κ_k + 1)π gilt, gibt es z.B. folgende zwei Möglichkeiten für die Bestimmung von κ im (k + 1)ten Zeitschritt:
• — Zunächst erfolgt eine Bestimmung von $${\kappa }_S=\left[\frac{{\left({\theta }^k\right)}^2}{{\pi }^2}\right]$$
  
• Dann erfolgt eine Berechnung von κ als die natürliche Zahl, die kleiner oder gleich der Quadratwurzel von κ_s ist, also $$\kappa =\left[\sqrt{{\kappa }_S}\right]$$

  
• Mit der so erhaltenen Konstanten κ lässt sich folgende Größe berechnen $$\epsilon =\frac{{\left({\theta }^k\right)}^2}{{\kappa }^2{\pi }^2}-1$$

  
• Diese gibt den relativen, also auf κ²π² bezogenen Rest an, der nach Subtraktion der maximal möglichen ganzzahligen Vielfachen von π² von dem Quadrat der Norm (θ^k)² noch übrig bleibt. Aufgrund der Konstruktion gilt in jedem Fall 0 ≤ ε. Mit diesem Ergebnis lässt sich dann der Repräsentant minimaler Norm durch Umformung der Gleichung (3.5) wie folgt konstruieren $$\theta {\text{'}}^k=\left(1-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa \sqrt{1+\epsilon }}\right){\theta }^k$$

  

  - Alternativ, um die (wenn auch nur approximative) Berechnung der Quadratwurzel der Konstanten κ_s zu vermeiden, können ganzzahlige Werte κ' über das Intervall 1 ≤ κ' < (κ_k + 2) iteriert werden, entweder vorwärts startend von der unteren oder rückwärts startend von der oberen Intervallgrenze. Dann kann mit dem aktuellen Wert κ' jeweils $${\epsilon }_{\kappa \text{'}}=\frac{{\left({\theta }^k\right)}^2}{\kappa {\text{'}}^2{\pi }^2}-1$$

  

  berechnet und der Vorzeichenwechsel von der Größe ε_κ' detektiert werden. Der größte Wert von κ', für den ε_κ' ≥ 0 gilt, ist dann der Wert für κ und ε_κ ≥ 0. Mit diesen Werten kann dann über die Gleichung (3.10) der Repräsentant mit minimaler Norm berechnet werden.
• Der Ausdruck (3.10) für die Berechnung des Repräsentant mit minimaler Norm erfordert die Berechnung der Quadratwurzel $\sqrt{1+\epsilon }.$

  

  Im Allgemeinen ist κπ ≤ θ^k+1 < (κ + 1)π, woraus folgt, dass ε dann die Bedingung $0\le \mathrm{\epsilon }<\frac{2}{\mathrm{\kappa }}+\frac{1}{{\mathrm{\kappa }}^2}$

  

  erfüllt. Somit ist ε nicht notwendiger Weise klein, so dass die Quadratwurzel nicht einfach durch ihr Taylorpolynom in fester Ordnung approximiert werden kann. Jedoch wird, wie oben bereits erwähnt, in der Praxis überwiegend der Fall auftreten, dass die Norm der relativen Winkelinkremente Δθ^k die Bedingung Δθ^k << π erfüllt. Dieses ist ohnehin eine Bedingung, die für die Anwendbarkeit des auf dem Baker-Campbell-Hausdorff-Formel beruhenden Algorithmus gegeben sein muss und stellt damit keine zusätzliche Einschränkung dar. Dann ergibt sich aus der Gleichung (3.6), dass die Norm des kumulierten Winkels θ^k auf θ^k ≤ θ^k-1 + Δθ^k ≤ π + Δθ^k beschränkt ist. Eine Berechnung von ε z.B. mittels (3.11) mit κ = κ' = 1 ergibt dann, dass ε folgendermaßen beschränkt ist $$0\le \epsilon \le 2\frac{\mathrm{\Delta }{\theta }^k}{\pi }+\frac{{\left(\mathrm{\Delta }{\theta }^k\right)}^2}{{\pi }^2}<<1$$

  
• Wenn nun also ein Maximum Δ̃θ für die Norm der relativen Winkelinkremente für alle Zeitpunkte t_k z.B. aus den Spezifikationen der konkreten Anwendung abgeleitet werden kann, man also Δ̃θ so wählt, dass $\mathrm{\Delta }{\theta }^k\left(1+\frac{\mathrm{\Delta }{\theta }^k}{2\pi }\right)\le \tilde{\Delta }\theta$

  

  für alle k gilt, dann ergibt sich als Begrenzung für ε die Bedingung $0\le \epsilon \le 2\frac{\tilde{\Delta }\theta }{\pi }.$

  

  Durch Erhöhung der Abtastrate können die Δθ^k und somit auch Δ̃θ verkleinert werden. Eine Approximation von $\sqrt{1+\epsilon }$

  

  durch ihr Taylorpolynom wird möglich. Die maximale Ordnung der Taylorentwicklung, also der Grad des Taylorpolynoms zur Erreichung einer gewünschten Genauigkeit kann bestimmt werden. Aus der Gleichung (3.5) ergibt sich $$\begin{array}{l}\theta {\text{'}}^k=\left(1-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa \sqrt{1+\epsilon }}\right){\theta }^k=\left(1-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\displaystyle \sum _{n=1}^{\infty }\frac{1}{2^{2n}}\left(\begin{array}{c}2n\\ n\end{array}\right)}{\left(-\epsilon \right)}^n\right){\theta }^k=\\ \left(1-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\displaystyle \sum _{n=1}^{\infty }\frac{\left(2n-1\right)!!}{\left(2n\right)!!}{\left(-\epsilon \right)}^n}\right){\theta }^k\\ =\left(1-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\displaystyle \sum _{n=1}^{n_0}\frac{\left(2n-1\right)!!}{\left(2n\right)!!}{\left(-\epsilon \right)}^n}\right){\theta }^k+R_{n_0+1}\left(\epsilon \right){\theta }^k,\end{array}$$

  

  wobei das Restglied unter Verwendung der Stirlings Formel für das asymptotische Verhalten der Fakultäten wie folgt approximiert werden kann $$\begin{array}{l}{R}_{n_0+1}\left(\epsilon \right)=-\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\displaystyle \sum _{n_0+1}^{\infty }\frac{\left(2n-1\right)!!}{\left(2n\right)!!}{\left(-\epsilon \right)}^n}\simeq -\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\displaystyle \sum _{n_0+1}^{\infty }\frac{1}{\sqrt{\pi n}}{\left(-\epsilon \right)}^n}\\ \simeq -\frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }\frac{{\left(-\epsilon \right)}^{n_0+1}}{1+\epsilon }\end{array}$$

  
• Der durch Verwendung des Taylorpolynoms vom Grad n₀ verursachte Fehler in der Norm lässt sich nun wie folgt abschätzen $$\left|R_{n_0+1}\left(2\frac{\widehat{\delta }\theta }{\pi }\right)\right|\le \frac{2\left[\frac{\kappa +1}{2}\right]}{\kappa }{\left(2\frac{\widehat{\delta }\theta }{\pi }\right)}^{n_0+1}\le 2{\left(2\frac{\widehat{\delta }\theta }{\pi }\right)}^{n_0+1}$$

  
• Die Berechnung des Repräsentanten minimaler Norm kann mit folgenden Verfahrensabläufen durchgeführt werden:

  
• Daraus lässt sich leicht n₀ bei vorgegebener Genauigkeit und vorgegebenem Δ̃θ berechnen.
• Bei Überprüfung der Norm des kumulierten Winkelinkrements und ggf. Ersetzung durch den Repräsentanten minimaler Norm in jedem Zeitschritt ist im Fall relativer Winkelinkremente, die die Bedingung 0 ≤ Δθ^k < π erfüllen, der Fall κ = 1 der einzige, in dem die Ersetzung des kumulierten Winkelinkrementes durch seinen Repräsentanten minimaler Norm durchzuführen ist. Daher lässt sich der vereinfachte Algorithmus 2 formulieren.
• Im Fall, dass die Norm Δθ^k nicht beschränkt ist, kann der allgemeinere Algorithmus 3 eingesetzt werden.

  

  
• Eine Alternative kann z.B. auch die als Algorithmus 4 gegebene Variante sein. Modifikationen der Berechnung von θ'^k, z.B. in Form einer Ersetzung von r durch ε = r - 1 und anschließender Auswertung der Taylorreihe aus der Gleichung (3.13) bis zur gewünschten Ordnung n₀, können ebenfalls bei Bedarf leicht vorgenommen werden.

  

  
• Nachfolgend wird die Approximation der transzendenten Funktion B(ζ) aus der Gleichung (1.3) erläutert.
• Die Auswertung der in Gleichung (1.3) definierten transzendenten Funktion B(ζ) ist zu jedem Zeitpunkt t_k, zu dem ein neues relatives Winkelinkrement Δθ^k vorliegt, für die Berechnung des kumulierten Winkelinkrements θ^k erforderlich. Im k-ten Zeitschritt ist C; durch die halbe Norm des Ergebnisses für das kumulierte Winkelinkrement aus dem vorhergehenden Zeitschritt gegeben, also durch $\zeta =\frac{{\theta }^{k-1}}{2}.$

  

  Die Anwendung des oben zur Renormierung der Winkelinkremente beschriebenen Algorithmus garantiert dabei, dass die Norm θ^k-1 die Bedingung 0 ≤ θ^k-1 ≤ π erfüllt.
• In dem folgenden Abschnitt werden numerische Approximationen für B(C;) hergeleitet, die es erlauben, B(ζ) effizient und mit der gewünschten Präzision über den gesamten erforderlichen Wertebereich $0\le \zeta \le \frac{\pi }{2}$

  

  zu berechnen.
• Die in der Gleichung (1.3) definierte Funktion besitzt die folgende Taylorreihenentwicklung: $$B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta \text{ cot }\zeta \right)=4{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n}}{\left(2n\right)!}{\left(2\text{i}\zeta \right)}^{2\left(n-1\right)}}$$

  

  wobei B_2n die Bernoulli-Zahlen sind. Diese lassen sich mittels der Relation $$B_{2n}={\left(-1\right)}^{n+1}\frac{2}{{\left(2\pi \right)}^{2n}}\left(2n\right)!\zeta \left(2n\right)$$

  

  durch die Riemannsche ζ-Funktion, ausgewertet für gerade Argumente 2n, n ∈ ℕ ausdrücken. Die Riemannsche ζ-Funktion ist definiert über die Reihe $$\zeta \left(s\right)={\displaystyle \sum _{n=1}^{\infty }\frac{1}{n^s}}$$

  
• Einsetzen der Gleichung (3.20) in die Gleichung (3.19) ergibt dann für die Taylorreihe von B(e;) $$\begin{array}{l}B\left(\zeta \right)=2{\displaystyle \sum _{n=1}^{\infty }\frac{\zeta \left(2n\right)}{{\pi }^{2n}}{\zeta }^{2\left(n-1\right)}=2}{\displaystyle \sum _{n=1}^{n_0+1}\frac{\zeta \left(2n\right)}{{\pi }^{2n}}{\zeta }^{2\left(n-1\right)}+R_{n_0+1}\left(\zeta \right)},\hfill \\ R_{n_0+1}\left(\zeta \right)=2{\displaystyle \sum _{n=n_0+2}^{\infty }\frac{\zeta \left(2n\right)}{{\pi }^{2n}}{\zeta }^{2\left(n-1\right)}}\hfill \end{array}$$

  

  wobei R_n0+1 (ζ) das Restglied der Approximation von B(ζ) durch das Taylorpolynom vom Grad 2n₀ in ζ ist. Es ist wichtig, dass ζ(2n), d.h. die Riemannsche ζ-Funktion, nicht mit dem Argument ζ verwechselt wird.
• Die erfindungsgemäß vorgeschlagene Pade-Approximation ist bei annähernd gleichem Rechenaufwand für die Berechnung erheblich genauer, als die Verwendung des Taylorpolynoms zur Approximation. Die Eignung der Pade-Approximation beruht darauf, dass aus der Definition in Gleichung (3.21) der Riemannschen ζ-Funktion für ihre führende Asymptotik für großes Argument s gilt $$\zeta \left(s\right)\stackrel{s>>1}{\to }1+\frac{1}{2^s}$$

  
• Mit den Bedingungen, dass die ersten n₀ +3 Terme des Gesamtausdrucks mit denen der Taylorentwicklung übereinstimmen und unter Berücksichtigung der Tatsache, dass nach der Gleichung (3.23) Quotienten aus zwei ζ-Funktionen mit großem Argument in niedrigster Ordnung sich asymptotisch dem Wert 1 (Eins) annähern, kann der Ausdruck für das Restglied zu allen Ordnungen im Argument ζ approximativ als geometrische Reihe resummiert werden. Damit wird das Restglied R erhalten $$\begin{array}{l}{R}_{n_0+1}\left(\zeta \right)=2{\displaystyle \sum _{n=0}^{\infty }\frac{\zeta \left(2n+2n_0+4\right)}{{\pi }^{2\left(n+n_0+2\right)}}{\zeta }^{2\left(n+n_0+1\right)}}\\ =2\frac{\zeta \left(2n_0+4\right)}{{\pi }^{2\left(n_0+2\right)}}{\zeta }^{2\left(n_0+1\right)}{\displaystyle \sum _{n=0}^{\infty }\frac{\zeta \left(2n+2n_0+4\right)}{{\pi }^{2n}\zeta \left(2n_0+4\right)}{\zeta }^{2n}}\\ \simeq {\tilde{R}}_{n_0+1}\left(\zeta \right)=2\frac{\zeta \left(2n_0+4\right)}{{\pi }^{2\left(n_0+2\right)}}{\zeta }^{2\left(n_0+1\right)}{\displaystyle \sum _{n=0}^{\infty }{\left(\frac{\zeta \left(2n_0+6\right)}{{\pi }^2\zeta \left(2n_0+4\right)}{\zeta }^2\right)}^{{}^n}}\\ =2\frac{\zeta \left(2n_0+4\right)}{{\pi }^{2\left(n_0+2\right)}}\frac{{\zeta }^{2\left(n_0+1\right)}}{1-\frac{\zeta \left(2n_0+6\right)}{{\pi }^2\zeta \left(2n_0+4\right)}{\zeta }^2}\end{array}$$

  
• Wenn nun in die Gleichung (3.22) das Restglied durch seine in der Gleichung (3.24) definierte Approximation R̃_n0+1 (ζ) ersetzt wird, dann wird folgende Approximation B̃_n0 (ζ) für B(ζ) erhalten: $$B\left(\zeta \right)\simeq {\tilde{B}}_{n_0}\left(\zeta \right)=2{\displaystyle \sum _{n=0}^{n_0}\frac{\zeta \left(2n+2\right)}{{\pi }^{2n+2}}{\zeta }^{2n}}+2\frac{\zeta \left(2n_0+4\right)}{{\pi }^{2\left(n_0+2\right)}}\frac{{\zeta }^{2\left(n_0+1\right)}}{1-\frac{\zeta \left(2n_0+6\right)}{{\pi }^2\zeta \left(2n_0+4\right)}{\zeta }^2}$$

  
• Die obige Approximation besitzt nur rationale Koeffizienten, obwohl sie scheinbar transzendente Zahlen wie π und ζ(2n) enthält, denn es gilt ζ(2n) = r(n)π²ⁿ, wobei r(n) ∈ ℚ eine rationale Zahl ist, die von n abhängt. Die ersten in der Gleichung (3.25) auftretenden Koeffizienten lauten explizit

  
• In 4 sind die relativen Fehler bei Verwendung der Approximation durch die Gleichung (3.25) für verschiedene Werte von n₀ und bei Verwendung des Taylorpolynoms mit verschiedenem Grad n in ζ² über den gesamten relevanten Wertebereich $0\le \zeta \le \frac{\pi }{2}$

  

  dargestellt. Es ist wichtig anzumerken, dass für n = n₀ +2 die beiden Approximationen einen vergleichbaren Rechenaufwand besitzen. Für alle Kombinationen von n = n₀ +2 ist die Approximation durch die Gleichung (3.25) genauer als die Approximation durch das entsprechende Taylorpolynom. So erreicht die Approximation durch die (3.25) bereits bei n₀ = 10 über fast den gesamten Wertebereich die in der Software Mathematica voreingestellte Präzision „$MachinePrecision“ = 15.9546 (Zahl der Dezimalstellen für die Darstellung von Machine-Precision Zahlen). Dieses ist mit der Approximation durch Taylorpolynome erst für n = 24 und somit mit ungefähr doppeltem Rechenaufwand möglich.
• Mit Hilfe des in der Gleichung (3.23) angegebenen asymptotischen Verhaltens der Riemannschen ζ-Funktion kann der Fehler der Approximation durch Gleichung (3.25) abgeschätzt und begründet werden, warum diese eine gegenüber der Approximation durch das entsprechende Taylorpolynom deutlich genauere und damit effizientere Berechnung der Funktion B(ζ) erlaubt. Die in der Gleichung (3.24) durchgeführte Approximation, auf der die Gleichung (3.25) basiert, verursacht bei Berücksichtigung auch des nächstführenden Terms der asymptotischen Entwicklung (3.23) der Riemannschen ζ-Funktion folgenden Näherungsfehler $$\begin{array}{c}\delta R_{n_0+1}\left(\zeta \right)=R_{n_{0+1}}\left(\zeta \right)-{\tilde{R}}_{n_0+1}\left(\zeta \right)\\ =-2\frac{\zeta \left(2n_0+4\right)}{{\pi }^2\left(n_0+2\right)}{\zeta }^{2\left(n_0+1\right)}{\displaystyle \sum _{n=0}^{\infty }\left(\frac{\zeta \left(2n+2n_0+4\right)}{{\pi }^{2n}\zeta \left(2n_0+4\right)}-{\left(\frac{\zeta \left(2n_0+6\right)}{{\pi }^2\zeta \left(2n_0+4\right)}\right)}^n\right){\zeta }^{2n}}\\ =\frac{1}{2{\pi }^2}{\displaystyle \sum _{n=0}^{\infty }\left(1+\left(\frac{3}{4}n-1\right)4^n\right){\left(\frac{{\zeta }^2}{4{\pi }^2}\right)}^{n+n_0+1}}\\ =\frac{1}{2{\pi }^2}{\left(\frac{{\zeta }^2}{4{\pi }^2}\right)}^{n_0+1}\left({\displaystyle \sum _{n=0}^{\infty }{\left(\frac{{\zeta }^2}{4{\pi }^2}\right)}^n}+\frac{3}{4}{\displaystyle \sum _{n=0}^{\infty }n}{\left(\frac{{\zeta }^2}{{\pi }^2}\right)}^n-{\displaystyle \sum _{n=0}^{\infty }{\left(\frac{{\zeta }^2}{{\pi }^2}\right)}^n}\right)\\ =-\frac{9}{2{\pi }^2}\frac{{\left(\frac{{\zeta }^2}{4{\pi }^2}\right)}^{n_0+3}}{{\left(1-\frac{{\zeta }^2}{{\pi }^2}\right)}^2\left(1-\frac{{\zeta }^2}{4{\pi }^2}\right)}.\end{array}$$

  
• Dieser Fehler ist als Funktion von ζ monoton wachsend und damit maximal an der oberen Intervallgrenze, also bei $\zeta =\frac{\pi }{2}$

  

  Dort ergibt sich also für den relativen Fehler der Approximation der transzendenten Funktion B(ζ) durch die Gleichung (3.25) $$\left|\frac{B\left(\frac{\pi }{2}\right)-\tilde{B}\left(\frac{\pi }{2}\right)}{B\left(\frac{\pi }{2}\right)}\right|=\frac{\left|\delta R_{n_0+1}\left(\frac{\pi }{2}\right)\right|}{B\left(\frac{\pi }{2}\right)}=\frac{{\pi }^2}{4}\frac{9}{2{\pi }^2}{\left(\frac{1}{16}\right)}^{n_0+3}\frac{256}{9\cdot 15}=\frac{1}{120}{\left(\frac{1}{16}\right)}^{n_0+1}$$

  
• Dieser Fehler tendiert für wachsendes n₀ mit exponentieller Geschwindigkeit gegen Null. Dagegen wir der maximale Fehler der Approximation durch das Taylorpolynom vom Grad n = n₀ +2 in ζ² durch das Restglied aus der Gleichung (3.22) nach Ersetzung n₀→n = n₀ +2 bestimmt. Für dieses lässt sich für große n nach Ersetzung der Riemannschen ζ-Funktion durch ihre konstante Asymptotik aus der Gleichung (3.23) eine untere Grenze wie folgt angeben $$R_{n+1}\left(\zeta \right)=2{\displaystyle \sum _{n\text{'}=n+2}^{\infty }\frac{\zeta \left(2n\text{'}\right)}{{\pi }^{2n\text{'}}}}{\zeta }^{2\left(n\text{'}-1\right)}\ge \frac{2}{{\pi }^2}{{\displaystyle \sum _{n\text{'}=n+2}^{\infty }\left(\frac{{\zeta }^2}{{\pi }^2}\right)}}^{n\text{'}-1}=\frac{2}{{\pi }^2}\frac{{\left(\frac{{\zeta }^2}{{\pi }^2}\right)}^{n+2}}{1-\frac{{\zeta }^2}{{\pi }^2}}$$

  
• Auch wird das Maximum an der Obergrenze des Intervalls von ζ angenommen, also bei $\zeta =\frac{\pi }{2}.$

  

  Der relative Fehler der Approximation von B(ζ) durch das Taylorpolynom vom Grad n beträgt daher $$\left|\frac{B\left(\frac{\pi }{2}\right)-\tilde{B}\left(\frac{\pi }{2}\right)}{B\left(\frac{\pi }{2}\right)}\right|\ge \frac{{\pi }^2}{4}\frac{2}{{\pi }^2}\frac{{\left(\frac{1}{4}\right)}^{n+2}}{\frac{3}{4}}=\frac{2}{3}{\left(\frac{1}{4}\right)}^{n+2}$$

  
• Dieser Fehler tendiert für wachsendes n zwar ebenfalls gegen Null, es ergibt sich für den Grad n = n₀ +2 mit vergleichbarem Rechenaufwand ein Fehler von $$\left|\frac{B\left(\frac{\pi }{2}\right)-\tilde{B}\left(\frac{\pi }{2}\right)}{B\left(\frac{\pi }{2}\right)}\right|\ge \frac{2}{3}{\left(\frac{1}{4}\right)}^{n_0+4}$$

  
• Dieser ist erheblich größer als der Fehler gemäß Gleichung (3.28). Möchte man einen vergleichbaren Fehler mit der Taylorapproximation erreichen, so findet man durch Gleichsetzen der relativen Fehler die Relation $${\left(\frac{1}{4}\right)}^{n+1}=\frac{4}{5}{\left(\frac{1}{16}\right)}^{n_0+2},$$

  

  woraus folgt, das n > 2n₀ +3 erforderlich ist. Zur Erreichung derselben Präzision wie der der Approximation (3.25) mittels einer Taylorreihe muss also die Anzahl der zu berechnenden Terme und damit der Rechenaufwand etwa doppelt so groß sein.
• Die Approximation der transzendenten Funktion B(e;) wird nun weiter verbessert. Dabei wird ausgenutzt, dass die Approximation durch die Gleichung (3.25) ein Spezialfall der Pade-Approximation ist, die eine Funktion durch eine rationale Funktion von zwei Polynomen im Funktionsargument approximiert, deren Funktionswert und ersten Ableitungen am Entwicklungspunkt mit denen der zu approximierenden Funktion übereinstimmt. Im Fall der Approximation durch die Gleichung (3.25) ist die Funktion und das Zählerpolynom und das Nennerpolynom besitzen die Grade n₀ +1 bzw. 1 in ζ².
• Im Allgemeinen ist die Pade-Approximation von B(ζ) mit einem Zählerpolynom vom Grad M und einem Nennerpolynom vom Grad N in ζ² durch $$\frac{{\displaystyle \sum _{m=0}^M{a}_m{\zeta }^{2m}}}{1+{\displaystyle \sum _{n=1}^N{b}_n}{\zeta }^{2n}}$$

  

  gegeben und stimmt nach geeigneter Wahl der Koeffizienten a_m und b_n mit dem Funktionswert B(0) und den ersten n = 1,..., M + N Ableitungen ${\frac{d^{n}B\left(\zeta \right)}{d{\left({\zeta }^2\right)}^n}|}_{\zeta =0}$

  

  überein.
• Mit der Zahl N₀ = n₀ +2 = M + N zeigt sich, dass bei gegebenen 0 ≤ N₀ ≤ 10 die Pade-Approximationen mit M = N₀ - k und N = k und verschiedenen k = 1,..., N₀ verschieden genau die Funktion B(ζ) approximieren. Für die Wahl $k=N_0-\left[\frac{N_0}{2}\right]=\left[\frac{N_0+1}{2}\right]$

  

  $$

  

  ergibt sich im Intervall $0\le \zeta \le \frac{\pi }{2}$

  

  die Pade-Approximation mit der insgesamt größten Genauigkeit. Diese Spezialfälle bei jedem N₀ können daher mit der folgenden Gleichung (3.34) für die Pade-Approximation bezeichnet werden $${\tilde{B}}_{\text{Pade}\text{'},N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\displaystyle \sum _{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m{\zeta }^{2m}}}{1+{\displaystyle \sum _{n=1}^{\left[\frac{N_0+1}{2}\right]}{b}_n{\zeta }^{2n}}}$$

  
• Die Koeffizienten lauten wie folgt:

  
• Die oben angegebenen ganzzahligen Zähler und Nenner werden durch Linearkombinationen von Produkten von ^-Funktionswerten erzeugt, wobei jeder Term in der Summe gleiche gerade Transzendentalität besitzt. Von jeder dieser Kombinationen fester Transzendentalität wird dann diese durch Division durch eine entsprechende Potenz von π entfernt. Die sich ergebenden ganzzahligen Zähler und Nenner sollten daher auch in den Beziehungen zwischen ζ-Funktionen auftreten, wie sie z.B. aus den Shuffle- und Stuffle-Relationen folgen.
• In 5 sind dargestellten relativen Fehler $\left|\frac{B\left(\zeta \right)-\tilde{B}\left(\zeta \right)}{B\left(\zeta \right)}\right|$

  

  bei Verwendung der Approximation von B̃_N0 (ζ) durch die Gleichung (3.25) für verschiedene Werte von N₀ (Diagramm a)) und bei Verwendung der Pade-Approximation von B_{Padé, N0} (ζ) mit Gleichung (3.34) für verschiedene Werte von N₀ (Diagramm b)) über den gesamten relevanten Wertebereich $0\le \zeta \le \frac{\pi }{2}$

  

  dargestellt. Für N₀ = n₀ +2 besitzen die beiden Approximationen die gleiche Anzahl an Termen. Der Rechenaufwand zur Berechnung von B̃_{Pade, N0+2}(ζ) ist sogar im Vergleich zu dem für B̃_reg,N0 (ζ) noch etwas geringer, da wegen der gleichmäßigeren Verteilung des Gesamtgrades N₀ auf Zähler- und Nennerpolynom weniger Potenzen von ζ² berechnet werden müssen.
• Für alle Kombinationen von N₀ = n₀ +2 ist die Approximation der Funktion B(ζ) im Vergleich zu der Approximation durch die Gleichung (3.25) sogar noch genauer und dabei auch noch mit etwas geringerem Rechenaufwand berechenbar als die Approximation durch die Gleichung (3.25). So erreicht die Approximation der Funktion B(ζ) bereits bei N₀ = 9 bereits über fast den gesamten Wertebereich die in der Software Mathematica voreingestellte Präzision „$MachinePrecision“ = 15.9546 (Zahl der Dezimalstellen für die Darstellung von Machine-Precision Zahlen), während die Approximation (3.25) erst bei N₀ = 10 annäherend diese Genauigkeit erreicht.
• Es sei noch angemerkt, dass anders als der exakte Ausdruck in der Gleichung (1.3) keine der drei Approximationen durch die Gleichungen (3.22), (3.25) und (3.2) einen scheinbaren Pol bei e; = 0 besitzt. Die Approximationen können daher über den gesamten Wertebereich $0\le \zeta \le \frac{\pi }{2}$

  

  verwendet werden.
• Nachfolgend wird die Approximation der transzendenten Funktion B⁽²⁾(ζ) erläutert.
• In der Berechnung des kumulierten Winkelinkrementes gemäß Gleichung (1.2) in der in (A.30) gegebenen expliziten Form tritt neben der transzendenten Funktion β (ζ) auch die Kombination $$B^{\left(2\right)}\left(\zeta \right)=\frac{1}{2}\frac{1}{{\zeta }^2}\left(1-3B\left(\zeta \right)\left(1-{\zeta }^{2}B\left(\zeta \right)\right)\right)$$

  

  auf. Diese besitzt wie auch schon B(ζ) in der Gleichung (1.3) selbst scheinbar einen Pol bei e; = 0. Allerdings ist die Funktion regulär und bei e; = 0 ergibt sich $$B^{\left(2\right)}\left(0\right)=\frac{2}{5}.$$

  
• Analog zum Vorgehen für die im vorhergehenden Abschnitt erläuterte Konstruktion von Approximationen von B(ζ) können die verschiedenen dort vorgestellten Approximationen für B⁽²⁾(ζ) konstruiert werden. So ist beispielsweise die Taylorreihe von B⁽²⁾(ζ) durch $$B^{\left(2\right)}\left(\zeta \right)=-24{\displaystyle \sum _{n=1}^{\infty }\frac{B_{2n+2}}{\left(2n+1\right)!}}{\left(2\text{i}\zeta \right)}^{2\left(n-1\right)}$$

  

  gegeben. Wenn die Pade-Approximation von B(ζ) aus der Gleichung (3.34) direkt in der Gleichung (3.36) verwendet wird, liefert dies für diese Approximationen das beste Ergebnis. Dieses ist nicht überraschend, denn die Taylorreihe in der Gleichung (3.37) von B⁽²⁾(ζ) besitzt eben nicht die für die hohe Präzision der Pade-Approximation verantwortliche Eigenschaft der Taylorreihe von B(ζ) in der Gleichung (3.22), dass die Koeffizienten für große n exponentiell gegen eine Konstante konvergieren (oder selbst exponentielles Verhalten zeigen). Jedoch ist bei direkter Verwendung von der Gleichung (3.34) in der Approximation nach Gleichung (3.36) der Einfluss des scheinbar vorhandenen Pols noch sichtbar. Von Nachteil ist aber insbesondere, dass bei Verwendung dieser Strategie eine getrennte Behandlung des scheinbar vorhandenen Pols, also des Wertes e; = 0, notwendig ist.
• Der zuvor genannte Nachteil lässt sich vermeiden und gleichzeitig lässt sich das beste Approximationsergebnis von allen untersuchten Approximationen konstruieren. Dieses wird erreicht, wenn auf Basis der Pade-Approximation durch Gleichung (3.34) ein Ausdruck für B⁽²⁾(ζ) berechnet wird, der durch Entfernen des scheinbar vorhandenen Pols regularisiert wird.
• Hierfür kann die Formulierung der Pade-Approximation (3.34) durch Einführung zweier Polynome $$P\left(\zeta \right)={\displaystyle \sum _{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m}{\zeta }^{2\left(m-1\right)},Q\left(\zeta \right)={\displaystyle \sum _{n=1}^{\left[\frac{N_0+1}{2}\right]}{b}_n}{\zeta }^{2\left(n-1\right)}$$

  

  wie folgt umgeschrieben werden: $${\tilde{B}}_{\text{Pade}\text{'},N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\zeta }^{2}P\left(\zeta \right)}{1+{\zeta }^{2}Q\left(\zeta \right)}$$

  
• Die Verwendung dieser Darstellung der Pade-Approximation von B(e;) in der Definition von B⁽²⁾(ζ) in der Gleichung (3.36) erlaubt es, aus der erhaltenen Approximation den scheinbar vorhandenen Pol bei e; = 0 folgendermaßen zu entfernen $$\begin{array}{l}{B}^{\left(2\right)}\left(\zeta \right)=\frac{1}{2{\zeta }^2}\left(1-3B\left(\zeta \right)\left(1-{\zeta }^{2}B\left(\zeta \right)\right)\right)=\frac{1}{2}\left(\frac{1-3B\left(\zeta \right)}{{\zeta }^2}+3B{\left(\zeta \right)}^2\right)\\ \simeq {\tilde{B}}_{\text{reg},N_0}^{\left(2\right)}\left(\zeta \right)=\frac{1}{2}\left(\frac{Q\left(\zeta \right)-P\left(\zeta \right)}{1+{\zeta }^{2}Q\left(\zeta \right)}+3B_{\text{Pade}\text{'},N_0}{\left(\zeta \right)}^2\right).\end{array}$$

  
• Hier bezeichnet $B_{reg,N_0}^{\left(2\right)}\left(\zeta \right)$

  

  die sich ergebene Approximation, die auf der Pade-Approximation von B(ζ) basiert aber keine numerisch problematischen scheinbaren Pole mehr enthält und somit regulär ist.
• In 5 sind die relativen Fehler bei direkter Verwendung der Pade-Approximation (3.25) zur Berechnung von B⁽²⁾(ζ) und bei Verwendung der regularisierten Approximation $B_{reg,N_0}^{\left(2\right)}\left(\zeta \right)$

  

  aus der Gleichung (3.40) für verschiedene Werte von N₀ über den gesamten relevanten Wertebereich $0\le \zeta \le \frac{\pi }{2}$

  

  dargestellt.
• Das Diagramm a) zeigt die relativen Fehler $\left|\frac{B^{\left(2\right)}\left(\zeta \right)-{\tilde{B}}^{\left(2\right)}\left(\zeta \right)}{B^{\left(2\right)}\left(\zeta \right)}\right|$

  

  wobei B̃⁽²⁾(ζ) durch das Einsetzen der Approximation von B̃_Pade,_N0 (ζ) aus der Gleichung (3.34) in die Gleichung (3.36) gegeben ist. Diagramm b) zeigt den relativen Fehler, wenn B̃⁽²⁾(ζ) durch die regularisierte Approximation B⁽²⁾ _reg,N0 (ζ) aus der Gleichung (3.40) gegeben ist.
• Insbesondere für N₀ > 8 und kleine ζ ist der regularisierte Ausdruck B⁽²⁾ _reg,N0 (ζ) aus der Gleichung (3.40) erheblich genauer als die Approximation, die auf das direkte Einsetzen der Pade-Approximation von B(ζ) in die Definition der Gleichung (3.36) von B⁽²⁾(ζ) beruht.
• Die Ergebnisse der Approximationen können zur Berechnung des kumulierten Winkels aus relativen Winkelinkrementen genutzt werden, um damit die Lage des Objektes zu bestimmen.
• Hierzu kann die in der Zeile 4 des Algorithmus 1 angegebene Operation verfeinert werden durch
• 4: Berechne ${\theta }_{out}^{I,k}\text{ aus }{\theta }_{out}^{I,k-1}\text{ und }{\theta }_{out}^k$

  

  gemäß Gleichung (3.4), d.h. in Gleichung (1.2) mit der expliziten Form nach Gleichungen (3.5a), (3.5). Diese Operation wertet die folgende Gleichung auf Basis der mit Hilfe der erläuterten Pade-Approximationen angepassten Gleichung (3.5) aus: $$\begin{array}{l}\gamma \left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)\simeq {\displaystyle \sum _{m_{\beta }=0}^2{\gamma }^{\left(m\beta \right)}}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)\hfill \\ {\gamma }^{\left(0\right)}\left(\frac{{\theta }_1}{2}.\frac{{\theta }_2}{2}\right)=\frac{{\theta }_1}{2}\hfill \\ {\gamma }^{\left(1\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)=\frac{{\theta }_2}{2}+\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}+{\tilde{B}}_{\text{Padé},N_0}\left(\frac{{\theta }_1}{2}\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_1}{2}-{\left(\frac{{\theta }_1}{2}\right)}^2\frac{{\theta }_2}{2}\right)\hfill \\ {\gamma }^{\left(2\right)}\left(\frac{{\theta }_1}{2},\frac{{\theta }_2}{2}\right)={\tilde{B}}_{\text{reg},N_0}^{\left(2\right)}\left(\frac{{\theta }_1}{2}\right)\left({\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)}^2-{\left(\frac{{\theta }_1}{2}\right)}^2{\left(\frac{{\theta }_2}{2}\right)}^2\right)\frac{{\theta }_1}{2}\hfill \\ +{\tilde{B}}_{\text{Padé},N_0}\left(\frac{{\theta }_1}{2}\right)\left(1-{\left(\frac{{\theta }_1}{2}\right)}^2{\tilde{B}}_{\text{Padé},N_0}\left(\frac{{\theta }_1}{2}\right)\right)\left(\left(\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\right)\frac{{\theta }_2}{2}-{\left(\frac{{\theta }_2}{2}\right)}^2\frac{{\theta }_1}{2}\right)\hfill \\ +{\tilde{B}}_{\text{Padé},N_0}\left(\frac{{\theta }_1}{2}\right)\frac{{\theta }_1}{2}\cdot \frac{{\theta }_2}{2}\frac{{\theta }_1}{2}\times \frac{{\theta }_2}{2}\hfill \end{array}$$

  
• Diese beinhaltet nun die folgenden Schritte, wobei angenommen wird, dass θ^k-1 entweder direkt vor Aufruf des folgenden Algorithmus 5 in diesem Zeitschritt oder am Ende der Prozedur im letzten Zeitschritt mittels beispielsweise des Algorithmus 2 durch den Repräsentanten mit minimaler Norm ersetzt werden.

  
• Es wird darauf hingewiesen, dass bei der obigen abkürzenden Schreibweise ${\theta }^k=\Delta {\theta }_{out}^{I,k}\text{ und }\Delta {\theta }^k=\Delta {\theta }_{out}^{i-1+k}$

  

  $$

  

  zu identifizieren sind. Die Genauigkeit wird durch den
• Abbruch der Entwicklung der Gleichung (3.41) nach der zweiten Ordnung in $\frac{\Delta {\mathrm{\theta }}^k}{2}$

  

  limitiert, was den in der Praxis vor allem relevanten Fall Δθ^k << 1 jedoch abdeckt. Bei Bedarf kann eine Erweiterung der Gleichung (3.41) mit den oben beschriebenen Methoden ohne konzeptionelle Hindernisse vorgenommen werden.
• Im Folgenden wird die Berechnung trigonometrischer Funktionen und der Rotationsmatrix erläutert.
• Die im Sensorpaket der IMU vorhandenen Beschleunigungssensoren messen je nach Messprinzip entweder direkt die Beschleunigung oder sie messen relative Geschwindigkeitsinkremente, die in Analogie zur Relation (1.1) für die relativen Winkelinkremente durch $$\Delta v_S^i={\displaystyle {\int }_{t_{i-1}}^{t_{i-1}+\Delta t_i}\begin{array}{c}\text{S}\left(\text{t}\text{'}\right)\\ \text{S}\left(\text{t}\text{'}\right),\text{I}\end{array}}\text{a}\left(t\text{'}\right)\text{d }t\text{'}$$

  

  gegeben sind. Hier bezeichnet $\begin{array}{c}S\left(t\right)\\ S\left(t\right),I\end{array}a\left(t\right)$

  

  die Beschleunigung des Beschleunigungssensors bezüglich eines Inertialsystems I, ausgedrückt im Bezugssystem S des Sensors am selben Zeitpunkt t, an dem der Wert der Beschleunigung a(t) vorliegt. Wenn der Sensor also Geschwindigkeitsinkremente misst, so sind diese also aus einer Integration der Beschleunigung im mitbewegten Bezugssystem S entstanden. Insbesondere sind in der Gleichung (1.1) daher auch Drehbewegungen enthalten, die der Sensor während der Integrationszeit ausführt. Die Geschwindigkeitsinkremente gemäß Gleichung (1.1) können daher nicht direkt zur Bestimmung der Geschwindigkeit und Position des Fahrzeugs gegenüber eines Inertialsystems oder anderen Referenz-Bezugssystems verwendet werden. Benötigt werden stattdessen Geschwindigkeitsinkremente, die durch Integration der Beschleunigung in einem zumindest im Integrationsintervall feststehenden Koordinatensystem berechnet worden sind. Das entsprechende Integral ist von der Form $${}_{\text{S},\text{I}}{}^{{\text{S}}_{\text{i}-1}}\Delta v^i={\displaystyle {\int }_{t_{i-1}}^{t_{i-1}+\Delta t_i}\begin{array}{r}\hfill \text{S}\left({\text{t}}_{\text{i}-1}\right)\\ \hfill \text{S}\left(\text{t}\text{'}\right)\end{array}}\mathcal{A}\begin{array}{r}\hfill \text{S}\left(\text{t}\text{'}\right)\\ \hfill \text{S}\left(\text{t}\text{'}\right),\text{I}\end{array}\text{a}\left(t\text{'}\right)\text{d }t\text{'},$$

  

  wobei die Rotationsmatrix $\begin{array}{l}S\left(t_{i-1}\right)\hfill \\ S\left(t\right)\hfill \end{array}A$

  

  die im Sensorsystem S zum Zeitpunkt t gemessene Beschleunigung in ein am Startzeitpunkt t_i-1 der Integration dem Sensorsystem entsprechendes aber nicht mit diesem mitrotierendes Koordinatensystem transformiert, d.h. es ist $\begin{array}{l}S\left(t_{i-1}\right)\hfill \\ S\left(t\right),I\hfill \end{array}a\left(t\right)=\begin{array}{l}S\left(t_{i-1}\right)\hfill \\ S\left(t\right)\hfill \end{array}A\begin{array}{l}S\left(t_i\right)\hfill \\ S\left(t\right),I\hfill \end{array}a\left(t\text{'}\right)$

  

  und somit ist $\begin{array}{l}{S}_{i-1}\hfill \\ S,I\hfill \end{array}\Delta v^i$

  

  das auf das Sensorbezugssystem S_i-1 = S(t_i-1) bezogene und in diesem auch durch Integration erhaltene relative Geschwindigkeitsinkrement. Um nun eine Anzahl n dieser relativen Geschwindigkeitsinkremente aufaddieren zu können, sind die individuellen Geschwindigkeitsinkremente wiederum in ein gemeinsames Bezugssystem, z.B. S_l-1 zum Startzeitpunkt t_I-1 der Summation, zu transformieren. Die Summation lautet (ohne Berücksichtigung von Gravitation) also $$\begin{array}{r}\hfill {\text{S}}_{\text{I}-1}\\ \hfill \text{S},\text{I}\end{array}\Delta v^{I,n}={\displaystyle \sum _{k=1}^n\begin{array}{r}\hfill {\text{S}}_{\text{I}-1}\\ \hfill {\text{S}}_{\text{I}-1+\text{k}}\end{array}\mathcal{A}}\begin{array}{r}\hfill {\text{S}}_{\text{I}-1+\text{k}}\\ \hfill \text{S},\text{I}\end{array}\Delta v^k.$$

  
• Hierin ist $\begin{array}{l}{\text{S}}_{I-1}\hfill \\ {\text{S}}_{I-1+k}\hfill \end{array}A$

  

  die mit dem vom Startzeitpunkt t_I-1 aus kumulierten relativen Winkelinkrement $\Delta {\mathrm{\theta }}_{out}^{I\mathrm{,0}}$

  

  assoziierte Rotationsmatrix. Dieses motiviert, warum die mit dem jeweiligen kumulierten Winkelinkrement $\Delta {\mathrm{\theta }}_{out}^{I\mathrm{,0}}$

  

  assoziierte Rotationsmatrix ${}_{S_{I-1+k}}^{S_{I-1}}A$

  

  benötigt wird.
• Der Algorithmus ist nicht nur für die oben beschriebene Anwendung, der Verarbeitung der Sensordaten einer IMU und der effizienten Berechnung von Position und Lage eines Fahrzeugs relevant. Darüber hinaus sollte insbesondere die hier vorgestellte effiziente Berechnung trigonometrischer Funktionen vielfältige Anwendungsmöglichkeiten in verschiedenen Bereichen haben, z.B. in Computer Vision.
• Es wird im Folgenden gezeigt, dass das für das Update des kumulierten Winkelinkrementes θ^k bereits berechnete Ergebnis für den Funktionswert der Funktion B(ζ) benutzt werden kann, um effizient und insbesondere ohne Approximation irgendwelcher weiterer transzendenter Funktionen die (trigonometrischen) Funktionen $\frac{sin2\zeta }{2\zeta }$

  

  und cos2ζ zu erhalten.
• Mit diesen Ergebnissen lässt sich dann direkt die assoziierte Rotationsmatrix ${}_F{}^G\text{A}{}^T$

  

  (θ) der aktiven Rotation zwischen den Koordinatensystemen F und G berechnen, die mittels der Drehung aufeinander abgebildet werden. Die aktiven Rotation transformiert die Basisvektoren $e_i^F$

  

  des Koordinatensystems F in die Basisvektoren $e_i^G={}_F{}^G\text{A}{}^T\left(2\zeta \right)$

  

  $$

  

  des nach Anwendung der durch den Winkel θ beschriebenen Drehung erhaltenen Koordinatensystems G. Die Kenntnis dieser Matrix ist ein wesentlicher Baustein für die Berechnung der translatorischen Geschwindigkeit und letztendlich der Position aus den von den Beschleunigungssensoren gemessenen Beschleunigungen oder relativen Geschwindigkeitsinkrementen.
• Es wird $\zeta =\frac{{\mathrm{\theta }}_1}{2}=\frac{{\mathrm{\theta }}^{k-1}}{2}$

  

  identifiziert, wobei ${\mathrm{\theta }}_1=\sqrt{{\theta }_1\cdot {\theta }_1}$

  

  ist und θ₁ das erste Argument in der dem Update des kumulierten Winkelinkrementes zugrundeliegenden Gleichung (3.41) bezeichnet und ${\mathrm{\theta }}^{k-1}=\sqrt{{\theta }^{k-1}\cdot {\theta }^{k-1}}$

  

  ist und $\frac{{\mathrm{\theta }}^{k-1}}{2}$

  

  der Repräsentant minimaler Norm des kumulierten Winkelinkrementes aus dem letzten Zeitschritt ist. Allerdings ist der folgende Algorithmus nicht auf die konkrete Anwendung auf die kumulierten Winkelinkremente θ^k-1 beschränkt. Vielmehr könnte θ₁ irgendein Winkel und $\zeta =\frac{{\mathrm{\theta }}_1}{2}$

  

  sein halber Betrag sein.
• Zunächst werden Ausdrücke für (trigonometrische) Funktionen hergeleitet, die von der halben Norm des kumulierten Winkelinkrementes abhängen, also von ζ, wobei die Abhängigkeit mittels B(ζ) ausgedrückt wird.
• Aus der Definition von B(ζ) in der Gleichung (1.3) folgt $$\text{cot }\zeta =\frac{1-{\zeta }^2\text{B}\left(\zeta \right)}{\zeta }$$

  

  womit nach Einsetzen in die bekannten Relationen $$\text{sin }\zeta =\frac{1}{\sqrt{1+{\text{cot}}^2\zeta }},\text{ cos }\zeta =\frac{\text{cot }\zeta }{\sqrt{1+{\text{cot}}^2\zeta }}$$

  

  die für Argumente $0\le \zeta \le \frac{\pi }{2}$

  

  gültig sind, die folgenden Ausdrücke hergeleitet werden können $$\frac{\text{sin }\zeta }{\zeta }=\frac{1}{\sqrt{{\zeta }^2+{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}},\text{ cos }\zeta =\frac{1-{\zeta }^{2}B\left(\zeta \right)}{\sqrt{{\zeta }^2+{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}}$$

  
• Das Argument $\zeta =\frac{{\mathrm{\theta }}_1}{2}=\frac{{\mathrm{\theta }}^{k-1}}{2}$

  

  ist lediglich die halbe Norm des kumulierten Winkelinkrementes. Es werden jedoch die trigonometrische Funktionen der Norm selbst benötigt. Diese können aus den bekannten Relationen zwischen trigonometrischen Funktionen mit einfachem und doppeltem Winkel erhalten werden. Dies führt zu $$\begin{array}{l}\frac{\text{sin }2\zeta }{2\zeta }=\frac{\text{sin }\zeta \text{ cos}\zeta }{\zeta }=\frac{1-{\zeta }^{2}B\left(\zeta \right)}{{\zeta }^2+{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}\\ \text{cos }2\zeta =1-2{\text{ sin}}^2\zeta =-\frac{{\zeta }^2-{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}{{\zeta }^2+{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}\end{array}$$

  
• Ein Vorteil der Tatsache, dass die Funktion B(ζ) nur von der halben Norm des Drehwinkels abhängt, ist also, dass durch die dann notwendige Verdoppelung des Argumentes der trigonometrischen Funktionen die Quadratwurzeln eliminiert werden und sich die Ausdrücke in Gleichung (3.45) ergeben, die rationale Funktionen in ζ² und B(ζ) sind. Diese enthalten in Form der Kombinationen $\frac{sin2\zeta }{2\zeta }$

  

  und cos2ζ, die gerade Potenzreihen in ζ sind, dann nicht einmal mehr eine Abhängigkeit von der Quadratwurzel aus der Berechnung der Norm selbst.
• Durch Verwendung der Approximationen gemäß Gleichung (3.25) oder (3.34) können also die trigonometrischen Funktionen durch rein rationale Ausdrücke sehr effizient approximiert werden. Bei Verwendung der Gleichung (3.34) wird $$\begin{array}{l}\frac{{\text{sin}}_{N_0}2\zeta }{2\zeta }=\frac{1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)}{{\zeta }^2+{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2},\\ {\text{cos}}_{N_0}2\zeta =\frac{{\zeta }^2-{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2}{{\zeta }^2+{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2}\end{array}$$

  

  definiert.
• In der 6 sind die relativen Fehler der Approximationen von sin θ (Diagramm a)) und cos θ (Diagramm b)) mit der Pade-Approximation nach Gleichung (3.34) in den Approximationen gemäß Gleichung (3.46) in dem relevanten Wertebereich 0 ≤ θ ≤ π dargestellt. In der 7 sind die relativen Fehler der Approximationen von sin θ (Diagramm a)) und cos θ (Diagramm b)) bei der direkten Approximation durch Taylorpolynome vergleichbaren Grades dargestellt. 8 zeigt die relativen Fehler der Approximationen von $2\text{ sin}\left(\frac{\theta }{2}\right)\text{ cos}\left(\frac{\theta }{2}\right)$

  

  (Diagramm a)) und $1-2{\text{ sin}}^2\left(\frac{\theta }{2}\right)$

  

  (Diagramm b)) bei der direkten Approximation durch Taylorpolynome vergleichbaren Grades dargestellt.
• Es ist zu erkennen, dass der relative Fehler der Approximationen von sin θ wegen der Nullstelle bei θ = π ansteigt.
• Die Verwendung der Pade-Approximation durch Gleichung (3.34) in den in diesem Abschnitt hergeleiteten Approximationen gemäß Gleichung (3.46) liefert gegenüber einer direkten Approximation durch ein Taylorpolynom vergleichbaren Grades über den gesamten Wertebereich die genauere Approximation. Die Approximation der trigonometrischen Funktionen mit halbem Winkel durch Taylorpolynome und die sich daran anschließende Berechnung der Funktionen mit dem ganzen Winkel als Argument ist eine deutliche Verbesserung gegenüber der direkten Approximation durch die entsprechenden Taylorpolynome. Allerdings ist die hier entwickelte Methode bei vergleichbarem Rechenaufwand insbesondere für kleinere Winkel erheblich genauer. Ein weiterer Vorteil der Approximation durch die Gleichung (3.46) ist, dass die Bedingung sin² θ + cos² θ = 1 unabhängig von der gewählten Ordnung N₀ erfüllt wird.
• Die mit einem vektorwertigen Drehwinkel θ assoziierte Rotationsmatrix der aktiven Rotation lautet

  

  wobei L ein aus der Matrixdarstellung $$L_1=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -1\\ 0& 1& 0\end{array}\right),L_2=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ -1& 0& 0\end{array}\right),L_3=\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)$$

  

  der drei Lie-Algebra-Generatoren der Drehgruppe SO(3) konstruierter Vektor L = (L₁, L₂, L₃)^T ist.
• Also wird für die Konstruktion der Rotationsmatrix gemäß Gleichung (3.47) insbesondere folgende Kombination benötigt: $$\frac{1-\text{cos}2\zeta }{{\zeta }^2}=\frac{2{\text{sin}}^2\zeta }{{\zeta }^2}=\frac{2}{{\zeta }^2+{\left(1-{\zeta }^{2}B\left(\zeta \right)\right)}^2}$$

  
• Damit ergibt sich dann für die Rotationsmatrix gemäß Gleichung (3.47) der folgende Ausdruck

  
• Auch dieser Ausdruck ist eine rationale Funktion in ζ² und B(ζ), $\zeta =\frac{\sqrt{\theta \cdot \theta }}{2}$

  

  und seine Auswertung erfordert keine Auswertung transzendenter Funktionen und er enthält nur gerade Potenzen von ζ, so dass auch die Quadratwurzel in $\zeta =\frac{\sqrt{\theta \cdot \theta }}{2}$

  

  nicht berechnet werden muss.
• Bei Anwendung der Rotation, um den Vektor v in einem anderen Koordinatensystem darzustellen, lässt sich die Gleichung (3.50) wie folgt modifizieren

  
• Der Algorithmus 6 berechnet die Rotationsmatrix selbst bei gegebenen θ, θ² und B (ζ), z.B. approximiert via Gleichung (3.34). Es ist hier anzumerken, dass darin die Berechnung der Kombination $C=1-{\left(\frac{\theta }{2}\right)}^{2}B\left(\frac{\theta }{2}\right)$

  

  nicht erforderlich ist, wenn der Algorithmus 5 zuvor ausgeführt wurde, denn dieser benötigt ebenfalls die Kombination C, so dass das Ergebnis einfach im Rahmen von Algorithmus 5 abgespeichert und dann verwendet werden kann.

  
• Der folgende Algorithmus 7 ist eine leichte Modifikation von Algorithmus 6, die zeigt, wie die Koordinaten ^Gv eines Vektors v basierend auf der Gleichung (3.51) direkt im Bezugssystem G aus den Koordinaten ^Fv im Bezugssystem F berechnet werden können, anstatt die Rotationsmatrix zu berechnen.

  

  
• Beide Algorithmen, 6 und 7, sind auf irgendwelche Winkel θ anwendbar und könnten z.B für die effiziente und präzise Berechnung trigonometrischer Funktionen und Drehungen in Echtzeitanwendungen eingesetzt werden, z.B. zum Berechnen von Drehungen im Bereich der Computergrafik / Computer Vision und CAD.
• 9 zeigt ein Flussdiagramm des Verfahrens, bei dem mindestens ein Gyroskop in x-Richtung, in y-Richtung und in z-Richtung eine Messung der Winkelinkremente $\mathrm{\Delta }{\phi }_{\mathrm{1,2,3}}^k$

  

  vornimmt. Diese können aus dem Integral der Winkelgeschwindigkeit ω_1,2,3(t) über eine Zeit von einem vorhergehenden Messzeitpunkt t_k-1 bis zu einem aktuellen Messzeitpunkt t_k gebildet werden. Der Einfachheit halber wurden nur die Gleichungen für die Winkelinkremente dargestellt. Für die Einbeziehung auch der Geschwindigkeitsinkremente ist lediglich jedes Messergebnis eines Gyroskops als um das Messergebnis des Geschwindigkeitsinkrements in der entsprechenden Achsrichtung erweitert zu denken und dann sind korrespondierende Schritte unter Verwendung der Algorithmen für das Geschwindigkeits- und/oder Positionsupdate in analog zu den gezeigten Schritten für die Winkelinkremente durchzuführen. Die Beschreibung der 9 im Folgenden ist der Einfachheit ebenfalls auf die Winkelinkrementupdates beschränkt.
• Die drei Messwerte für die Winkelinkremente $\mathrm{\Delta }{\phi }_{\mathrm{1,2,3}}^k$

  

  können vektorisiert werden, um einen Drehwinkelinkrement-Vektor $\mathrm{\Delta }{\overrightarrow{\phi }}^k$

  

  zu erhalten.
• Damit können dann die einzelnen Drehwinkelinkremente Δθ^k rekursiv aus dem Exponenten $\mathrm{\Delta }{\theta }^k\left(\mathrm{\Delta }{\overrightarrow{\phi }}^{k-1},\mathrm{\Delta }{\overrightarrow{\phi }}^k\right)$

  

  mit der Verzögerung (Delay) Z^-1 zur Einbeziehung des vorhergehenden Drehwinkelinkrement-Vektors $\mathrm{\Delta }{\overrightarrow{\phi }}^{k-1}$

  

  berechnet werden.
• Die Berechnung der kumulierten Drehwinkelinkremente $\mathrm{\Delta }{\theta }_{out}^{I,k}$

  

  erfolgt mit dem Exponenten $\mathrm{\Delta }{\theta }_{out}^{I,k}=2\gamma \left(\frac{\mathrm{\Delta }{\theta }_{out}^{I,k-1}}{2},\frac{\mathrm{\Delta }{\theta }^k}{2}\right)$

  

  rekursiv mit einer Verzögerung (Delay) Z^-1 zur Einbeziehung des vorher berechneten kumulierten Drehwinkelinkrementes Δ $\mathrm{\Delta }{\theta }_{out}^{I,k-1}$

  

  aus den voher bestimmten einzelnen Drehwinkelinkrementen Δθ^k.
• Das letzte so berechnete kumulierte Drehwinkelinkrement bildet dann den kumulierten Drehwinkel $\mathrm{\Delta }{\theta }_{out}^I$

  

  als Messergebnis.

Claims (13)

1. Verfahren zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehwinkel des Objektes, wobei relative Drehwinkelinkremente Δθ^k an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen und für die jeweiligen Messzeitpunkte t_k zu kumulierten Drehwinkelinkrementen θ^k zusammengefasst werden, gekennzeichnet durch Berechnen der kumulierten Drehwinkelinkremente θ^k für die Messzeitpunkte t_k iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθ^k, die an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ^k-1 und der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

   

   $$

   

   mit ζ als Funktion des kumulierten Drehwinkelinkrementes θ^k-1 zum jeweils vorhergehenden Messzeitpunkt t_k-1, wobei ein Approximieren der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

   

   durch Padé-Approximationen erfolgt.
2. Verfahren nach Anspruch 1, gekennzeichnet durch Berechnen der kumulierten Drehwinkelinkremente θ^k für einen jeweiligen Messzeitpunkt t_k als Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3), wobei iterativ kumulierte Drehwinkelinkremente θ^k für die Messzeitpunkte t_k mit k = 1,..., n jeweils aus der Funktion des approximierten Exponenten der Baker-Campbell-Hausdorff-Formel mit jeweils einem zum Messzeitpunkt t_k gemessenen relativen Drehwinkelinkrement Δθ^k und dem für den vorhergehenden Zeitpunkt t_k-1 bestimmten kumulierten Drehwinkelinkrement θ^k-1 oder einem vorgegebenen Start-Drehwinkelinkrement θ⁰ als Eingangsgrößen bestimmt werden.
3. Verfahren nach Anspruch 2, gekennzeichnet durch Bestimmen eines kumulierten Drehwinkels $\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I,k}$

   

   für den Messzeitpunkt t_k durch das das zum letzten Zeitpunkt t_n des Zeitbereichs I, der die Messzeitpunkte k = 1, ..., n umfasst, bestimmte letzte kumulierte Drehwinkelinkrement θⁿ bestimmt ist.
4. Verfahren nach einem der Ansprüche 1 bis 3, dadurch gekennzeichnet, dass das kumulierte Drehwinkelinkrement θ^k für den Messzeitpunkt t_k mit der Funktion des approximierten Exponenten $\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

   

   der Baker-Hausdorff-Formel iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθ^k, die an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen wurden, und dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ^k-1 berechnet wird, wobei der Beitrag zum Exponenten in nullter Ordnung ${\mathrm{\gamma }}^{\left(0\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{{\mathrm{\theta }}^{k-1}}{2},$

   

   der Beitrag zum Exponenten in erster Ordnung ${\mathrm{\gamma }}^{\left(1\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}+\frac{{\mathrm{\theta }}^{k-1}}{2}\times \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}+\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^2\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)$

   

   $$

   

   und der Beitrag zum Exponenten in zweiter Ordnung ${\mathrm{\gamma }}^{\left(2\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)=\frac{1}{2}{\left(\frac{2}{{\mathrm{\theta }}^{k-1}}\right)}^2(1-3B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\left(1-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^{2}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right))\right)$

   

   $\left({\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^2{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\right)\frac{{\mathrm{\theta }}^{k-1}}{2}+\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)$

   

   $\left(1-{\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)}^{2}B\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\right)\left(\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\cdot \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}-{\left(\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\right)}^2\frac{{\mathrm{\theta }}^{k-1}}{2}\right)+\text{B}\left(\frac{{\mathrm{\theta }}^{k-1}}{2}\right)\frac{{\mathrm{\theta }}^{k-1}}{2}.$

   

   $\frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}\frac{{\mathrm{\theta }}^{k-1}}{2}\times \frac{\mathrm{\Delta }{\mathrm{\theta }}^k}{2}$

   

   lauten, mit der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}(1-\zeta cot\zeta )\simeq {\tilde{B}}_{Pade,N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\displaystyle {\sum }_{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m{\zeta }^{2m}}}{1+{\displaystyle {\sum }_{n=1}^{\left[\frac{N_0+1}{2}\right]}{b}_n{\zeta }^{2n}}},$

   

   $$

   

   wobei N₀ ein vorgegebener ganzzahliger Wert, [ ] die Floor-Funktion und die Parameter a_m und b_m von dem Wert N₀ abhängig vorgegebene Koeffizienten der Pade-Approximation sind, und sich das für den Zeitpunkt t_k kumulierte Drehwinkelinkrement $\frac{{\mathrm{\theta }}^k}{2}=\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\Delta {\mathrm{\theta }}^k}{2}\right)$

   

   $$

   

   jeweils aus der Summe der Beiträge ${\mathrm{\gamma }}^{\left({\text{m}}_{\mathrm{\alpha }}\right)}\left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\Delta {\mathrm{\theta }}^k}{2}\right)$

   

   zum Exponenten in der nullten, ersten und mindestens zweiten Ordnung m_α = 0, 1 und 2 ergibt.
5. Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, dass die Funktion des Exponenten $\gamma \left(\frac{{\mathrm{\theta }}^{k-1}}{2},\frac{\Delta {\mathrm{\theta }}^k}{2}\right)$

   

   die Summe der Argumente der Approximationen in der nullten, ersten, zweiten und dritten Ordnung m_α = 0, 1, 2 und 3 umfasst.
6. Verfahren nach einem der vorhergehenden Ansprüche, dadurch gekennzeichnet, dass zur Bestimmung des kumulierten Drehwinkels Δθ^I _out für den Zeitbereich $\left[t_{i-1};t_{i-1}+{\displaystyle {\sum }_{k=1}^n\Delta t_k}\right]$

   

   mit den gemessenen Drehwinkelinkrementen Δθ¹ bis Δθⁿ für die jeweiligen Zeitpunkte t_k mit k = 1,..., n als Eingangsgrößen das Start-Drehwinkelinkrement $\mathrm{\Delta }{\mathrm{\theta }}_{out}^{I\mathrm{,0}}$

   

   für k = 0 mit dem Wert Null vorgegeben wird und die kumulierten Drehwinkelinkremente θ^k iterativ mit k = 1 bis n aus dem für den vorhergehenden Iterationsschritt k-1 bestimmten oder vorgegebenen kumulierten Drehwinkelinkrement θ^k-1 und dem zum Zeitpunkt t_k, das zum jeweiligen Iterationsschritt k gehörig ist, gemessenen Drehwinkelinkrement Δθ^k als die transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

   

   $$

   

   umfassenden Funktion des Exponenten der Baker-Campbell-Hausdorff-Formel für die Lie-Rotationsgruppe SO(3) mittels einer Pade-Approximation der tranzendenten Funktion B(ζ) berechnet wird.
7. Verfahren nach einem der vorhergehenden Ansprüche, gekennzeichnet durch Berechnung einer Rotationsmatrix der Rotation zwischen einem ersten Koordinatensystem F des Objektes und einem nach Anwendung einer durch den kumulierten Drehwinkel θ^k beschriebenen Drehung erhaltenen zweiten Koordinatensystem G mit den Ausdrücken der mittels der Funktion B(ζ) und insbesondere deren Pade-Approximationen berechneten trigonometrischen Funktionen und Berechnung der translatorischen Geschwindigkeit und/oder Position des Objektes aus den in das zweite Koordinatensystem transformierten kumulierten Drehwinkeln θ^k.
8. Verfahren zur automatischen rechnergestützten Berechnung von trigonometrischen Funktionen von Winkeln θ^k oder Drehungen aus gegebenen Winkeln θ^k, gekennzeichnet durch Berechnen der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

   

   mit ζ als Funktion des gegebenen Winkels θ^k, wobei ein Approximieren der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

   

   durch Pade-Approximationen erfolgt, und Erhalten der trigonometrischen Funktion oder Drehung als Ergebnis der Berechnung der transzendenten Funktion.
9. Verfahren nach Anspruch 8, dadurch gekennzeichnet, dass die Approximation der transzendenten Funktion durch Auswerten der Gleichung ${\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)=\frac{1}{3}\frac{1+{\displaystyle {\sum }_{m=1}^{\left[\frac{N_0}{2}\right]}{a}_m{\zeta }^{2m}}}{1+{\displaystyle {\sum }_{m=1}^{\left[\frac{N_0}{2}\right]}{b}_m{\zeta }^{2n}}}$

   

   $$

   

   erfolgt, wobei N₀ ein vorgegebener ganzzahliger Wert und die Parameter a_m und b_m von dem Wert N₀ abhängig vorgegebene Koeffizienten der Pade-Approximation sind.
10. Verfahren nach Anspruch 9, gekennzeichnet durch Erhalten der trigonometrischen Funktionen als Ergebnis der Pade-Approximation der transzendenten Funktion durch $$\frac{{\text{sin}}_{N_0}2\zeta }{2\zeta }=\frac{1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)}{{\zeta }^2+{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2},$$

    

    und/oder $${\text{cos}}_{N_0}2\zeta =\frac{{\zeta }^2-{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2}{{\zeta }^2+{\left(1-{\zeta }^2{\tilde{B}}_{\text{Padé},N_0}\left(\zeta \right)\right)}^2}.$$

    
11. Messeinrichtung zum Bestimmen der Lage eines Objektes durch Ermittlung der Drehraten des Objektes, wobei die Messeinrichtung eine Sensoreinheit zum Messen von relativen Drehwinkelinkrementen Δθ^k für aufeinanderfolgende Zeitpunkte t_k und eine Datenverarbeitungseinheit hat, die eingerichtet ist, die an aufeinanderfolgenden Zeitpunkten t_k gemessenen relativen Drehwinkelinkremente Δθ^k zu einem kumulierten Drehwinkelinkrement θ^k für die jeweiligen Messzeitpunkte t_k zu kumulierten Drehwinkelinkrementen θ^k zusammenzufassen, dadurch gekennzeichnet, dass die Datenverarbeitungseinheit zum Berechnen der kumulierten Drehwinkelinkremente θ^k für die Messzeitpunkte t_k iterativ aus einer Anzahl n von relativen Drehwinkelinkrementen Δθ^k, die an aufeinanderfolgenden Messzeitpunkten t_k mit k = 1, ..., n gemessen wurden, dem jeweils für das vorhergehende Zeitintervall bestimmten kumulierten Drehwinkelinkrement θ^k-1 und der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

    

    mit ζ als Funktion des kumulierten Drehwinkelinkrementes θ^k-1 zum jeweils vorhergehenden Messzeitpunkt t_k-1 und zum Approximieren der transzendenten Funktion $B\left(\zeta \right)=\frac{1}{{\zeta }^2}\left(1-\zeta cot\zeta \right)$

    

    durch Pade-Approximationen eingerichtet ist.
12. Messeinrichtung zum Bestimmen der Lage eines Objektes und/oder zur automatischen rechnergestützten Berechnung von trigonometrischen Funktionen von Winkeln θ^k oder Drehungen aus gegebenen Winkeln θ^k mit einer Datenverarbeitungseinheit, dadurch gekennzeichnet, dass die Datenverarbeitungseinheit zur Durchführung der Schritte der Verfahren nach einem der Ansprüche 1 bis 10 eingerichtet ist.
13. Computerprogramm, umfassend Befehle, die bei der Ausführung des Computerprogramms durch eine Datenverarbeitungseinheit diese veranlassen, die Schritte des Verfahrens nach einem der Ansprüche 1 bis 10 auszuführen.

Images