航空機の運動方程式
何よりもまず結論から示すと、
並進運動
$
\begin{pmatrix} F_x \\ F_y \\ F_z \end{pmatrix}
= m
\begin{pmatrix} \dot{u}+qw-rv \\ \dot{v}+ru-pw \\ \dot{w}+pv-qu \end{pmatrix}
$
回転運動
$
\begin{pmatrix} \bar{L} \\ \bar{M} \\ \bar{N} \end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p}-(I_y-I_z)qr-I_{xy}(rp-\dot{q})-I_{yz}(q^2-r^2)-I_{zx}(\dot{r}+pq) \\
I_y \dot{q}-(I_z-I_x)rp-I_{xy}(qr-\dot{p})-I_{yz}(\dot{r}+pq)-I_{zx}(r^2-p^2) \\
I_z \dot{r}-(I_x-I_y)pq-I_{xy}(p^2-q^2)-I_{yz}(\dot{q}+rp)-I_{zx}(\dot{p}-qr)
\end{pmatrix}
$
以下は、この式の導出から変形等を示す。
使用する記号等
機体固定座標系
x軸：機首方向 y軸：右翼方向 z軸：機体の腹方向 の右手系
重心：C
地面固定座標 $x_e$ 軸 , $y_e$ 軸 , $z_e$ 軸 earth
原点：O
オイラー角 $(\phi ,\theta ,\psi)$
$\phi$ ： $x_e$ 軸周りの角度 roll角
$\theta$ ： $y_e$ 軸周りの角度 pitch角
$\psi$ ： $z_e$ 軸周りの角度 yaw角
機体の質量 $m$
機体の一部の微小質量 $∆ m$
重心の位置ベクトル $r_c$
微小質量の位置ベクトル $r_m$
微小質量の重心からの位置ベクトル $r$
速度 $V= (u ,v ,w)$
角速度（機体の重心C周りの回転ベクトル） $\omega = (p ,q ,r)$
全外力 $F= (X ,Y ,Z)$
モーメント $(L ,M ,N)$
角運動量 $H= (H_x ,H_y ,H_z )$
x軸,y軸,z軸に関する慣性能率 $(I_{xx} ,I_{yy} ,I_{zz} )$
x-y面,y-z面,z-x面に関する慣性乗積 $(I_{xy} ,I_{yz} ,I_{zx} )$
迎角 $\alpha$：$\tan \alpha = \dfrac{w}{u}$
横滑り角 $\beta$：$\sin \beta =\dfrac{v}{V}$
飛行経路角 $\gamma$：$\gamma=\theta-\alpha$
上反角 $\Gamma$：$\Gamma=\dfrac{-z_e'}{ \sqrt{ {x_e'}^2 + {y_e'}^2 } }$
舵角 $(\delta_a ,\delta_e ,\delta_r )$ : エルロン角、エレベータ角、ラダー角
$u',v',w'=f( u ,v ,w, p ,q ,r, \theta ,\phi,\psi , m ,T ,\rho ,V ,S , C_x ,C_y ,C_z )$
$p',q',r'=f( p ,q ,r, I_x ,I_y ,I_z ,I_{xz} , ρ ,V ,S ,c ,b , C_l ,C_m ,C_n )$
回転座標系の時間微分の仕方
通常、固定座標系においては、
$x_e$を微分すると、$\dfrac{dx_e}{dt}$ となるが、
角速度 $\omega$ で回転する回転座標系を基準に考える場合、
$x_e$を微分すると、$\dfrac{dx_e}{dt}=\dfrac{dx}{dt}+ \omega \times x$ となる。
航空機の運動方程式の導出（並進運動方程式）
重心の位置ベクトル $r_c$
微小質量の位置ベクトル $r_m$
重心から微小質量の位置ベクトル $r$ の関係は、
$r_m=r_c+r$
微分すると
$\dfrac{dr_m}{dt}=\dfrac{dr_c}{dt}+\bigl( \dfrac{dr}{dt}+\omega \times r \bigr)$
$r$は時間によって変化することはない定数なので、$\dfrac{dr}{dt}=0$
$\dfrac{dr_m}{dt}=\dfrac{dr_c}{dt}+ \omega \times r $
$\dfrac{dr_c}{dt}$ は機体重心の速度ベクトルであり、機体座標系の成分で表すと
$\dfrac{dr_c}{dt}=v_c=( u ,v ,w )$ なので、
$\dfrac{dr_m}{dt}=v_c+ \omega \times r $
（速度の式）
もう一回微分すると
$
\begin{align}
\dfrac{d^2 r_m}{dt^2}&=( \dfrac{dr_c}{dt}+\omega \times r) + \Bigl( \dfrac{d}{dt} (\omega \times r)+\omega \times (\omega \times r) \Bigr) \nonumber \\
&=( \dfrac{dr_c}{dt}+\omega \times r)+\Bigl( \dfrac{d \omega}{dt} \times r+ \omega \times \dfrac{dr}{dt} \Bigr) + \omega \times ( \omega \times r ) \nonumber \\
&=\dfrac{dr_c}{dt}+\omega \times v_c + \dfrac{d \omega}{dt} \times r +\omega \times ( \omega \times r ) \nonumber
\end{align}$
（加速度の式）
微小質量に係る外力を $∆F$ とすると、
ニュートンの第2法則 $F=m \alpha $ より、
$
\begin{align}
∆F&=∆m \dfrac{d^2r_m}{dt^2} \nonumber \\
&=∆m \Bigl(\dfrac{dr_c}{dt}+\omega \times v_c + \dfrac{d \omega}{dt} \times r +\omega \times ( \omega \times r ) \Bigr) \nonumber
\end{align}
$
積分して、機体全体の質量にして考えると
$
\begin{align}
F&=\int dm \Bigl(\dfrac{dr_c}{dt}+\omega \times v_c + \dfrac{d \omega}{dt} \times r +\omega \times ( \omega \times r ) \Bigr) \nonumber \\
&=\int dm \Bigl(\dfrac{dr_c}{dt}+\omega \times v_c \Bigr) +\dfrac{d \omega}{dt} \times \int r dm+\omega \times ( \omega \times \int r dm ) \nonumber
\end{align}$
rは時間によって変化することはない定数なので、$\int r dm=0$
よって、
$
\begin{align}
F&=\int dm \Bigl(\dfrac{dr_c}{dt}+\omega \times v_c \Bigr) \nonumber \\
&=m \Bigl(\dfrac{dr_c}{dt}+\omega \times v_c \Bigr) \nonumber
\end{align}$
成分表示をすると
$
\begin{align}
\begin{pmatrix} F_x \\ F_y \\ F_z \end{pmatrix}
&= m \left\{ \dfrac{d}{dt}
\begin{pmatrix} u \\ v \\ w \end{pmatrix}
+
\begin{pmatrix} p \\ q \\ r \end{pmatrix}
\times
\begin{pmatrix} u \\ v \\ w \end{pmatrix}
\right\}
\nonumber \\
&=m \left\{
\begin{pmatrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{pmatrix}
+
\begin{pmatrix} qw-rv \\ ru-pw \\ pv-qu \end{pmatrix}
\right\}
\nonumber \\
&=m
\begin{pmatrix} \dot{u}+qw-rv \\ \dot{v}+ru-pw \\ \dot{w}+pv-qu \end{pmatrix}
\nonumber \\
\end{align}$
（並進運動の方程式）
航空機の運動方程式の導出（回転運動方程式）
回転運動のモーメントを考える。
微小質量$∆m$にかかるモーメントは、
$
\begin{align}
r \times ∆F &= r \times ∆m \bigl( \dfrac{dv_c}{dt}+\omega \times v_c+ \dfrac {d \omega}{dt} \times r + \omega \times (\omega \times r) \bigr) \nonumber \\
&=∆m \Bigl(r \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + r \times \bigl( \dfrac{d \omega}{dt}\times r \bigr) + r \times \bigl( \omega \times ( \omega \times r ) \bigr) \Bigr) \nonumber
\end{align}$
行列の外積の計算の仕方 $[A \times (B \times C)=(A・C)・B-(A・B)・C ]$ なので
$
\begin{align}
r \times ∆F &= ∆m \Bigl(r \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + (r・r)\dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+r \times \bigl( (\omega・r)\omega - (\omega・\omega)r \bigr) \Bigr) \nonumber \\
&=∆m \Bigl(r \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + (r・r)\dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+(\omega・r)・(\omega \times r)-(\omega ・\omega)・(r \times r) \nonumber
\end{align}$
$r \times r = 0$ なので、
$r \times ∆F = ∆m \Bigl(r \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + (r・r)\dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt} \bigr)r+(\omega・r)(\omega \times r) \Bigr) $
積分して、機体全体の質量にして考えると
$
\begin{align}
r \times F &= \int dm \Bigl(r \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + (r・r)\dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+(\omega・r)(\omega \times r) \Bigr) \nonumber\\
&= \int r dm \times \bigl( \dfrac{dv_c}{dt}+\omega \times v_c \bigr) + \int dm \Bigl( (r・r) \dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+(\omega・r)(\omega \times r) \Bigr) \nonumber\\
&= 0 + \int dm \Bigl( (r・r) \dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+(\omega・r)(\omega \times r) \Bigr) \nonumber\\
&= \int dm \Bigl( (r・r) \dfrac{d \omega}{dt}- \bigl( r・\dfrac{d \omega }{dt}\bigr)r+(\omega・r)(\omega \times r) \Bigr) \nonumber
\end{align}$
成分表示をすると
$
\begin{pmatrix} \bar{L} \\ \bar{M} \\ \bar{N} \end{pmatrix}
= \int dm \left\{
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
･
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
\right\}
\dfrac{d}{dt}
\begin{pmatrix} p \\ q \\ r \end{pmatrix}
-
\left\{
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
･
\dfrac{d}{dt}
\begin{pmatrix} p \\ q \\ r \end{pmatrix}
\right\}
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
+
\left\{
\begin{pmatrix} p \\ q \\ r \end{pmatrix}
･
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
\right\}
\left\{
\begin{pmatrix} p \\ q \\ r \end{pmatrix}
\times
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
\right\}
$
$
= \int dm \left\{
(x^2+y^2+z^2)
\begin{pmatrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{pmatrix}
-
(x \dot{p}+y \dot{q}+z \dot{r})
\begin{pmatrix} x \\ y \\ z \end{pmatrix}
+
(px+qy+rz)
\begin{pmatrix} yr-zq \\ zp-xr \\ xq-yp \end{pmatrix}
\right\}
$
$
= \int dm \left\{
\begin{pmatrix} (x^2+y^2+z^2) \dot{p} \\ (x^2+y^2+z^2) \dot{q} \\ (x^2+y^2+z^2) \dot{r} \end{pmatrix}
-
\begin{pmatrix} x^2 \dot{p}+xy \dot{q}+zx \dot{r} \\ xy \dot{p}+y^2 \dot{q}+yz \dot{r} \\ zx \dot{p}+yz \dot{q}+z^2 \dot{r} \end{pmatrix}
+
\begin{pmatrix}
xyrp-zxpq+y^2 qr-yzq^2 +yzr^2 -z^2 qr\\
zxp^2 -x^2 rp+yzpq-xyqr+z^2 rp-zxr^2 \\
x^2 pq-xyp^2 +xyq^2 -y^2 pq+zxqr-yzrp
\end{pmatrix}
\right\}
$
$
= \int dm
\begin{pmatrix}
(y^2+z^2) \dot{p} +(y^2-z^2) qr-xy(rp- \dot{q})-yz(q^2-r^2)-zx( \dot{r}+pq) \\
(z^2+x^2) \dot{q} +(z^2-x^2) rp-xy(qr+ \dot{p})-yz( \dot{r}-pq)-zx(r^2-p^2) \\
(x^2+y^2) \dot{r} +(x^2-y^2) pq-xy(p^2-q^2)-yz( \dot{q}+rp)-zx( \dot{p}-qr)
\end{pmatrix}
$
$
=
\begin{pmatrix}
\int (y^2+z^2) dm \dot{p} +\int (y^2-z^2) dm qr-\int xy dm (rp- \dot{q})-\int yz dm (q^2-r^2)-\int zx dm ( \dot{r}+pq) \\
\int (z^2+x^2) dm \dot{q} +\int (z^2-x^2) dm rp-\int xy dm (qr+ \dot{p})-\int yz dm ( \dot{r}-pq)-\int zx dm (r^2-p^2) \\
\int (x^2+y^2) dm \dot{r} +\int (x^2-y^2) dm pq-\int xy dm (p^2-q^2)-\int yz dm ( \dot{q}+rp)-\int zx dm ( \dot{p}-qr)
\end{pmatrix}
$
ここで、
\begin{cases}
\int (y^2+z^2) dm =I_x ，\int (z^2+x^2) dm =I_y ，\int (x^2+y^2) dm =I_z \\ \\
\int xy dm =I_{xy}=I_{yx} ，\int yz dm =I_{yz}=I_{zy} ，\int zx dm =I_{zx}=I_{xz}
\end{cases}
とすると、
$
\begin{pmatrix} \bar{L} \\ \bar{M} \\ \bar{N} \end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p} + (I_y-I_z) qr-I_{xy} (rp- \dot{q})-I_{yz} (q^2-r^2)-I_{zx} ( \dot{r}+pq) \\
I_y \dot{q} + (I_z-I_x) rp-I_{xy} (qr+ \dot{p})-I_{yz} ( \dot{r}-pq)-I_{zx} (r^2-p^2) \\
I_z \dot{r} + (I_x-I_y) pq-I_{xy} (p^2-q^2)-I_{yz} ( \dot{q}+rp)-I_{zx} ( \dot{p}-qr)
\end{pmatrix}
$
（回転運動の方程式）
左右対称な機体の場合の回転運動方程式
もしも機体が 左右対称形 であった場合 $I_{xy}=I_{yz}=0$ となる。
そのため、回転運動の方程式は
$
\begin{pmatrix} \bar{L} \\ \bar{M} \\ \bar{N} \end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p} + (I_y-I_z) qr-I_{zx} ( \dot{r}+pq) \\
I_y \dot{q} + (I_z-I_x) rp-I_{zx} (r^2-p^2) \\
I_z \dot{r} + (I_x-I_y) pq_{zx} ( \dot{p}-qr)
\end{pmatrix}
$
この時、
$\bar{M}$ の第２項の $(I_z-I_x) rp$ を「ロール・ヨー・カップリング」と言い、
$\bar{N}$ の第２項の $(I_x-I_y) pq$ を「ロール・ピッチ・カップリング」と言う。
「ロール・ヨー・カップリング」
ロール運動中にヨー運動が加わるとピッチングモーメントが発生する。
「ロール・ピッチ・カップリング」
ロール運動中にピッチ運動が加わるとヨーイングモーメントが発生する
外力と（並進運動方程式の左辺）について
航空機に働く外力は （重力）と（推力）と（揚力）と（抗力）です。
（揚力）と（抗力）は、合わせて（空気力）と言うことができます。
つまり、航空機に働く外力は （重力）と（推力）と（空気力）です。
１ 重力
$W=(-mg \sin \theta , mg \cos \theta \sin \phi , mg \cos \theta \cos \phi)$
２ 推力
$T=(T \cos i_T , 0 , -T \sin i_T)$
３ 空気力
$( \dfrac{1}{2}C_x \rho SV^2 , \dfrac{1}{2}C_y \rho SV^2 , \dfrac{1}{2}C_z \rho SV^2 )$
よって、外力Fは
$
\begin{pmatrix} F_x \\ F_y \\ F_z \end{pmatrix}
=
\begin{pmatrix}
-mg \sin \theta + T \cos i_T + \dfrac{1}{2}C_x \rho SV^2 \\
mg \cos \theta \sin \phi + \dfrac{1}{2}C_y \rho SV^2 \\
mg \cos \theta \cos \phi -T \sin i_T + \dfrac{1}{2}C_z \rho SV^2
\end{pmatrix}
$
なので、航空機の並進運動の運動方程式
$
\begin{pmatrix} F_x \\ F_y \\ F_z \end{pmatrix}
= m
\begin{pmatrix} \dot{u}+qw-rv \\ \dot{v}+ru-pw \\ \dot{w}+pv-qu \end{pmatrix}
$は、
$
\begin{pmatrix}
-mg \sin \theta + T \cos i_T + \dfrac{1}{2}C_x \rho SV^2 \\
mg \cos \theta \sin \phi + \dfrac{1}{2}C_y \rho SV^2 \\
mg \cos \theta \cos \phi -T \sin i_T + \dfrac{1}{2}C_z \rho SV^2
\end{pmatrix}
= m
\begin{pmatrix} \dot{u}+qw-rv \\ \dot{v}+ru-pw \\ \dot{w}+pv-qu \end{pmatrix}
$
と書くことができ、これを整理すると
$
\begin{pmatrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{pmatrix}
= m
\begin{pmatrix}
-qw+rv-g \sin \theta + \dfrac{T}{m} \cos i_T + \dfrac{1}{2m}C_x \rho SV^2 \\
-ru+pw+g \cos \theta \sin \phi + \dfrac{1}{2m}C_y \rho SV^2 \\
-pv+qu+g \cos \theta \cos \phi -\dfrac{T}{m} \sin i_T + \dfrac{1}{2m}C_z \rho SV^2
\end{pmatrix}
$
これも（並進運動方程式）
モーメント（回転運動方程式の左辺）について
航空機に働くモーメントは、
（空気力によるモーメント）と（エンジンジャイロモーメント）に分けられます。
１ 空気力によるモーメント
平均空力翼弦を $\bar{c}$ 、翼幅を $b$ 、空力モーメント係数を$(C_l , C_m , C_n)$ とすると、
（空気力によるモーメント）は $\left( \dfrac{1}{2}C_l \rho SV^2b , \dfrac{1}{2}C_m \rho SV^2 \bar{c} , \dfrac{1}{2}C_n \rho SV^2b \right)$
２ エンジンジャイロモーメント
エンジンによる慣性モーメントを $I_R$、エンジンによる回転角速度を $ω_R$とすると、
（エンジンジャイロモーメント）は $\left(0 , -I_R \omega_R･r , I_R \omega_R･q \right)$
航空機に働くモーメントは
$
\begin{pmatrix} \bar{L} \\ \bar{M} \\ \bar{N} \end{pmatrix}
=
\begin{pmatrix}
\dfrac{1}{2}C_l \rho SV^2b \\
\dfrac{1}{2}C_m \rho SV^2 \bar{c} -I_R \omega_R･r \\
\dfrac{1}{2}C_n \rho SV^2b + I_R \omega_R･q
\end{pmatrix}
$なので、
航空機の回転運動の運動方程式
$
\begin{pmatrix}
\bar{L} \\
\bar{M} \\
\bar{N}
\end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p}-(I_y-I_z)qr-I_{xy}(rp-\dot{q})-I_{yz}(q^2-r^2)-I_{zx}(\dot{r}+pq) \\
I_y \dot{q}-(I_z-I_x)rp-I_{xy}(qr-\dot{p})-I_{yz}(\dot{r}+pq)-I_{zx}(r^2-p^2) \\
I_z \dot{r}-(I_x-I_y)pq-I_{xy}(p^2-q^2)-I_{yz}(\dot{q}+rp)-I_{zx}(\dot{p}-qr)
\end{pmatrix}
$
は、
$
\begin{pmatrix}
\dfrac{1}{2}C_l \rho SV^2b \\
\dfrac{1}{2}C_m \rho SV^2 \bar{c} -I_R \omega_R･r \\
\dfrac{1}{2}C_n \rho SV^2b + I_R \omega_R･q
\end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p}-(I_y-I_z)qr-I_{xy}(rp-\dot{q})-I_{yz}(q^2-r^2)-I_{zx}(\dot{r}+pq) \\
I_y \dot{q}-(I_z-I_x)rp-I_{xy}(qr-\dot{p})-I_{yz}(\dot{r}+pq)-I_{zx}(r^2-p^2) \\
I_z \dot{r}-(I_x-I_y)pq-I_{xy}(p^2-q^2)-I_{yz}(\dot{q}+rp)-I_{zx}(\dot{p}-qr)
\end{pmatrix}
$
機体が左右対称であれば、$I_{xy}=I_{yz}=0$ なので、
$
\begin{pmatrix}
\dfrac{1}{2}C_l \rho SV^2b \\
\dfrac{1}{2}C_m \rho SV^2 \bar{c} -I_R \omega_R･r \\
\dfrac{1}{2}C_n \rho SV^2b + I_R \omega_R･q
\end{pmatrix}
=
\begin{pmatrix}
I_x \dot{p} + (I_y-I_z) qr-I_{zx} ( \dot{r}+pq) \\
I_y \dot{q} + (I_z-I_x) rp-I_{zx} (r^2-p^2) \\
I_z \dot{r} + (I_x-I_y) pq_{zx} ( \dot{p}-qr)
\end{pmatrix}
$
と書くことができ、これを整理すると
$
\begin{pmatrix} I_x \dot{p} \\ I_y \dot{q} \\ I_z \dot{r} \end{pmatrix}
=
\begin{pmatrix}
(I_y-I_z)qr+I_{zx}(\dot{r}+pq)+\dfrac{1}{2}C_l \rho SV^2b \\
(I_z-I_x)rp+I_{zx}(r^2-p^2)+\dfrac{1}{2}C_m \rho SV^2 \bar{c} -I_R \omega_R･r \\
(I_x-I_t)pq+I_{zx}(\dot{p}-qr)+\dfrac{1}{2}C_n \rho SV^2b + I_R \omega_R･q
\end{pmatrix}
$
ここで、
$
\begin{pmatrix}
(I_y-I_z)qr+I_{zx}pq+\dfrac{1}{2}C_l \rho SV^2b \\
(I_z-I_x)rp+I_{zx}(r^2-p^2)+\dfrac{1}{2}C_m \rho SV^2 \bar{c} -I_R \omega_R･r \\
(I_x-I_t)pq-I_{zx}qr+\dfrac{1}{2}C_n \rho SV^2b + I_R \omega_R･q
\end{pmatrix}
=
\begin{pmatrix} L \\ M \\ N \end{pmatrix}
$
とおくと、
$
\begin{pmatrix} I_x \dot{p} \\ I_y \dot{q} \\ I_z \dot{r} \end{pmatrix}
=
\begin{pmatrix} L+I_{zx} \dot{r} \\ M \\ N+I_{zx} \dot{p} \end{pmatrix}
$
さらに整理すると、
$
\begin{pmatrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{pmatrix}
=
\begin{pmatrix}
\dfrac{\dfrac{L}{I_x}+\dfrac{I_{zx}}{I_x}･\dfrac{N}{I_z}}{1-\dfrac{{I_{zx}}^2}{I_xI_z}} \\ \\
\dfrac{M}{I_y} \\ \\
\dfrac{\dfrac{N}{I_z}+\dfrac{I_{zx}}{I_z}･\dfrac{L}{I_x}}{1-\dfrac{{I_{zx}}^2}{I_xI_z}}
\end{pmatrix}
$
これも（回転運動方程式）