高云峰

(清華大學(xué)航天航空學(xué)院， 北京 100084)

本系列文章介紹了很多利用數(shù)值計(jì)算處理問題的方法和過程。由于存在計(jì)算誤差，以及人為的輸入錯誤，計(jì)算得到的大量數(shù)據(jù)可能有誤。作為Matlab求解理論力學(xué)的系列文章，十分強(qiáng)調(diào)計(jì)算的可靠性，所以前面曾經(jīng)介紹過利用系統(tǒng)的守恒量來驗(yàn)證數(shù)值計(jì)算的精度，或者利用解析解進(jìn)行對比測試。如果系統(tǒng)不存在守恒量，可以考慮退化等特殊情況(這時往往有解析解或守恒量)，對退化情況檢驗(yàn)之后，再進(jìn)行后續(xù)的計(jì)算。

某些問題可能在不同的坐標(biāo)系中處理更方便，那如何判斷其是否正確？本文介紹一種數(shù)據(jù)轉(zhuǎn)換的方法，可以把不同坐標(biāo)系的結(jié)果轉(zhuǎn)換到希望的坐標(biāo)系中來觀察，并驗(yàn)證其準(zhǔn)確性。

為此本文設(shè)計(jì)了這樣兩個問題：(1)在非慣性系中看，隨手拋出的一個物體，其軌跡很復(fù)雜，我們?nèi)绾蝸砼袛嘤?jì)算是正確的呢？(2)在地球某點(diǎn)觀察太陽，太陽的視運(yùn)動軌跡是什么曲線？

1 不同坐標(biāo)系中的數(shù)據(jù)轉(zhuǎn)換關(guān)系

設(shè)有矢量ρ和兩個坐標(biāo)系Oxiyizi和Oxjyjzj(圖1)，基矢量分別為ei?(exi eyi ezi)T和ej?(exj eyj ezj)T。若用(ρ)i表示ρ在坐標(biāo)系Oxiyizi中的分量，有

圖1 矢量和坐標(biāo)系

利用基矢量，矢量ρ可以表示為

矢量與坐標(biāo)系沒有關(guān)系，但是它在不同坐標(biāo)系中的分量就與坐標(biāo)系有關(guān)系了。因此一個問題就是：(ρ)i與(ρ)j有什么關(guān)系？

將式(2) 改寫，有

由于兩個單位矢量的點(diǎn)積等于其夾角的余弦，所以Aij也稱為方向余弦矩陣，從而有[1]

從式(5) 類似有(ρ)i=Aik(ρ)k，從而得到坐標(biāo)系轉(zhuǎn)換矩陣的傳遞關(guān)系

以式(5)和式(6)為基礎(chǔ)可以得到更一般的關(guān)系。假設(shè)Ox1y1z1是基準(zhǔn)參考坐標(biāo)系(如動力學(xué)中的慣性系)，ox′y′z′是平動坐標(biāo)系(坐標(biāo)系始終與基準(zhǔn)坐標(biāo)系平行)，ox2y2z2是任意的動系(如動力學(xué)中的非慣性系)，ox2y2z2初始時刻與ox′y′z′重合，轉(zhuǎn)動后用三個角度(如歐拉角ψ，θ，φ) 來表示；P是動點(diǎn)，在基準(zhǔn)坐標(biāo)系中矢徑為R，相對運(yùn)動的矢徑為ρ，動系原點(diǎn)的矢徑為RO(圖2)。

圖2 空間一般坐標(biāo)系

根據(jù)運(yùn)動學(xué)關(guān)系，始終有[2]

式(7) 可以在任意坐標(biāo)系中表示，根據(jù)式(5)，則有

式(7) 給出了一般情況下坐標(biāo)轉(zhuǎn)換關(guān)系，其中坐標(biāo)轉(zhuǎn)換矩陣還是比較復(fù)雜的，如果用歐拉角表示，利用式(6) 的關(guān)系后，有

2 不同坐標(biāo)系中的拋物問題

兩位觀察者，A在地面上(慣性坐標(biāo)系)，B在勻速轉(zhuǎn)動的轉(zhuǎn)盤邊緣(非慣性系)。B隨手拋出一物體，求兩位觀察者看到的物體運(yùn)動軌跡(圖3，不考慮空氣阻力)。

圖3 拋物示意圖

圖 4 是運(yùn)動學(xué)關(guān)系示意圖， 假設(shè) (R)1=(a b c)T，轉(zhuǎn)盤半徑為r，轉(zhuǎn)角為φ=ωt，B的位置在動系中為(r)2=(0r0)T，出手初速度在動系中為(v)2=(vx0vy0vz0)T。

圖4 運(yùn)動學(xué)關(guān)系

注意物體在不同坐標(biāo)系中運(yùn)動的動力學(xué)方程是不同的，根據(jù)理論力學(xué)知識，在非慣性系中動力學(xué)方程為

展開后的方程及初始條件為

給予具體參數(shù)

可以計(jì)算出結(jié)果，如圖5 (空間三視圖) 和圖6 (xy平面俯視圖)。

圖5 曲線的空間視圖

圖6 曲線的俯視圖

看起來曲線很復(fù)雜，我們怎么判斷計(jì)算結(jié)果可靠呢？

眾所周知，在慣性系中拋出一個物體其軌跡是拋物線，這是很明確的結(jié)論。因此可以這樣處理：設(shè)慣性系中算出曲線為κ1，非慣性系中算出的曲線為κ2，把慣性系中的曲線κ1數(shù)據(jù)轉(zhuǎn)換到非慣性系中得到曲線κ3，如果κ3能與κ2完全重合(兩者誤差小于給定的誤差，如積分的允許誤差)，就認(rèn)為非慣性系中的結(jié)果是可靠的。

從慣性系的角度看，動力學(xué)方程有

展開后的動力學(xué)方程及初始條件為

積分后結(jié)果如圖7 和圖8。我們有理由認(rèn)為這些結(jié)果是合理可靠的，如曲線看起來是拋物線，俯視圖是一根直線，表示拋物線在一個平面內(nèi)。當(dāng)然還可以利用運(yùn)動過程中的機(jī)械能守恒來判斷(略)，總之我們認(rèn)為慣性系中這些曲線合理可靠。

圖7 曲線的空間視圖

圖8 曲線的俯視圖

現(xiàn)在的問題是如何利用慣性系的結(jié)果去證明非慣性系的結(jié)果正確？

式(7) 和式(8) 給出了不同坐標(biāo)系之間數(shù)據(jù)的轉(zhuǎn)換關(guān)系。具體到本問題中，轉(zhuǎn)換矩陣很簡單，z軸是平行的，退化為

如果想把慣性系的結(jié)果轉(zhuǎn)換到非慣性系，式(7)具體為

在Matlab 中，數(shù)據(jù)轉(zhuǎn)換的源代碼如圖9 所示。

圖9 數(shù)據(jù)轉(zhuǎn)換的源代碼

其中 (x1(i)y1(i)z1(i))T表示慣性系中第i個點(diǎn)，(x3(i)y3(i)z3(i))T是轉(zhuǎn)換后得到的非慣性系中第i個點(diǎn)。在Matlab 中，A12’ 表示對A12 進(jìn)行轉(zhuǎn)置。

數(shù)據(jù)轉(zhuǎn)換后，可以看到每個圓圈(非慣性系中直接算出的結(jié)果)都包含一個圓點(diǎn)(從慣性系中轉(zhuǎn)換得到的結(jié)果)，說明兩者的精度都很高(如圖10 和圖11所示)。當(dāng)然在比較的時候要注意，不同坐標(biāo)系中的積分要同步，用等步長積分正好可以實(shí)現(xiàn)這一點(diǎn)。

圖10 非慣性系中的空間圖

如果想把非慣性系的結(jié)果轉(zhuǎn)換到慣性系，結(jié)論依然很可靠(具體略)。綜合起來，說明圖5 和圖6的結(jié)果雖然有點(diǎn)出乎意料，卻是正確的。

根據(jù)這一事實(shí)，未來處理非慣性系的動力學(xué)問題時，除了直接處理非慣性系中的問題(動力學(xué)方程復(fù)雜，結(jié)果也復(fù)雜) 之外，還多了一種選擇：先在慣性系中進(jìn)行分析(動力學(xué)方程簡單，結(jié)論正確性容易驗(yàn)證)，然后再利用坐標(biāo)系之間的轉(zhuǎn)換關(guān)系，得到非慣性系中相應(yīng)的結(jié)果。

3 太陽的視運(yùn)動問題

人人都知道地球繞太陽的運(yùn)動很簡單，但是如果做個調(diào)查，你在當(dāng)?shù)乜吹降奶栠\(yùn)動軌跡(視運(yùn)動軌跡) 是什么曲線？有什么特點(diǎn)？估計(jì)絕大部分人一時答不上來。而利用數(shù)據(jù)轉(zhuǎn)換關(guān)系可以得出既簡單又出乎意料的結(jié)論。

假設(shè)地球繞太陽為圓周運(yùn)動，為簡單設(shè)365 天。在日心坐標(biāo)系OXY Z中，地軸(z軸) 的單位矢量為(0，-sinθ，cosθ)T，θ=23.5°(圖12)。設(shè)時間為春分點(diǎn)之后第t小時，地球球心運(yùn)動至

其中AU 為天文單位。設(shè)跟隨地心運(yùn)動的地心坐標(biāo)系為oxyz，則轉(zhuǎn)換矩陣為

下標(biāo)S 表示太陽坐標(biāo)系，E 表示地球坐標(biāo)系。

在地球表面P點(diǎn)，建立當(dāng)?shù)氐臇|北天坐標(biāo)系PENZ(圖12)。在地心坐標(biāo)系中，P點(diǎn)的表達(dá)式為其中Re為地球半徑，φ為當(dāng)?shù)鼐暥龋?2πt/24 為地球自轉(zhuǎn)角度。

東北天坐標(biāo)系相對地心坐標(biāo)系的轉(zhuǎn)換矩陣為

下標(biāo)G 表示地理坐標(biāo)系。

綜合圖12 和圖13，仍有運(yùn)動學(xué)關(guān)系R=Ro+ρ，在地理坐標(biāo)系中分解為

圖12 地球相對太陽的運(yùn)動

圖13 當(dāng)?shù)貣|北天坐標(biāo)系

根據(jù)式(6)、式(17) 和式(19)，有

我們看到的太陽視運(yùn)動， 就是式 (21) 中的-(R)G。式(21) 右邊每項(xiàng)都為已知函數(shù)，代入數(shù)據(jù)后太陽視運(yùn)動軌跡見圖14，顯示了全年特殊日期的太陽視運(yùn)動軌跡：夏至當(dāng)天軌跡偏北方，冬至當(dāng)天軌跡偏南方，其他時間視運(yùn)動軌跡會在夏至和冬至軌跡之間緩慢地變化，春分當(dāng)天軌跡居中間。我們看不到水平面下方的太陽視運(yùn)動軌跡，但仍畫出。注意每一天的視運(yùn)動軌跡并不封閉，因此全年的太陽視運(yùn)動軌跡則很像一個彈簧，但是這個彈簧兩端半徑微小，中間半徑略大，不過差距沒有圖中彈簧那么明顯(圖15)，另外視運(yùn)動軌跡很密集，在夏至與冬至之間有365 圈。

圖14 太陽某一天的視運(yùn)動

圖15 全年太陽視運(yùn)動示意圖

4 小結(jié)

本文介紹了不同坐標(biāo)系之間數(shù)據(jù)轉(zhuǎn)換的一般方法，利用這一方法，可以把問題先在某個坐標(biāo)系中處理好，然后再轉(zhuǎn)換到另一個坐標(biāo)系中。

在非慣性系拋體問題中，直接處理的方程比較復(fù)雜，結(jié)果也很復(fù)雜，難以判斷其正確性。但是在慣性系中，動力學(xué)方程很簡單，結(jié)果的物理含義也很清楚。利用不同坐標(biāo)系之間的數(shù)據(jù)轉(zhuǎn)換關(guān)系，可以很容易證明兩種動力學(xué)方程的計(jì)算結(jié)果是正確的。

而在太陽視運(yùn)動的問題中，從慣性系看地球的運(yùn)動很簡單，但在地球上看太陽的運(yùn)動，則比較復(fù)雜，出乎很多人的意料。利用數(shù)據(jù)轉(zhuǎn)換關(guān)系，可以把幾個簡單的運(yùn)動(地球繞太陽的運(yùn)動是圓周運(yùn)動、地球上某點(diǎn)相對地球的運(yùn)動也很簡單) 直接轉(zhuǎn)換為需要的結(jié)果。