四元数微分方程的推导
四元数微分方程的推导
由于载体的运动，四元数Q 是变量，即0123,,,q q q q 是时间的函数。
刚体绕瞬时转轴转过σ角，其角速度为：
.t
tb n ωσ--
=（式1）
设这个运载体坐标系（b 系）和地理坐标系（t 系）之间的变换四元数的三角形式为：
cos
sin
2
2
Q n σ
σ
-
=+（式2）
对式2求导可得：
.
.
11sin cos sin
22222dQ d n
n dt dt
σσσσσ-
-
=-⋅+⋅+（式3） 因为：
.0t tb d n n n n dt
ωσ-
----
=⨯=⨯=（式4） 1n n --
⋅=-（式5）
则有：
.
..
.
11sin cos 2222
(cos
sin
)
2
2
2
Q n n
n σσσσ
σ
σ
σ
--
-
=-⋅+⋅=⋅+（式6）
将式1和式2代入式6得：
.
12
t tb Q Q ω-
=（式7）
由于捷联惯性导航系统的惯性器件是直接固联在运载体上的，所以陀螺测量得到的角速度是沿运载体坐标系的绝对角速度，因此应用式7不方便，需要进行进一步变换。
因为：
*t b tb tb Q Q ωω=（式8）
**1Q Q Q Q ⋅=⋅=（式9）
式中b tb ω是沿运载体的角速度. 将式8、9代入式7得：
.
1
2b tb Q Q ω-
=⋅（式10）
将式10写成矩阵形式为：
001
231103222301332100b x b
y
b z q q q q q q q q q q q q q q q q q q q q ωωω---⎡⎤⎡⎤⎡⎤
⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦（式11） 00112233001020b b b x y
z b b b x
z
y b b b y
z
x b b b z
y
x
q q q q q q q q ωωωωωωωωωωωω⎡⎤---⎡⎤
⎡⎤
⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥
⎢⎥-⎢⎥⎢
⎥⎢⎥⎣⎦⎣⎦⎣⎦（式12） 式中，,,b b b x y z ωωω分别表示载体坐标系相对于地理坐标系沿各个轴向的角速度分量。
四元数解算姿态完全解析及资料汇总
四元数完全解析及资料汇总本文原帖出自匿名四轴论坛，附件里的资源请到匿名论坛下载：感谢匿名的开源分享，感谢群友的热心帮助。
说什么四元数完全解析其实都是前辈们的解析，小弟真心是一个搬砖的，搬得不好希望大神们给以批评和指正，在此谢过了。
因为本人是小菜鸟一枚，对，最菜的那种菜鸟······所以对四元数求解姿态角这么一个在大神眼里简单的算法，小弟我还是费了很大劲才稍微理解了那么一点点，小弟搬砖整理时也是基于小弟的理解和智商的，有些太基础，有些可能错了，大牛们发现了再骂过我后希望能够给与指正哈。
好，废话到此为止，开始说主体。
四元数和姿态角怎么说呢？先得给和我一样的小菜鸟们理一理思路，小鸟我在此画了一个“思维导图”（我承认我画的丑），四元数解算姿态首先分为两部分理解：第一部分先理解什么是四元数，四元数与姿态角间的关系；第二部分要理解怎么由惯性单元测出的加速度和角速度求出四元数，再由四元数求出欧拉角。
图1 渣渣思维导图在讲解什么是四元数时，小弟的思维是顺着说的，先由四元数的定义说起，说到四元数与姿态角间的关系。
但在讲解姿态解算时，小弟的思维是逆向的，就是反推回来的，从欧拉角一步步反推回到惯性器件的测量数据，这样逆向说是因为便于理解，因为实际在工程应用时和理论推导有很大差别。
实际应用时正确的求解顺序应该为图1中序号顺序，即1->2->3->…….但在笔者讲解姿态求解时思路是如图2的。
图2 逆向讲解思路大家在看四元数时最好结合着代码一块看，小弟看的是匿名四轴的代码，感觉写的非常好也非常清晰，粘出来大家一块观摩。
红色部分是核心代码，总共分为八个步骤，和图1中的八个步骤是一一对应的。
讲解介绍时也是和代码对比起来讲解的。
代码可以去匿名官网上下载，都是开源的，不是小弟的，所以小弟不方便加在附件中。
好的，下面搬砖开始！。
嘿咻嘿咻关于四元数的定义摘自秦永元的《惯性导航》，里面有非常好的讲解，大家可以直接看绪论和第九章就可以。
四元数法求解姿态动力学方程
四元数法求解姿态动力学方程四元数法是一种常用于求解机器人、飞行器等复杂机械系统的姿态动力学方程的数学方法。
在四元数法中，姿态状态由一个四元数表示，称为方向余弦四元数或旋转四元数。
四元数的基本概念是，它是一个由实部和三个虚部组成的数，形式为q=s+xi+yj+zk。
其中，s是实部，i、j、k是三个独立的虚部，满足如下关系：
i²=j²=k²=ijk=-1
四元数可以表示三维空间中的旋转。
当一个向量v围绕一个单位向量u旋转θ角时，对应的旋转四元数为：
q=cos(θ/2)+usin(θ/2)
四元数可以进行加减、乘除等运算，并具有一些特殊的性质，如四元数的模长为1，单位四元数可以表示任意旋转等。
在使用四元数法求解姿态动力学方程时，首先需要将姿态状态表示为一个四元数q。
然后，根据四元数的性质和基本运算，可以推导出四元数微分方程，用于描述姿态状态的变化：
dq/dt=0.5*q*ω
其中，ω表示机体角速度矢量，表示机体围绕三个坐标轴的旋转速度。
该方程可以用来计算姿态状态在时间上的变化。
通过四元数法求解姿态动力学方程，可以有效地避免欧拉角表示法中出现的万向锁现象和奇异性问题，从而提高了姿态控制的精度和可靠性。
03.四元数
03.四元数四元数单位四元数可以⽤来描述三维空间中的旋转，它既是紧凑的，也没有奇异性。
⼀个四元数描述唯⼀⼀个空间旋转，但是⼀个空间旋转可以由互为相反数的两个四元数表⽰。
1.四元数四元数的定义四元数是由四个元构成的数Q(q0,q1,q2,q3)=q0+q1i+q2j+q3kq0,q1,q2,q3是实数，i,j,k是满⾜如下条件的单位向量i⊗i=−1,j⊗j=−1,k⊗k=−1 ⾃⼰与⾃⼰进⾏四元数乘结果为-1。
i⊗j=k,j⊗k=i,k⊗i=j按照i→j→k→i...的顺序两两作四元数乘可以得到下⼀个。
j⊗i=−k,k⊗j=−i,i⊗k=−j按照i→j→k→i...相反的顺序两两作四元数乘可以得到下⼀个的负。
⊗表⽰四元数乘法四元数的表达⽅式⽮量式Q=q0+q其中q0称四元数Q的标量部分，q称四元数Q的⽮量部分。
q是三维空间中的⼀个向量。
复数式Q=q0+q1i+q2j+q3k可视为⼀个超复数，Q的共轭复数记为:Q∗=q0−q1i−q2j−q3k三⾓式Q=cos θ2+u sinθ2式中，u为单位向量，即旋转轴。
θ为实数，即绕单位向量u的旋转⾓度。
指数式Q=e u θ2式中，u和θ同上。
矩阵式Q=q0 q1 q2 q3四元数的⼤⼩—范数四元数的⼤⼩⽤四元数的范数来表⽰：||Q||=q20+q21+q22+q23||Q||=1，则Q成为规范化四元数。
描述旋转的四元数成为规范化四元数。
规范化四元数参与旋转运动时要作归⼀化。
四元数的运算加法和减法：对应实部和虚部相加/减[]Q=q0+q1i+q2j+q3kP=p0+p1i+p2j+p3k则Q±P=(q0±p0)+(q1±p1)i+(q2±p2)j+(q3±p3)k乘法:合并同类项(根据向量i,j,k向量相乘规则)标量乘a Q=aq0+aq1i+aq2j+aq3k其中a为标量四元数乘P⊗Q=(q0+q1i+q2j+q3k)⊗(p0+p1i+p2j+p3k)=(p0q0−p1q1−p2q2−p3q3)+(p0q1+p1q0+p2q3−p3q2)i+(p0q2+p2q0+p3q1−p1q3)j+(p0q3+p3q0+p1q2−p2q1)k =r0+r1i+r2j+r3k写成矩阵形式r0r1 r2 r3=p0−p1−p2−p3p1p0−p3p2p2p3p0−p1p3−p2p1p0q0q1q2q3=M(P)Q或r0r1 r2 r3=q0−q1−q2−q3q1q0q3−q2q2−q3q0q1q3q2−q1q0p0p1p2p3=M′(Q)P注意：四元数乘法不满⾜交换律P⊗Q=M(P)Q≠M′(P)Q=Q⊗P 四元数乘法满⾜分配律和结合律P⊗（Q+R）=P⊗Q+P⊗RP⊗Q⊗R=(P⊗Q)⊗R=P⊗(Q⊗R)此外，还有(P⊗Q)∗=Q∗⊗P∗证明，略。
P17四元数 导航原理 教学课件
22
2
2
cos sin n
22
四元数既表示了转轴方向，又表示了转角大小（转动四元数）
四元数 表示
四元数的表示
由一个实单位和三个虚数单位 i, j, k 组成的数
q 1 P1i P2 j P3k
或者省略 1，写成
q P1i P2 j P3k
i, j, k 服从如下运算公式：
四元数 组成部分
i, j, k 服从如下运算公式
i i j j k k 1 i j ji k
k
j k k j i k i i k j
i
j
q P1i P2 j P3k
λ 称作标量部分， P1i P2 j P3k
四元数的另一种表示法 q , P
称作矢量部分 P 泛指矢量部分
q q2 q1
cos 2
ห้องสมุดไป่ตู้sin
2
(cos
i
sin
j ) cos
2
sin
2
k
cos cos sin cos i sin sin j cos sin k
22
22
22
22
求方向余弦 映象方式1
以瞬时转轴映象形式给出 转动四元数的表达式并求 出合成转动四元数 第一次转时，映象形式的 q1 和非映象形式的 q1 是 一致的：
提示：四元数与刚体转动的关系
四元数基本性质 加减法
q P1i P2 j P3k
M v 1i 2 j 3k
1．四元数加减法
qM ( v) (P1 1 )i (P2 2 ) j (P3 3 )k
或简单表示为
q M v, P
四元数基本性质 乘法
2．四元数乘法
导航原理之四元数
Lecture 11 -- Quaternion
13
2.4 共扼和范数
共扼四元数的定义 -----两个四元数的标量部分相同，向 量部分相反
q P i P2 j P3 k 1
15
Outline
四元数的定义 四元数的运算 利用四元数进行旋转变换
利用四元数进行旋转合成
捷联惯导系统基于四元数的姿态微分方程求解
Lecture 11 -- Quaternion
16
3.1 矢量的旋转
如果矢量 R 相对固定坐标系旋转, 得到新矢 量 R’，且该旋转可以用四元数 q 描述， 则从 R=[Rx Ry Rz]T 到 R’=[R’x R’y R’z]T 的 变换可以用下述四元数运算表示:
2
3 2 P1 1 P2 v P3
M L (q)
M R (h)
Lecture 11 -- Quaternion
11
2.3 四元数乘法自定义函数
function [q1]=qmul(q, m) lm=q(1); p1=q(2); p2=q(3); p3=q(4); q1=[lm -p1 -p2 -p3 p1 lm -p3 p2 p2 p3 lm -p1 p3 -p2 p1 lm]*m; >> a=[1 2 2 3]'; >> b=[2 4 2 3]'; >> q=qmul(a,b)
“Here as he walked by on the 16th of October 1843, Sir William Rowan Hamilton in a flash of genius discovered the fundamental formula for quaternion multiplication
关于四元数的一些理解第二部分
距离上一部分关于四元数的理解写完已经有段不短的日子里，这次终于理解了所谓θ2的真谛。
这里我就一一道来。
首先经常看到网络中的关于四元数与旋转的关系中看到这样的描述;q=(cosθ2,℮n sinθ2)其中℮n为旋转轴的单位向量θ表示旋转的角度;有的甚至在后面描述为“四元数表达式的形式跟其旋转的角度θ以及旋转轴℮n有一定的关系”。
这种模糊不清的关系令人存疑。
笔者认为这只是大神们在解读解算程序语言时对于姿态解算中部分程序的误读造成的。
正如笔者在上一部分中提到的：若存在两个向量A,B其中B为单位向量且A,B的夹角为θ则有：B= −rcosθ+r℮n sinθ A−1=−1rr cosθ+r℮n sinθ−1−A∗=cosθ+℮n sinθ(−A∗)其中(−A∗)表的是一个方向与A相同的单位向量，而式子前段的cosθ+℮n sinθ部分，℮n很显然表示一个垂直于A，B的向量，θ角表示A，B两个向量之间的夹角。
其中℮n确实表示旋转轴且这个向量是一个单位向量，而θ的的确确是旋转的角度，这一点也能从推导过程中看出端倪。
那么问题究竟出在哪里呢？其实问题就在相乘的顺序上，通过阅读网上的博文我们注意到一般的文章中对于四元数表达的方式是这样的：存在一个单位四元数q=(cosθ2,℮n sinθ2);P是一个没有实部的单位四元数，表述为(0,v)(P的模为1);p′=qpq−1;好了，到了这里细心地读者已经发现一点问题了。
虽然这里，向量的表示换成了没有实部的四元数（例如p，p′）。
但是在这里我们发现，在变换向量的时候不但用q左乘向量p同时也用q−1右乘。
到这里笔者又要啰嗦些关于四元数的性质（尽可能详细的解读）。
如果您对这一部分有所了解同时时间有限的话可以直接略过！虽然笔者的上一篇文档提到了四元数的一部分内容，不过那只是笔者对于其本质在空间几何中的一些猜想。
在这里我就重新梳理一遍关于四元数的前世今生(大雾)。
在正式说明之前首先我希望明确一点，不论是用空间向量的方式表述还是用复数表达，他们都表达的是同一事物。
四元数——精选推荐
四元数四元数四元数四元数是最简单的超复数。
复数是由实数加上元素 i 组成，其中i^2 = -1 \,。
相似地，四元数都是由实数加上三个元素 i、j、k 组成，⽽且它们有如下的关系： i^2 = j^2 = k^2 = ijk = -1 \, 每个四元数都是 1、i、j 和 k 的线性组合，即是四元数⼀般可表⽰为a + bi + cj + dk \,。
⽬录四元数(Quaternions)是由威廉·卢云·哈密顿(William Rowan Hamilton, 1805-1865)在1843年爱尔兰发现的数学概念。
四元数的乘法不符合交换律(commutative law)，故威廉·卢云·哈密顿它似乎破坏了科学知识中⼀个最基本的原则。
明确地说，四元数是复数的不可交换延伸。
如把四元数的集合考虑成多维实数空间的话，四元数就代表著⼀个四维空间，相对於复数为⼆维空间。
四元数是除法环的⼀个例⼦。
除了没有乘法的交换律外，除法环与场是相类的。
特别地，乘法的结合律仍旧存在、⾮零元素仍有唯⼀的逆元素。
四元数形成⼀个在实数上的四维结合代数（事实上是除法代数），并包括复数，但不与复数组成结合代数。
四元数（以及实数和复数）都只是有限维的实数结合除法代数。
四元数的不可交换性往往导致⼀些令⼈意外的结果，例如四元数的 n-阶多项式能有多於 n 个不同的根。
四元数就是形如ai+bj+ck+d 的数a、b、c、d是实数i^2=j^2=k^2=-1ij=k ji=-k jk=i kj=-i ki=j ik=-j(a^2+b^2+c^2+d^2)的平⽅根称为四元数的模.例⼦假设：x = 3 + i \,y = 5i + j - 2k \,那么：x + y = 3 + 6i + j - 2k \,xy = \left( {3 + i} \right)\left( {5i + j - 2k} \right) = 15i + 3j - 6k + 5i^2 + ij - 2ik= 15i + 3j - 6k - 5 + k + 2j = - 5 + 15i + 5j - 5k \, 群旋转象在四元数和空间转动条⽬中详细解释的那样，⾮零四元数的乘法群在R3的取实部为零的拷贝上以共轭作⽤可以实现转动。
四元数,矢量运算规则,场论基础,并矢,算符,场量的Taylor展开,正交曲线坐标系,Delta函数
四元数quaternions复数对四则运算,代数运算,极限自封.四元数是复数的扩展.四元数有四个单元:k ,j ,i ,1.四元数定义dk cj bi a +++=α,其中R d ,c ,b ,a ∈ 另一四元数R 'd ,'c ,'b ,'a ,k 'd j 'c i 'b 'a ∈+++=β,则四元数加减法定义对应分量相加减;四元数乘法定义为)k 'd j 'c i 'b 'a )(dk cj bi a (++++++=αβk )'cb 'bc 'da 'ad (j )'bd 'db 'ca 'ac (i )'dc 'cd 'ba 'ab ()'dd 'cc 'bb 'aa (-+++-+++-+++---=四元数的单元间的运算规则: j ik ki ,i kj jk ,k ji ij ,1k j i 222=-==-==-=-===四元数加法适合结合律,交换律;,即)()(βγαγαβ=而一般βααβ≠.(βααβα=⇒∈R ) 对实数有效的运算规则对复数总有效,但对复数有效的运算规则对四元数不总有效,(如上述的乘法的交换律)！！！ 四元数的共轭: dk cj bi a :---=α,若dk cj bi a +++=α 性质:αβαβ=四元数的迹: R a 2:)(S ∈=+=ααα性质: )(S )(S )(S βαβα+=+四元数的模:R d c b a :)(N 2222∈+++==ααα性质: )(N )(N )(N βααβ⋅=,0)(N 0=⇔=αα证明: 0,oder ,00==⇒=βααβ00)(N 00)(N 00,und ,0=⇒⎭⎬⎫≠⇒≠=⇒=⇒≠=βααβαααβαααβ,即00,und ,0=⇒≠=βααβ,同理00,,0=⇒≠=αβαβund证明:若α是方程0)(N )(S x 2=+-αα的根,则α也是其根.因为,α是方程0)(N )(S x 2=+-αα的根0)()(2=++-⇒ααααα⇒=++-⇒0)()(2αααααα也是其根)四元数域内二次方程一般不止两个根,如最简单的方程1x 2-=就最少存在k ,j ,i ±±±6个根,实际上1x 2-=有无穷多个根,因为使1r q p 222=++成立的实数r ,q ,p 有无穷多个,而1)r q p ()rk qj pi (2222-=++-=++Halmiton 四元素体;第一个非交换体,1843 年 W.R.Hamilton 为建立三维复数空间,把复数x+iy 作为有序偶的实数,并定义规则,使i 在有明确意义: 4阶实方阵集H内方阵型如⎪⎪⎪⎪⎭⎫ ⎝⎛------a b c d b a d c c d a b d c b a ,令⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎪⎪⎭⎫ ⎝⎛--=⎪⎪⎪⎪⎭⎫ ⎝⎛=0110011010011001011001101111K ;J ,E ,I ,则集H 内任意方阵可唯一表为dK cJ bE aI +++,即}R d ,c ,b ,a |dK cJ bE aI {H ∈+++=,H 对矩阵减法封闭;且I K J E -===222,;J KE ,E JK ,K EJ ===J EK ,E KJ ,K JE -=-=-=,矩阵乘法在H 内封闭,故H 对矩阵加,乘法构成环;H 的元素个数>1;I 是H 的单位元,又因I )d c b a ()dK cJ bE aI )(dK cJ bE aI (2222+++=---+++,且当0≠+++dK cJ bE aI 时,d ,c ,b ,a 不全为零,故02222>+++d c b a ,所以H 中非零元在H 内存在逆元,综上所述H 是非交换体,常称H 为四元数环,称H 内的元为四元数Quaterion : t+xi+yj+zk,其中t为数量部分/纯量部分,xi+yj+zk 为向量部分.四元数系构成了以实数域为系数域的有限维可除代数,是向量代数和向量分析基础. 矢量运算规则两矢的内积:)b ,a cos(|b ||a |b a ∧=⋅ R V ,V →两矢的外积: )b ,a sin(|b ||a ||b a |∧=⨯, b ,a )b a ( ⊥⨯ V )V ,V (→ 物理意义: b ,a 两矢内积是功; b ,a 两矢外积的模是以b ,a两矢的为边平行四边形的面积. 故内积可交换,外积可反交换外积和内积的关系:)b ,a (sin |b ||a |))b ,a (cos 1(|b ||a |)b a (|b ||a ||b a |2222222222∧∧=-=⋅-=⨯ 即)b ,a sin(|b ||a ||b a |∧=⨯推论 22b a b b a a )b a ()b a (;b a )b a ()b a (-=⋅-⋅=+⋅-⨯=+⨯-四元数和两重积间的联系:两四元数k a j a i a 321++=α,k b j b i b 321++=β;两矢量)a ,a ,a (a 321=,)b ,b ,b (b 321= 间关系βα↔↔b ,a 两矢内积和四元数间的关系:两量积)Re()Re()(21)(21b a αββαβαβααββα-==+=+=⋅ ,即两矢内积b ,a 对应于四元数βα的实部.两矢外积和四元数间的关系:矢量内积)Im()Re()(21b a αβαβαβαβαβ=-=-=⨯ ,即两矢外积b a ⨯对应于四元数αβ的非实部.两矢内积,外积和四元数间的关系:αβ=⋅-⨯b a b a三矢内积)c b (a c )b a (:]c ,b ,a [ ⨯⋅≡⋅⨯=,R V ,V ,V → 物理意义: c ,b ,a 三矢的内积是以c ,b ,a三矢为边的平行六面体的体积性质:b )c a (a )b c (c )a b (b )a c (a )c b (c )b a (⋅⨯-=⋅⨯-=⋅⨯-=⋅⨯=⋅⨯=⋅⨯推论:0]q c ,p b ,r a []p c ,r b ,q a []r c ,q b ,p a [=⨯⨯⨯+⨯⨯⨯+⨯⨯⨯2]c ,b ,a []a c ,c b ,b a [ =⨯⨯⨯ 三矢外积c )b a ()c a (b )c b (a⋅-⋅=⨯⨯V )V ,V ,V (→c)b a (b )c a (c c c )b a b a b a (b b b )c a c a c a (c )b a b a b a (b )c a c a c a (c )b a b a b a (b )c a c a c a (c )b a b a b a (b )c a c a c a ()c b c b (a )c b c b (a )c b c b (a )c b c b (a )c b c b (a )c b c b (a c b c b c b c b c b c b a a a )c b (a 321332211321332211333221133322112332211233221113322111332211233223113112211233233113312212122131132332321⋅-⋅=⎪⎪⎪⎭⎫ ⎝⎛++-⎪⎪⎪⎭⎫ ⎝⎛++=⎪⎪⎪⎭⎫⎝⎛++-++++-++++-++=⎪⎪⎪⎭⎫ ⎝⎛---------=⎪⎪⎪⎭⎫ ⎝⎛---⨯⎪⎪⎪⎭⎫ ⎝⎛=⨯⨯推论0)b a (c )a c (b )c b (a=⨯⨯+⨯⨯+⨯⨯四矢内积:)c b )(d a ()d b )(c a (db cb d a ca )d c ()b a (⋅⋅-⋅⋅=⋅⋅⋅⋅=⨯⋅⨯R )V ,V ,V ,V (→)c b )(d a ()d b )(c a (a )d )c b (c )d b ((a ))d c (b ()d c ()b a (三矢外积三矢内积 ⋅⋅-⋅⋅=⋅⋅-⋅=⋅⨯⨯=⨯⋅⨯四矢外积:a ]b ,d ,c [b ]a ,d ,c [d ]c ,b ,a [c ]d ,b ,a [)d c ()b a (-=-=⨯⨯⨯ V )V ,V ,V ,V (→ a ]b ,d ,c [b ]a ,d ,c [a )b )d c ((b )a )d c (()d c ()b a (;d ]c ,b ,a [c ]d ,b ,a [d )c )b a ((c )d )b a (()d c ()b a (三矢外积三矢外积 -=⋅⨯-⋅⨯=⨯⨯⨯-=⋅⨯-⋅⨯=⨯⨯⨯推论c ]c ,b ,a []d ,b ,a [b ]c ,b ,a []c ,d ,a [a ]c ,b ,a []c ,b ,d [d 0]c ,b ,a [ ++=→≠ )d a )(c b ()c a )(d b ()}d c (b {a ⨯⋅-⨯⋅=⨯⨯⨯流线 等X 面/线 通量 环流量 散度 旋度 方向导数 梯度为形象描述矢量场)z ,y ,x (f 定义)z ,y ,x (f 的流线f.为形象描述标量场)z ,y ,x (ϕ定义)z ,y ,x (ϕ的等X 面/线.S d 为开/闭有向曲面S 上一面元,矢量f 在面元S d 上的元通量S d f d f⋅=Φ,面积分得矢量场)z ,y ,x (f 在曲面S 上的通量(标)⎰⋅=ΦSf S d fl d 为开/闭有向曲线l 上一面元,矢量f 在线元l d 上的元环量l d f d f⋅=Θ,线积分得矢量场)z ,y ,x (f 在曲线l 上的环量(标)⎰⋅=Θlfl d f 矢量场)z ,y ,x (f的散度(标):描述有源场源/汇强度. 正/负/零散度对应于源/汇/无源无汇闭合曲面S 包围体积V ∆,0V →∆时f 在S上的通量与V ∆比的极限称为矢量场)z ,y ,x (f 的散度V /S d f lim f div S 0V ∆⋅=⎰→∆矢量场)z ,y ,x (f的旋度(矢):描述有旋场旋涡强度和旋涡法矢方向. 旋度的法向分量的模的大小顺比于涡旋场旋涡程度.闭合曲线l 包围有向曲面S ∆,0S |S |→∆=∆ 时f 在l 上的环量与S ∆的比的极限称为矢量场)z ,y ,x (f 的旋度f rot 沿S∆法向的分量S /l d f lim )f rot (l 0S n ∆⋅=⎰→∆ 等效于0S →∆时S )f rot (l d f ∆⋅=⋅⎰标量场)z ,y ,x (ϕ的梯度(矢):描述标量场各点空间变化率及方向.某场点的梯度的方向是标量场变化最快的方向,其模是标量场单位长度的变化率.场)z ,y ,x (ϕ沿l d 向改变ϕd ,称dl d ϕ为ϕ沿l d 向的方向导数,dl d ϕ等于ϕ的梯度的l d 向分量l d )grad (d dld )grad (l ⋅=⇔=ϕϕϕϕ积分变换公式Gauss 定律: ⎰⎰∂⋅=⋅∇VV dS f n dV f(f 的散度对体积V 体积分 ←转换→ f 对V 的包面的闭面积分)Stokes 定律: ⎰⎰∂⋅=⋅⨯∇SS l d f dS n )f ((f 对有向曲线S ∂的闭线积分 ←转换→ f 的旋度对以S∂为边的有向曲面S 的面积分)Green 恒等式:⎰⎰⎰⎰⎰∂∂∂∂∂=⋅∇=⋅∇=∇⋅∇=∇+∇⋅∇V V V V V 2dS ndS n )(S d )(dV )(dV )(ψφψφψφψφψφψφ (n ∂∂:外法线方向导数)Green 定理:⎰⎰⎰⎰⎰∂∂∂∂∂-∂∂=⋅∇-∇=⋅∇-∇=∇-∇⋅∇=∇-∇VV V V V 22dS )n n (dS n )ˆˆ(S d )ˆˆ(dV )ˆˆ(ˆdV )ˆˆ(φψψφφψψφφψψφφψψφφψψφ ⎰⎰⨯∇=⨯∂V V dV A A S d (⎰⎰⎰⎰⎰⎰⨯∇⋅=⨯⋅⇔⨯∇⋅=⨯⋅∇=⨯⋅=⨯⋅∂∂∂VV V V V V dV )A (C )A S d (C dV )A (C dV )C A ()C A (S d )A S d (C) ⎰⎰∇=∂VVdV S d ψψ (⎰⎰⎰⎰⎰⎰∇⋅=⋅⇔∇⋅=⋅∇=⋅=⋅∂∂∂VVVVVVdV C )S d (C dV C dV )C (S d C )S d (C ψψψψψψ)⎰⎰∂=∇⨯SSl d S d ψψ(⎰⎰⎰⎰⎰⎰∂∂⋅=⋅∇⨯⇔⋅=⋅⨯∇=⋅⨯∇=⋅∇⨯SSSSSSl d C C )S d (l d C S d )C (S d )C (C )S d (ψψψψψψ)并矢 及其运算标量是零阶张量,矢量是一阶张量,二阶三维张量借助于直角坐标转动矩阵定义:矢量i T 在坐标架转动满足变换关系i mi m T R T =,坐标转动矩阵miR 即二阶张量.二阶张量ij T 满足变换关系ij nj mi mn T R R T =.由两矢B A ,并列放置且之间无运算则构成并矢B A,含9个分量,记为j i B A ,由于i A 和j B 分别满足:i mi m A R A =,j nj m B R B =,故并矢B A满足j i nj mi B A R R B A = ,故并矢是二阶张量的一种形式,显然三个矢量的并矢具有三阶张量的变换关系. ()⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛=≡><≡⊗≡z z z z z z y y y y y y x x x x x x z y x z y x v v g f g f g f g f g f g f g f g f g f g g g f f f g f g f g f g f || , 单位并矢(单位二阶张量)ij ij kk jj ii r δ=I =++=∇=I,性质:X X X =⋅I =I ⋅;(X 为矢量或算符); 2:∇=∇∇ I ; ϕϕ∇=⋅∇)(I ; )(:T Spur T T I ii ==; g f g f I ⋅=:;并矢-矢量点乘区分左右:右点乘p g f p g f )(⋅=⋅;左点乘)(f p g f p g⋅=⋅,这样,三矢外积可用并矢表示)()(p g g p f p g f -⋅=⨯⨯两二阶张量B A ,间的双点乘:ji ij B A B A = :(或))(()(:)(q f p g q p g f⋅⋅=)双点乘得到的标量是两矩阵积的迹.并矢的微分运算要注意是对那个张量进行的,一般需加括号.g f f g g f)()()(∇⋅+⋅∇=⋅∇g f g f g f)()()(∇⨯-⨯∇=⨯∇f r f r r f ⋅+⋅∇=⋅∇2)(22f r f r f +⋅∇=⋅∇)()(f r r f r r f r r f++⋅∇=⋅∇)()( ⎰⎰⎰∂∂⋅I =I ⋅=I ⋅∇VV V S d dS n dV⎰⎰∂I ⨯=I ⨯∇VV dS n dV⎰⎰∂=∇V V dS f n dV f⎰⎰∂⋅=⋅∇VVS d g f dV g f)()(根据以上矢量运算定理,可把Gauss 定理⎰⎰∂=∇VV n S d dV 和Stokes 定理⎰⎰∂=∇⨯SS l d n dS 的运算推广到对标量,矢量,张量的各种运算∇算符具有:矢量性和算符性.∇对矢量左/右/点/叉乘不可交换,矢量运算规则也适于∇,但需调整∇在结果中的位置,使等式左右量同型. f )g (g )f ()g f ( ⋅⨯∇-⋅⨯∇=⨯⋅∇ (第一式点和叉换位,取正;第二式交换第一式中的两矢量次序,取负)f )f ()f (2 ∇-⋅∇∇=⨯∇⨯∇ (按c )b a ()c a (b )c b (a⋅-⋅=⨯⨯写结果,再调整次序,使右端得矢量)ψϕϕψϕψ∇+∇=∇ )( f f )()f ( ⋅∇+⋅∇=⋅∇ϕϕϕf f )()f (⨯∇+⨯∇=⨯∇ϕϕϕ 2f 21f )f ()f f (f )f (f )f ( ∇-⋅∇=⋅∇-⋅∇=⨯⨯∇]g )f ()f (g []f )g ()g (f [)g f (∇⋅+⋅∇-∇⋅+⋅∇=⨯⨯∇ )f (g )g (f f )g (g )f ()g f (⨯∇⨯+⨯∇⨯+∇⋅+∇⋅=⋅∇3r =⋅∇ ; I r r =⊗∇≡∇; r e r r r ==∇;0e r r =⨯∇=⨯∇; r 2e r =⋅∇ 2r 3re r r r 1 -=-=∇ r e dr df r ˆdr df )r (f ˆ =∇=∇ ⎩⎨⎧=∞≠=-=-∇=-∇=∇)0(,)0(,0)(4ˆˆ1ˆ232r r r r e r r r r πδ Coulomb 定理的微分式:Green 函数|'r r |141)'r ,r (G ),r (4r e ˆ02r -==∇πεπδ标量场的梯度场无旋0)(≡∇⨯∇ϕ无旋场必可表为一标量场的梯度ϕ∇=⇒=⨯∇f 0f矢量场的旋度场无源0)f (≡⨯∇⋅∇ 涡旋场必可表为一矢量场的旋度A f 0f⨯∇=⇒=⋅∇a,0E ,k 为常矢a r )a ()r a ( =∇⋅=⋅∇ r r a r a ⋅=∇⋅ r a]e )e a (a [r 1e )a (r r r ⊥=⋅-=∇⋅ 533r r )r a (3r a r r )a ( ⋅-=∇⋅ 0a )r (r )a ()r a (=⋅⨯∇-⋅⨯∇=⨯⋅∇a 2r )a ()r (a )a (r a )r ()r a ()r a ()r a (r a =∇⋅-⋅∇+⋅∇-∇⋅=⨯⨯∇+⨯⨯∇=⨯⨯∇ 5333333r r )r a (3r a a )r r ()a (r r r r )a ()r r (a )r r a ( ⋅-=∇⋅+⨯∇⨯+∇⋅+⨯∇⨯=⋅∇ 0a )rr (r r )a ()r r a (333=⋅⨯∇-⋅⨯∇=⨯⋅∇ 35333333r a r )r a (3r r )a (a )r r ()a (r r r r )a ()r r (a )r r a ( -⋅=∇⋅-=∇⋅+⋅∇-∇⋅-⋅∇=⨯⨯∇33rr a )r r (a r 1a ⋅-=-⋅=∇⋅ r k i 0e E E ⋅=:E k i E ⋅=⋅∇;E k i )k i (e E )]r k i (e [E e E )e E (E r k i 0r k i 0r k i 0r k i 0⨯=⨯-=⋅∇⨯-=⨯∇-=⨯∇=⨯∇⋅⋅⋅⋅其中1:=⋅=r r r r e e I e e 故r r e e I r r I r r :'':''22==.Taylor 展开:...)(...'...''!)(...)(''!)('|'|,...,,+∂∂∂∂-++∂∂∂+∂∂-=-∑∑∑kj i k j i nk j i ji j i j i ii i rx x x x x x n r x x x x r x x r r r 111211112 其中k j i ,...,,取1,2,3; i x 代表直角坐标系的三个分量,注意:1 上式是对'r 展开; 2 对'r 的展开和对'r的展开相差一个负号. 曲线正交坐标系(Krummlinigen Koordinaten)三维空间里确定一点P 的位置需3个坐标321u ,u ,u .若P 点坐标在直角坐标系中表为)u ,u ,u (z ),u ,u ,u (y y ),u ,u ,u (x x 321321321===,则)z ,y ,x (u u ),z ,y ,x (u u ),z ,y ,x (u u 332211===,两坐标系等价.=i u 常数)3,2,1i (=的曲面是坐标面,他们的单位法向矢量为)3,2,1i (,e i =,其指向为iu 增加的方向.当过P 点的三坐标曲面两两垂直时,三坐标面的三交线也两两垂直,称此类坐标系为正交曲线坐标系.正交条件)j i (,0)u z )(u z ()u y)(u y ()u x )(u x (h ji j i j i 2ij≠=∂∂∂∂+∂∂∂∂+∂∂∂∂=.由i 31i i i 31i i i 31i i du u z dz ,du u ydy ,du u x dx ∑∑∑===∂∂=∂∂=∂∂=得2222)dz ()dy ()dx ()ds (++=232322222121j,i j i 2ij )du (h )du (h )du (h du du h ++==∑,其中2i 2i 2i 2ii 2i )u z ()u y()u x (h h ∂∂+∂∂+∂∂==,称)3,2,1i (,h i =为Lame 系数或度量因子.Delta 函数定义)a (f )x (f )a x (=-δ⎰性质1 偶函数)x ()x (δ=-δ2 采样性)a (f )x (f )a x (=-δ⎰3 函数下的面积⎩⎨⎧∉∈=-δ⎰])b ,c [a (,])b ,c [a (,dx )a x (bc 01; ⎩⎨⎧=∞≠=-δ)a x (,)a x (,)a x (04 缩放 )x (|a |)ax (δ=δ1证明)x (|a |)ax ()a (,|a |dz )z (|a |dz )z ()a (,a dz)z (dx )ax (|a ||a |a a a a δ=δ⇒⎪⎩⎪⎨⎧<δ=-δ>δ=δ⎰⎰⎰⎰∆∆-∆∆-∆∆-∆∆-100 5 若)(x f 为连续函数,且∆为包含a 电的任意长度区间,则a )x (g |)]x ('g /)x (f [dx ]a )x (g [)x (f =∆=-δ⎰证明dy a y g g dx dy x g dx dx x g dy a x g y a y g x a x g y )](['1)('1)(')()()(11+=→⎪⎭⎪⎬⎫=→=→-=+=→-=-- 若)x (f 为单值连续函数,且有N 个过零点N ,...,,i ,x i 321=,则a x g x g x f a g g a g f dy a y g g y a y g f dx a x g x f =--∆--∆==++=-⎰⎰)(1111|)('1)()](['1)]([)](['1)()]([])([)(δδ 6 复合函数∑-δ=δii i )x x (|x /)x (df |)]x (f [1 证明: )x (f y =单值连续,则)x (f 在每个过零点的邻域内可逆;且)x (g 为任意品优(gutartig)函数,则dy )]y (f ['f 'x dy )x ('f 'x dy )]'y (f ['x )y (f x )x (f y 11111---=→=→=→=→=⎰∑∑∑∑⎰⎰+∞∞---∆+∆---+∞∞--δ===δ=δii i ii i iix x dx )x x (|)x ('f |)x (g |)x ('f |)x (g )](f ['f )](f [g dy )]y (f ['f ]y [)]y (f [g dx )]x (f [)x (g i i1101011111 被积函数须相等,再由)x (g 的任意性,得∑-δ=δii i )x x (|x /)x (df |)]x (f [1 6 导数x)a x (:)a x ('∂-δ∂=-δ则)a ('f dx )a x (')x (f -=-δ⎰∆证明)a ('f dx )x ('f )a x ()x (df )a x (|)]a x ()x (f [)a x (d )x (f dx )a x (')x (f =-δ-=-δ--δ=-δ=-δ⎰⎰⎰⎰∆∆=+∆-∆∆∆ 07)a (f )(dx )a x ()x (f )n (n )n (1-=-δ⎰∆证明)a (f )(dx )x (f )a x ()(dx )x (f )a x ()(...dx )x (''f )a x ()x ('df )a x (|)]a x ()x ('f [)a x (d )x ('f dx )x ('f )a x ()x (df )a x (|)]a x ()x (f [)a x (d )x (f dx )a x ()x (f )n (n )n (n)m ()m n (m)n ()n ()n ()n ()n ()n ()n ()n ()n (1112202211011-=-δ-=-δ-==-δ=-δ+-δ-=-δ-=-δ-=-δ--δ=-δ=-δ⎰⎰⎰⎰⎰⎰⎰⎰⎰∆∆-∆-∆-=+∆-∆-∆-∆-∆-=+∆-∆-∆-∆⎰∆=-=-δa)x (g |]})x ('g )x (f [dx d )x ('g {dx ]a )x (g [')x (f 1 三维δ函数⎩⎨⎧=≠=-δ-δ-δ=-δ)r r (,)r r (,)z z ()y y ()x x ()r r (00000010曲线系下的三维δ函数3213020103020100h h h )u u ()u u ()u u (|)u ,x (J |)u u ()u u ()u u ()r r (i i -δ-δ-δ=-δ-δ-δ=-δ ,(其中)u ,x (J i i 为Jacobi 行列式)柱坐标系下)z z ()()()r r (00001-δϕ-ϕδρ-ρδρ=-δ球坐标系下)()()r r (sin r )r r (000201ϕ-ϕδθ-θδ-δθ=-δ注意:n 维δ函数的量纲为n m -,即n -米δ函数的逼近钟形曲线: 2201xa a lim )x (a +π=δ→ Gauss 曲线;)x n exp(n lim )x (n 220π-=δ→sinc 函数: )kx (c sin k lim x kx sin lim )x (k k π=π=δ∞→∞→1sinc 函数平方: )kx (c sin k lim kx kx sin lim )x (k k 2221π=π=δ∞→∞→ 复指函数:⎰⎰+∞∞-+∞∞-π=±π=δdk )kx cos(dk )ikx exp()x (121盒子函数:∑+∞-∞==δn )L /inx exp(L )x (21*********************.cn。
四元数微分方程 runge kutta
四元数微分方程 Runge-Kutta一、概述四元数微分方程是描述四元数随时间变化的数学模型，它在多个领域中都有着重要的应用。
在数值计算中，为了求解四元数微分方程的数值解，常常需要使用Runge-Kutta方法。
本文将探讨四元数微分方程和Runge-Kutta方法的理论基础和应用。
二、四元数微分方程1.四元数的定义四元数是由一个实部和三个虚部组成的超复数。
一般地，四元数可以表示为q=w+xi+yj+zk，其中w、x、y、z分别代表实部和三个虚部的系数，i、j、k分别是四元数单位虚数，并满足以下性质：i^2=j^2=k^2=ijk=-12.四元数微分方程的一般形式四元数微分方程的一般形式可以表示为：dq/dt=f(q,t)其中q表示四元数，t表示时间，f(q,t)表示四元数随时间t的变化率。
求解四元数微分方程通常需要利用数值方法。
三、Runge-Kutta方法Runge-Kutta方法是一种常用的数值求解微分方程的方法，它通过递推的方式来逼近微分方程的解。
对于四元数微分方程，在使用Runge-Kutta方法时，需要考虑四元数的乘法和加法运算。
1.四阶经典Runge-Kutta方法四阶经典Runge-Kutta方法是最常用的一种Runge-Kutta方法，其递推公式为：k1=h*f(q_n,t_n)k2=h*f(q_n+0.5*k1,t_n+0.5*h)k3=h*f(q_n+0.5*k2,t_n+0.5*h)k4=h*f(q_n+k3,t_n+h)q_{n+1}=q_n+1/6*(k1+2*k2+2*k3+k4)其中h为步长，f(q_n,t_n)为四元数在时刻t_n处的变化率。
2.四元数的乘法和加法运算在使用Runge-Kutta方法时，需要考虑四元数的乘法和加法运算。
四元数的乘法运算满足结合律和分配律，但不满足交换律。
在计算k1、k2、k3和k4时，需要考虑四元数乘法的性质。
四、应用实例为了说明四元数微分方程和Runge-Kutta方法的应用，我们可以考虑一个具体的实例，比如飞行器姿态控制系统中的四元数微分方程。
P17四元数微分方程求解
姿态实时计算 概述
四元数的实时计算
b b b ω Eb = ω ib − ω iE 因假定“数学平台”跟踪地理坐标系，因 因假定“数学平台”跟踪地理坐标系，
此 四元数微分方程： 四元数微分方程：
ɺ q (t ) = (Ω b − Ω b )q (t ) ib iE
增量算法 四元数
2、四元数微分方程的计算： 、四元数微分方程的计算：
T X (t + T ) = X (t ) + [ K 1 + K 2 ] 2
数值积分 2阶 四元数 阶
（2）四元数微分方程 ）
K 1 = Ω b (t )q(t )
Y = q(t ) + TK1
K 2 = Ω b (t + T )Y T q (t + T ) = q (t ) + ( K 1 + K 2 ) 2
四元数的实时计算因假定数学平台跟踪地理坐标系因其中i为单位四元数四阶龙格库塔积分四元数微分方程1一阶龙格库塔法rungekutta一个矩阵微分方程
四元数精确解1 姿态计算 四元数精确解
三、四元数微分方程式及其解 由第一章，四元数微分方程式： 由第一章，四元数微分方程式：
ɺ q = Ωb q
对
Ω b 的类似处理
ɺ X (t ) = f [ X (t ), ω (t )]
当初始条件已知，其一阶龙格 库塔的解为 库塔的解为: 当初始条件已知，其一阶龙格-库塔的解为
X (t + T ) = X (t ) + Tf [ X (t ), ω (t )]
数值积分 1阶 四元数 阶
（2）四元数微分方程 ）
或
ɺ q = Ωb q ɺ − ω x / 2 − ω y / 2 − ω z / 2 λ λ 0 ɺ ω / 2 ω z / 2 − ω y / 2 P1 0 P1 = x ɺ P2 ω y / 2 − ω z / 2 ω x / 2 P2 0 ɺ ωz / 2 ωy / 2 −ωx / 2 0 P3 P3 q(t + T ) = q(t ) + TΩ b (t )q(t )