本系列文章介绍了很多利用数值计算处理问题的方法和过程。由于存在计算误差,以及人为的输入错误,计算得到的大量数据可能有误。作为Matlab求解理论力学的系列文章,十分强调计算的可靠性,所以前面曾经介绍过利用系统的守恒量来验证数值计算的精度,或者利用解析解进行对比测试。如果系统不存在守恒量,可以考虑退化等特殊情况(这时往往有解析解或守恒量),对退化情况检验之后,再进行后续的计算。
某些问题可能在不同的坐标系中处理更方便,那如何判断其是否正确?本文介绍一种数据转换的方法,可以把不同坐标系的结果转换到希望的坐标系中来观察,并验证其准确性。
为此本文设计了这样两个问题：(1)在非惯性系中看,随手抛出的一个物体,其轨迹很复杂,我们如何来判断计算是正确的呢?(2)在地球某点观察太阳,太阳的视运动轨迹是什么曲线?
1 不同坐标系中的数据转换关系
设有矢量${\mathbf \rho }$和两个坐标系$Ox_{i} y_{i} z_{i} $和 $Ox_{j} y_{j} z_{j}$ (图1),基矢量分别为${\mathbf e}_{i}\triangleq \left({\mathbf e}_{xi}\ \ {\mathbf e}_{yi}\ \ {\mathbf e}_{zi}\right)^{\rm T}$和${\mathbf e}_{j}\triangleq \left({\mathbf e}_{xj}\ \ {\mathbf e}_{yj}\ \ {\mathbf e}_{zj}\right)^{\rm T}$。若用$(\mathbf \rho )_{i} $表示${\mathbf \rho}$在坐标系$Ox_{i} y_{i} z_{i} $中的分量,有
$\begin{eqnarray} (\mathbf \rho )_{i} =\left(x_{i}\ \ y_{i}\ \ z_{i}\right)^{\rm T}, \ \ (\mathbf \rho )_{j} =\left(x_{j}\ \ y_{j}\ \ z_{j}\right)^{\rm T} \end{eqnarray}$
利用基矢量,矢量${\mathbf \rho }$可以表示为
$\begin{eqnarray} {\mathbf \rho }={\mathbf e}_{i}^{\rm T} (\mathbf \rho )_{i} ={\mathbf e}_{j}^{\rm T} (\mathbf \rho)_{j} \end{eqnarray}$
矢量与坐标系没有关系,但是它在不同坐标系中的分量就与坐标系有关系了。因此一个问题就是：$(\mathbf \rho)_{i} $与$(\mathbf \rho )_{j} $有什么关系?
图1
$\begin{eqnarray} (\rho )_{i} ={\mathbf e}_{i} \cdot {\mathbf e}_{j}^{\rm T} (\rho )_{j} \end{eqnarray}$
定义坐标转换矩阵$\mathbf A_{ij} $为
$\begin{eqnarray} \mathbf A_{ij}\triangleq {\mathbf e}_{i} \cdot {\mathbf e}_{j}^{\rm T} =\left[ {{\begin{array}{*{20}c} {{\mathbf e}_{xi} \cdot {\mathbf e}_{xj} } & {{\mathbf e}_{xi} \cdot {\mathbf e}_{yj} } & {{\mathbf e}_{xi} \cdot {\mathbf e}_{zj} } \\[-1mm] {{\mathbf e}_{yi} \cdot {\mathbf e}_{xj} } & {{\mathbf e}_{yi} \cdot {\mathbf e}_{yj} } & {{\mathbf e}_{yi} \cdot {\mathbf e}_{zj} } \\[-1mm] {{\mathbf e}_{zi} \cdot {\mathbf e}_{xj} } & {{\mathbf e}_{zi} \cdot {\mathbf e}_{yj} } & {{\mathbf e}_{zi} \cdot {\mathbf e}_{zj} } \\[-1mm] \end{array} }} \right] \end{eqnarray}$
由于两个单位矢量的点积等于其夹角的余弦,所以$\mathbf A_{ij}$也称为方向余弦矩阵,从而有[1]
$\begin{eqnarray} (\mathbf \rho )_{i} =\mathbf A_{ij} (\mathbf \rho )_{j} \end{eqnarray}$
从式(5)类似有$(\mathbf \rho )_{i} =\mathbf A_{ik} (\mathbf \rho )_{k}$,从而得到坐标系转换矩阵的传递关系
$\begin{eqnarray} \mathbf A_{ik} =\mathbf A_{ij} \mathbf A_{jk} \end{eqnarray}$
以式(5)和式(6)为基础可以得到更一般的关系。假设$Ox_{1} y_{1} z_{1}$是基准参考坐标系(如动力学中的惯性系),$ox'y'z'$是平动坐标系(坐标系始终与基准坐标系平行),$ox_{2}y_{2} z_{2} $是任意的动系(如动力学中的非惯性系),$ox_{2} y_{2} z_{2}$初始时刻与$ox'y'z'$重合,转动后用三个角度(如欧拉角$\psi,\theta,\varphi$)来表示;$P$是动点,在基准坐标系中矢径为${\mathbf R}$,相对运动的矢径为${\mathbf \rho }$,动系原点的矢径为${\mathbf R}_{O}$ (图2)。
图2
$\begin{eqnarray} {\mathbf R}={\mathbf R}_{O} +{\mathbf \rho } \end{eqnarray}$
式(7)可以在任意坐标系中表示,根据式(5),则有
$\begin{eqnarray} (\mathbf R)_{1} =(\mathbf R_{O} )_{1} +\mathbf A_{12} (\mathbf \rho )_{2} \end{eqnarray}$
式(7)给出了一般情况下坐标转换关系,其中坐标转换矩阵还是比较复杂的,如果用欧拉角表示,利用式(6)的关系后,有
$\begin{eqnarray} &&\mathbf A_{12} =\left( {{\begin{array}{*{20}c} {\cos \psi } & {-\sin \psi } & 0 \\[-1mm] {\sin \psi } & {\cos \psi } & 0 \\[-1mm] 0 & 0 & 1 \\[-1mm] \end{array} }} \right)\left( {{\begin{array}{*{20}c} 1 & 0 & 0 \\[-1mm] 0 & {\cos \theta } & {-\sin \theta } \\[-1mm] 0 & {\sin \theta } & {\cos \theta } \\ \end{array} }} \right)\cdot \\&&\qquad \left( {{\begin{array}{*{20}c} {\cos \varphi } & {-\sin \varphi } & 0 \\[-1mm] {\sin \varphi } & {\cos \varphi } & 0 \\[-1mm] 0 & 0 & 1 \\ \end{array} }} \right) \end{eqnarray}$
2 不同坐标系中的抛物问题
两位观察者,$A$在地面上(惯性坐标系),$B$在匀速转动的转盘边缘(非惯性系)。$B$随手抛出一物体,求两位观察者看到的物体运动轨迹(图3,不考虑空气阻力)。
图3
图4是运动学关系示意图,假设$(\mathbf R)_{1} =( a\ \ b\ \ c)^{\rm T}$,转盘半径为$r$,转角为$\varphi =\omega t$,$B$的位置在动系中为$(\mathbf r)_{2}=(0\ \ r\ \ 0)^{\rm T}$,出手初速度在动系中为$(\mathbf v)_{2}=\left({v_{x0} } \ \ {v_{y0} } \ \ {v_{z0} } \right)^{\rm T}$。
图4
注意物体在不同坐标系中运动的动力学方程是不同的,根据理论力学知识,在非惯性系中动力学方程为
$\begin{eqnarray} m{\tilde{{\ddot{{\mathbf r}}}}}=m{\mathbf g}-m\mathbf \omega \times (\mathbf \omega \times \mathbf r)-2m\mathbf \omega \times \tilde{{\dot{{\mathbf r}}}} \end{eqnarray}$
$\begin{eqnarray} \left. {\begin{array}{l} \tilde{{\ddot{{x}}}}=\omega^{2}x+2\omega \tilde{{\dot{{y}}}}, \\ \tilde{{\ddot{{y}}}}=\omega^{2}y-2\omega \tilde{{\dot{{x}}}}, \\ \tilde{{\ddot{{z}}}}=-g, \\ \end{array}} \right. \ \ \left. {\begin{array}{l} x=0, \\ y=r, \\ z=0, \\ \end{array}} \right.\ \ \left. {\begin{array}{l} \tilde{{\dot{{x}}}}_{0} =v_{x0} \\ \tilde{{\dot{{y}}}}_{0} =v_{y0} \\ \tilde{{\dot{{z}}}}_{0} =v_{z0} \\ \end{array}} \right. \end{eqnarray}$
$ \begin{eqnarray*} &&m=1 {\rm kg},\ g=9.8 {\rm m/s}^{2},\ \omega =10 {\rm rad/s}\\[-2mm]&& r=1 {\rm m} a=1 {\rm m},\ b=2 {\rm m},\ c=3 {\rm m},\ v_{x0} =0\\[-2mm]&&v_{y0} =-1 {\rm m/s},\ v_{z0} =30 {\rm m/s} \end{eqnarray*}$
可以计算出结果,如图5 (空间三视图)和图6 ($xy$平面俯视图)。
图5
图6
众所周知,在惯性系中抛出一个物体其轨迹是抛物线,这是很明确的结论。因此可以这样处理：设惯性系中算出曲线为$\kappa_{1} $,非惯性系中算出的曲线为$\kappa_{2} $,把惯性系中的曲线$\kappa_{1}$数据转换到非惯性系中得到曲线$\kappa_{3} $,如果$\kappa_{3} $能与$\kappa_{2}$完全重合(两者误差小于给定的误差,如积分的允许误差),就认为非惯性系中的结果是可靠的。
$\begin{eqnarray} m{\ddot{{\mathbf r}}}=m{\mathbf g} \end{eqnarray}$
$\begin{eqnarray} \left. {\begin{array}{l} m\ddot{{x}}=0, \\ m\ddot{{y}}=0, \\ m\ddot{{z}}=-mg, \\ \end{array}} \right.\ \ \left. {\begin{array}{l} x=a, \\ y=b+r, \\ z=c, \\ \end{array}} \right. \ \ \left. {\begin{array}{l} \dot{{x}}=v_{x0} -r\omega \\ \dot{{y}}=v_{y0} \\ \dot{{z}}=v_{z0} \\ \end{array}} \right. \end{eqnarray}$
积分后结果如图7和图8。我们有理由认为这些结果是合理可靠的,如曲线看起来是抛物线,俯视图是一根直线,表示抛物线在一个平面内。当然还可以利用运动过程中的机械能守恒来判断(略),总之我们认为惯性系中这些曲线合理可靠。
图7
图8
现在的问题是如何利用惯性系的结果去证明非惯性系的结果正确?
式(7)和式(8)给出了不同坐标系之间数据的转换关系。具体到本问题中,转换矩阵很简单,$z$轴是平行的,退化为
$\begin{eqnarray} \mathbf A_{12} =\left( {{\begin{array}{*{20}c} {\cos \varphi } & {-\sin \varphi } & 0 \\ {\sin \varphi } & {\cos \varphi } & 0 \\ 0 & 0 & 1 \\ \end{array} }} \right) \end{eqnarray}$
如果想把惯性系的结果转换到非惯性系,式(7)具体为
$\begin{eqnarray} ({\rho_{x} } \ \ {\rho_{y} } \ \ {\rho_{z} } )_{2}^{\rm T}=\mathbf A_{12}^{\rm T}( {R_{x} -a} \ \ {R_{y} -b} \ \ {R_{z} }-c)_{0}^{\rm T} \end{eqnarray}$
图9
其中$\left({x_{1} (i)} \ \ {y_{1} (i)} \ \ {z_{1} (i)} \right)^{\rm T}$表示惯性系中第$i$个点,$\left({x_{3} (i)} \ \ {y_{3} (i)} \ \ {z_{3} (i)}\right)^{\rm T}$是转换后得到的非惯性系中第$i$个点。在Matlab中,A12'表示对A12进行转置。
数据转换后,可以看到每个圆圈(非惯性系中直接算出的结果)都包含一个圆点(从惯性系中转换得到的结果),说明两者的精度都很高(如图10和图11所示)。当然在比较的时候要注意,不同坐标系中的积分要同步,用等步长积分正好可以实现这一点。
图10
图11
如果想把非惯性系的结果转换到惯性系,结论依然很可靠(具体略)。综合起来,说明图5和图6的结果虽然有点出乎意料,却是正确的。
根据这一事实,未来处理非惯性系的动力学问题时,除了直接处理非惯性系中的问题(动力学方程复杂,结果也复杂)之外,还多了一种选择：先在惯性系中进行分析(动力学方程简单,结论正确性容易验证),然后再利用坐标系之间的转换关系,得到非惯性系中相应的结果。
3 太阳的视运动问题
人人都知道地球绕太阳的运动很简单,但是如果做个调查,你在当地看到的太阳运动轨迹(视运动轨迹)是什么曲线?有什么特点?估计绝大部分人一时答不上来。而利用数据转换关系可以得出既简单又出乎意料的结论。
假设地球绕太阳为圆周运动,为简单设365天。在日心坐标系$OXYZ$中,地轴($z$轴)的单位矢量为$(0, -\sin \theta, \cos \theta)^{\rm T}$,$\theta=23.5^{\circ}$ (图12)。设时间为春分点之后第$t$小时,地球球心运动至
$\begin{eqnarray} \left.\begin{array}{l} \left( {R_{o} } \right)_{S} =\mbox{AU}\left( {\cos \alpha } \ \ {\sin \alpha } \ \ 0\right)^{\rm T}\\ \alpha =2\pi t/(365\ast 24) \end{array}\right\} \end{eqnarray}$
其中AU为天文单位。设跟随地心运动的地心坐标系为$oxyz$,则转换矩阵为
$\begin{eqnarray} \mathbf A_{\rm SE} =\left( {{\begin{array}{*{20}c} 1 & 0 & 0 \\ 0 & {\cos \theta } & {-\sin \theta } \\ 0 & {\sin \theta } & {\cos \theta } \\ \end{array} }} \right) \end{eqnarray}$
在地球表面$P$点,建立当地的东北天坐标系$PENZ$ (图12)。在地心坐标系中,$P$点的表达式为
$\begin{eqnarray} \left( \mathbf \rho \right)_{\rm E} =\left( {{\begin{array}{*{20}c} {R_{\rm e} \cos \varphi \cos \beta } \\ {R_{\rm e} \cos \varphi \sin \beta } \\ {R_{\rm e} \sin \varphi } \\ \end{array} }} \right) \end{eqnarray}$
其中$R_{\rm e} $为地球半径,$\varphi $为当地纬度,$\beta =2{\pi}t/24$为地球自转角度。
$\begin{eqnarray} \mathbf A_{\rm EG} =\left( {{\begin{array}{*{20}c} {-\sin \beta } & {-\sin \varphi \cos \beta } & {\cos \varphi \cos \beta } \\ {\cos \beta } & {-\sin \varphi \sin \beta } & {\cos \varphi \sin \beta } \\ 0 & {\cos \varphi } & {\sin \varphi } \\ \end{array} }} \right) \end{eqnarray}$
综合图12和图13,仍有运动学关系${\mathbf R}={\mathbf R}_{o} +{\mathbf \rho}$,在地理坐标系中分解为
$\begin{eqnarray} (\mathbf R)_{\rm G} =\mathbf A_{\rm GS} (\mathbf R_{o} )_{\rm S} +\mathbf A_{\rm GE} (\mathbf \rho )_{\rm E} \end{eqnarray}$
图12
图13
$\begin{eqnarray} (\mathbf R)_{\rm G} =\mathbf A_{\rm GE} \mathbf A_{\rm ES} (\mathbf R_{o} )_{\rm S} +\mathbf A_{\rm EG}^{\rm T}(\mathbf \rho )_{\rm E} \end{eqnarray}$
我们看到的太阳视运动,就是式(21)中的$-(\mathbf R)_{\rm G}$。式(21)右边每项都为已知函数,代入数据后太阳视运动轨迹见图14,显示了全年特殊日期的太阳视运动轨迹：夏至当天轨迹偏北方,冬至当天轨迹偏南方,其他时间视运动轨迹会在夏至和冬至轨迹之间缓慢地变化,春分当天轨迹居中间。我们看不到水平面下方的太阳视运动轨迹,但仍画出。注意每一天的视运动轨迹并不封闭,因此全年的太阳视运动轨迹则很像一个弹簧,但是这个弹簧两端半径微小,中间半径略大,不过差距没有图中弹簧那么明显(图15),另外视运动轨迹很密集,在夏至与冬至之间有365圈。
图14
图15
4 小结
本文介绍了不同坐标系之间数据转换的一般方法,利用这一方法,可以把问题先在某个坐标系中处理好,然后再转换到另一个坐标系中。
在非惯性系抛体问题中,直接处理的方程比较复杂,结果也很复杂,难以判断其正确性。但是在惯性系中,动力学方程很简单,结果的物理含义也很清楚。利用不同坐标系之间的数据转换关系,可以很容易证明两种动力学方程的计算结果是正确的。
而在太阳视运动的问题中,从惯性系看地球的运动很简单,但在地球上看太阳的运动,则比较复杂,出乎很多人的意料。利用数据转换关系,可以把几个简单的运动(地球绕太阳的运动是圆周运动、地球上某点相对地球的运动也很简单)直接转换为需要的结果。
参考文献
李俊峰, 张雄. 理论力学, 第2版. 北京: 清华大学出版社, 2010
[本文引用: 1]
1
1987
... 由于两个单位矢量的点积等于其夹角的余弦,所以$\mathbf A_{ij}$也称为方向余弦矩阵,从而有[1] ...
