[心得] 星體運動軌跡動畫之模擬設計

看板Mathematica作者harry901 (forcing to A cup)時間13年前 (2011/12/16 00:49)推噓1(1推 0噓 0→)

留言1則, 1人參與討論串1/2 (看更多)

今天突發奇想不知道mathematica能不能拿來設計遊戲... 可惜畢竟是數學軟體... 這文章的目的是寫出地球對於一天體之運行軌跡的影響，以動畫方式呈現。 ===========================天文物理背景知識================================== (1)一個質點受到一個向某定點方向直線的力量而作運動，稱為中心力運動(central- force motion)。而中心力通常是由於電場或重力場所引起。 (為了方便，我們以2D座標來計算) 為了考慮這個質點P，m質量，以及受力F，我們必須定義一個好用的座標系統，所以使用r，θ極座標來表示。其運動狀態如下圖: ↖切線方向軌道 ╲ ￣￣￣/╲ / ╲◎質點P /dθ↙﹨ / ／ F \ / ／θ 其中質點到定點之距離定義為r ⊙定點----------------------------- Σ Fr=-F=m*ar where ar=r方向之加速度 =m*[r''-r(θ')^2] ...eq(1) 說明:Fr代表整個系統在此時，r方向之受力，由於質點與定點的關係，所以定義r的方向為往外為'正'，故這時中心力(向心力)=-F=m*ar Σ F_θ=0=m*aθ=m(rθ''+2r'θ') ...eq(2) (2)這時我們可以得到兩個式子eq1,2 -F==m*[r''-r(θ')^2] 0==m(rθ''+2r'θ') 為了方便解聯立微分方程，必須將簡化之微分式做繁式，也就是將時間條件還原可以得到新的兩條式子 ╓ (d^2)r dθ ╕ -F==m║ -------- - r(-----)^2 ║=======>(a)式 ╙ dt^2 dt ╜ ╓ (d^2)θ dr dθ ╕ 0==m║ r------- + 2---- ---- ║=======>(b)式 ╙ dt^2 dt dt ╜ 由於(b)式中m=\=0，所以裡面那一陀要==0 ╓ (d^2)θ dr dθ ╕ 1 ╓ (d^2)θ 2r*dr dθ ╕ ║ r------- + 2---- ---- ║=0= ---║r^2------- + ---- ---- ║ ╙ dt^2 dt dt ╜ r ╙ dt^2 dt dt ╜ ^^^^^^^^^^^^^^^^^^^^^^^^^ f*g' + f'g =(fg)' 可見得f=r^2，g=dθ/dt 1 ╓ ╕ 1 ╓ dθ ╕' 1 d ╓ dθ ╕ 所以(b)式== ---║(fg)' ║ == ---║ r^2*----- ║ ==--- ----║ r^2*----- ║==0 r ╙ ╜ r ╙ dt ╜ r dt ╙ dt ╜ 此時，由(b)式最後衍生出來的微分式將其積分可得r^2*(dθ/dt)=h where h is a constant. 由於dA=0.5r^2*dθ 將上式除以時間因素可得 dA/dt=0.5r^2*(dθ/dt)=0.5h this is Kepler's law in the formulation. (3) 最後再將所有定義過的東西統合，再把速度條件v，位置條件r，初始速度條件v0 初始位置條件r0加入可以得到其軌道方程式 1/r=1/r0*[1-GMe/r0/v0^2]*cosθ+GMe/r0^2/v0^2 離心率e=1/r0*[1-GMe/r0/v0^2]*h^2/GMe 若e==0 質點運動==圓 e==1 拋物線 e<1 橢圓 e>1 雙曲線 ========================================================================== 我們來寫一個描述天體運行軌跡的程式吧以下程式基於： 1. 地球為引力主體，忽略其他天體 2. 地心為座標原點，使用卡式座標 3. 計算對象為對於任何天體(如隕石，衛星)之質量均可忽略 4. 初始條件有 a.天體對地心之速度r' b.天體對地心之角速度theta' c.天體之初始位置(r,theta) 基於以上我們在程式裡面使用數值方式求解利用Kepler's 2nd law 關於Kepler's law, refer to http://tinyurl.com/cloyesc ===========================Code Start===================================== endtime = 365*1;(*設定對象天體運行天數*) IniDis = 90000;(*初始距離地球km*) IniAngVel = 0.5;(*設定初始角速度,單位為每天相對於地球地心的徑度量*) IniRadVel = 1;(*設定對地心之徑向速度,單位為km/天*) M = 5.97*10^24;(*地球質量kg*) G = 6.67*10^-11;(*萬有引力常數m^3kg^-1s^-2*) GM = G*M*10^-9*(1/86400)^2;(*改成以天為單位,距離為km*) KEPeq1 := r''[t] - r'[t]*(theta'[t])^2 == -GM/(r[t])^2;(*Kepler D.E.EQ*) KEPeq2 := r[t]*theta''[t] + 2*r'[t]*theta'[t] == 0; NDSolve[{KEPeq1, KEPeq2, r[0] == IniDis, theta[0] == 0, r'[0] == IniRadVel, theta'[0] == IniAngVel}, {r[t], theta[t]}, {t, 0, endtime}]; tempr = Part[Part[%, 1], 1]; temptheta = Part[Part[%, 1], 2]; earth = Graphics[Disk[{0, 0}, 6370], PlotRange -> All, Axes -> True];(*標示地球,地球半徑6370km*) Animate[Show[{earth, ParametricPlot[{(r[t] /. tempr)* Sin[(theta[t] /. temptheta)], (r[t] /. tempr)* Cos[(theta[t] /. temptheta)]}, {t, 0, time}, Axes -> True, PlotRange -> All]}], {time, 0.1, endtime}, AnimationRate -> endtime/(100*IniAngVel), AnimationRunning -> False](*畫動畫*) ===================================End Code============================= 以上，改變IniAngVel(建議數值範圍0~5,負的只是方向相反) IniRadVel(0~100000) 會得到不同軌跡... 程式之修改，有興趣的人可以玩玩... 0. 可以利用軌跡方程式直接模擬天體運行軌跡 1. 可以作太空垃圾之模擬，產生n個太空垃圾，試圖以亂數產生其初始條件等等然後計算被地球擄獲的機率 2. 可以將地球質量等參數改為太陽，在不同引力下對天體運行的影響 3. 可以模擬八大行星運行的軌跡，並且將二維座標以三維方式呈現 4. 可以利用差分方式求得軌跡，這會更精確... 最後...其實有更好的方式來寫...只是每次碰到那個 /.或% 就很討厭=_= 如果版友發現程式內有什麼好玩的或是程式可以修改得更漂亮的話記得上來分享 ※ 編輯: harry901 來自: 122.120.6.78 (12/16 01:25)