高中研發人才培育學生研究計畫

利用網路上免費軟體進行天體運行的計算

作者: 黃茂誠, 許鈞筌

指導老師: 游大立

研究計劃題目：

利用網路上免費軟體進行天體運行的計算

壹、研究動機：

在偶然中看見陳威尹(第40屆全國科展地科作品第一名)的作品,想運用其原理深入研究天體長期運行的演化。在現在發達的網路世界中,有許多免費的自由軟體(freeware)供我們使用,那我們就不必花錢去買商業軟體來解決這些科學問題,同時善加運用網路的免費的自由軟體, 改寫原有的程式來進行研究。

貳、研究目的：

驗證在Linux作業系統下,用Fortran語言來執行與陳威尹(第40屆全國科展)在Windows98下用Visual Basic語言執行是否一致。並且運用改寫後程式, 深入研究天體長期運行的演化。

參、研究設備及器材
- 一、硬體

個人電腦(PC) 3部

二、軟體
1、Debian Linux 作業系統
2、G77

FORTRAN 77版, 簡易的科學程式

3、Gnuplot

科學繪圖軟體

肆、研究問題：
1. 雙星系統中行星的運行軌跡。
2. 三恆星系統的運動。
3. 行星系統中,小行星是否能穩定的存在。
4. 拉格朗日點附近小行星的運動軌跡。
5. 太陽系小行星長期軌道演化。
伍、研究步驟：
1. 將所需作業系統、軟體及相關套件安裝完畢。
2. 學習並熟悉所使用的Fortran語法。
3. 撰寫程式進行模擬運算。
4. 討論計算結果的意義。
5. 結果是否和陳威尹(第40屆全國科展地科作品第一名作者)執行一樣。
6. 將改寫的程式推廣至研究太陽系小行星長期軌道演化計算。
陸、研究方法：

利用一般個人電腦,搭配Linux作業系統,學習Fortran,並利用物理觀念撰寫程式,探討天體運動。

系統安裝

本研究使用Debian Linux作業系統.

科學計算的程式編譯器為G77

計算結果的繪圖程式為GNUPlot

成果展示在以Apache架設的網站中, 並利用PmWiki套件加以快速管理。

以上的軟體均為網路上可以自由且免費取得的軟體。

Linux的簡介

Linux的核心是由Linus Torvalds在1991年的時候發出來的，並且丟到網路上供人下載，後來大家覺得這個小東西不錯(Linus Kernel),相當的小而精巧，所以慢慢的就有相當多的人投入這東西的研究領域裡了。

Linux的核心是『Kernel』，這是作業系統最底層的東西，由其掌管整個硬體資源的工作狀態，因此它們各有不同核心。其實Linux就是一個Kernel！它提供了一個作業系統當中最底層的硬體控制與資源管理的完整架構，這個架構是沿襲Unix良好的傳統而來，所以相當的穩定而功能強大。

Linux的優點有：

穩定的系統：不易當機。
免費或少許費用：合於專題研究主題。
安全性、漏洞的修補：資料安全性高。
多工、多使用者：在研究時可以讓不同的使用者同時登入主機，同時使用主機的資源。
相對不耗資源的系統：只要p-100以上等級的電腦就可以安裝並使用。

G77的簡介

G77就是Fortran的編譯器,在現在的LINUX套件中都會包含GNU的Fortran編譯器G77,那Fortran是1954年IBM公司的一項計畫,它可以把使用接近的文字翻譯成機械語言,它是第一個被廣泛使用的高階語言,Fortran的名字來自於Formula Translator,這兩個字正是數學公式翻譯器的意思,Fortran的出現,使得電腦語言進入了高階語言的時代,後來在1966年由美國國家標準局的前身(American Standards Associaton)制定了Fortran的統一標準,視為Fortran66,後來又在1977年制定了新的Fortran語言標準,由美國國家標準局(American National Standards Institute)正式公佈,就是我們現在在用的Fortran77。

GNUPlot的簡介

gnuplot是一個命令導向的交談式跨平台數學程式繪圖程式。它會依照使用者輸入的每一項命令逐步修改或設定繪圖環境，後以圖形表達數據或函數，讓我們可藉由圖形做數據分析。它可以處理至少一百萬行的資料，這是大部分的繪圖軟體做不到的，而且這對於會產生大量數據的科學實驗幫助極大，尤其是天文方面的位置數據，如果繪出可以更容易顯現出週期性變化。它可產生多種開放性的圖形檔案格式的輸出，供編譯軟體使用。最重要的是它可以與其他軟體保持良好的溝通，例如：awk、sed、……等。所以在Shell命令稿中輸入了awk命令以及gnuplot的指令以加速圖檔的繪製。以下是例子：分割檔案至可讀之最大極限awk檔,包含繪圖指令的命令稿。

Apache/PmWiki套件

Apache 是很常見的網站伺服器軟體, 屬於自由軟體。目前常使用來與PHP網頁程式搭配使用。我們使用這樣子的搭配, 是為了能使用PmWiki網站管理程式, 以方便我們撰寫研究的報告資料。

PmWiki 是 PatrickMichaud 使用 PHP 程式語言開發出來的一套類似 WikiWikiWeb 的程式。 PmWiki 的主要目的是希望能夠讓多人同時且很輕鬆地製作和維護網站。 PmWiki 的散佈須受 the General Public License (GPL) 授權條款的限制。

在 PmWiki 裡建立的網頁並不需要特別的工具來編輯, 您所需要的只是一個標準的網站瀏覽器和一套很簡單的標示語言。透過這樣簡單的介面, 任何人都能輕輕鬆鬆地建立動態展示的網站, 群組軟體系統, 文件管理系統以及知識庫。 PmWiki 的特色包括 (這裡只列出一部份): 很容易安裝, 維護, 設定, 和昇級。

網頁可以在任何標準的網頁瀏覽器裡編輯 -- 只需按一下在任何網頁上的 "Edit Page" (編輯本頁)連結。簡單但很有用的文件格式規則在網頁之間自動連結, 包括 WikiWord 連結和 free links 連結。密碼-可保護網頁, 限制瀏覽和編輯組織網頁構成Wiki 群組的功能, 在群組內的網頁間, 共用標題和存取控管追蹤網頁更新的軌跡, 輕易地取回和回復任一網頁先前的版本 WikiTrails 可輕鬆地建立網頁間的循序連結檔案上傳和附件。

計算公式的推導:

F為萬有引力,A向心加速度,M,m分別為兩顆天體的質量,R為他們的距離,(X1,Y1,Z1),(X2,Y2,Z2)分別為他們的初始位置,時間間隔為 $\Large\Delta t$

由牛頓的第二運動定律和萬有引力定律可知:

$\Large\longrightarrow F=mA=\frac {GMm} {R^2}$

$\Large\longrightarrow A=\frac {GM} {R^2}$

那我們再把A拆解為x方向分量,Y方向分量,Z方向分量:

$\Large AX_0,AY_0,AZ_0$ 為初始的加速度

$\Large\longrightarrow AX=AX_0+A\time\frac {X1-X2} {R}$

$\Large\longrightarrow AY=AY_0+A\time\frac {Y1-Y2} {R}$

$\Large\longrightarrow AZ=AZ_0+A\time\frac {Z1-Z2} {R}$

那我們又可以由加速度算出後來的速度: $\Large\Delta V=A \time \Delta t$ ,所以後來的速度 $\Large V=$ 初始速度 $\Large +\Delta V$

$\Large VX_0,VY_0,VZ_0$ 為初始的速度

$\Large\longrightarrow VX=VX_0+AX\time{\Delta}t$

$\Large\longrightarrow VY=VY_0+AY\time{\Delta}t$

$\Large\longrightarrow VZ=VZ_0+AZ\time{\Delta}t$

我們又可以由速度推得後來的位置:位移 $\Large\Delta S=V\time\Delta t$ ,所以後來的位置 $\Large X=$ 初始位置 $\Large +\Delta S$

$\Large X_0,Y_0,Z_0$ 為初始的位置

$\Large\longrightarrow X=X_0+VX\time{\Delta}t$

$\Large\longrightarrow Y=Y_0+VY\time{\Delta}t$

$\Large\longrightarrow Z=Z_0+VZ\time{\Delta}t$

柒、研究結果
說明:在分析這些問題中,蝴蝶效應和重力彈弓,常是影響結果的因素。
蝴蝶效應:一些初始條件的些微不同,會導致後來的計算結果。
這兩張圖小行星的初始速度完全一樣,只有初始位置有些許的不同,最後導致結果完全不同。
這張圖小行星的位置是在恆星-行星連線夾 $\Large180^\circ$ 的地方

這張圖小行星的位置是在恆星-行星連線夾 $\Large180.0001^\circ$ 的地方

由這兩張圖可以看出蝴蝶效應對計算結果的影響,從上面兩張圖我們發現小行星的位置差距非常微小(角度改變量 $\Large\frac 1 {10000}$ ),但計算出來的結果完全不同。
重力彈弓:

在這個計算中,我們把小行星擺在和恆星-行星連線夾 $\Large350^\circ$ 的地方,在小行星靠近行星的地方,小行星受到行星力量的加速,所以小行星就被甩出去了。

計算程式碼

一、雙星系統中行星的運行軌跡

在這個雙星系統中,有兩顆質量相同的恆星和一顆質量微小的行星(故可以忽略)
在這個實驗中的G值(重力常數)=1
初始條件:雙星的質量相同,行星的質量為零。
固定變因:雙星的質量、位置和初速度
操縱變因:行星位置、初速度

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	可變	可變	可變	可變	可變	可變

在這個計算中由於行星的質量和恆星的質量差距甚大,所以行星的重力場對兩顆恆星的影響不大,故可以忽略不計,所以在程式中我們把行星的質量視為0,並非行星沒有質量。

在這個計算中,我們把 $\Large\Delta t$ 取為 $1\time10^{-2}$ ,計算時間到 $1\time10^{2}$
$\Large A=\frac {GM} {R^2}=\frac {V^2} R$
恆星 $\Large V=\sqrt {AR}=\sqrt {{\frac {GM} {(2R^2)}} R}= \frac 1 2$

其中的X代表行星的初始位置在(X,0,0),下面的四張圖為行星被恆星甩出去
Result of X=1.5

從這一張圖中,我們可以發現在越靠近雙星系統的行星,其行星軌道就越混亂,越不穩定。

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	1.5	0	0	0	$\Large\sqrt {\frac 2 {1.5}}$	0

Result of X=3

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	3	0	0	0	$\Large \sqrt {\frac 2 {3}}$	0

Result of X=-2

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	-2	0	0	0	$\Large \sqrt {\frac 2 {2}}$	0

下面的圖為"行星形成外行星"(繞著雙星的質星做運轉)
Result of X=5

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	5	0	0	0	$\Large \sqrt {\frac 2 {5}}$	0

Result of X=-1.5

下面的兩張圖為"行星形成內行星"(繞著其中一顆恆星做運轉)

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	-1.5	0	0	0	$\Large\sqrt {\frac 2 {1.5}}$	0

Result of X=-0.5

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	0.5	0
恆星2	1	-1	0	0	0	-0.5	0
行星	$\Large\longrightarrow 0$	-0.5	0	0	0	$\Large \sqrt {\frac 2 {0.5}}$	0

由上面這幾張圖可以了解:如果行星的初始位置和速度隨便給定的話,在雙星附近幾乎無法存在行星,也就是說會被甩出去,所以有則形成內行星或外行星兩種模式,在我們的計算時間內,它們的軌跡是穩定的狀況,但是我們並不能保證它們在計算時間外能夠一直很穩定的存在。

二、三恆星系統的運動

在這個系統中同時存在著三顆質量相同的天體(故都不可忽略)
在這個實驗中的G值(重力常數)=1

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	可變	可變	可變	可變	可變	可變
恆星2	1	可變	可變	可變	可變	可變	可變
恆星3	1	可變	可變	可變	可變	可變	可變

在這個計算中,我們把 $\Large\Delta$ t取為 $1\time10^{-2}$ 計算時間到 $1\time10^2$
恆星 $\Large V= \sqrt {\frac {G M} R}$

在這個計算中我們試著設計出一個巢狀雙星(其中的兩個恆星先形成小雙星,再和另一個恆星形成大雙星),但是因為在三星系統中的重力場是不穩定的,所以它們的軌跡並非一個平衡的狀態。
下面的這張圖是我們使用與陳威尹作品中的巢狀雙星的相同初始值,但是因為蝴蝶效應的影響,導致我們在模擬時,破壞了平衡,最後就都甩出去了。

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	-2.125	0	0	-5.23386590874456E-16	-1.03107632948017	0
恆星2	1	-3.125	0	0	-5.23386590874456E-16	0.383137232892924	0
恆星3	1	3.5	0	0	6.97848787832608E-16	0.431959397724831	0

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	2	0	0	0	1	0
恆星2	1	0.5	0	0	0	$\Large\sqrt {\frac 2 {0.5}}$	0
恆星3	1	-0.5	0	0	0	- $\Large\sqrt {\frac 2 {0.5}}$	0

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	2	0	0	0	1	0
恆星2	1	-1.5	0	0	0	$\Large\sqrt {\frac 2 {0.5}}$	0
恆星3	1	-2.5	0	0	0	- $\Large \sqrt {\frac 2 {0.5}}$	0

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	1	0
恆星2	1	-0.5	0	0	0	$\Large\sqrt {\frac 2 {0.5}}$	0
恆星3	1	-1.5	0	0	0	- $\Large\sqrt {\frac 2 {0.5}}$	0

	質量	X	Y	Z	Vx	Vy	Vz
恆星1	1	1	0	0	0	V	0
恆星2	1	$\Large cos(120^\circ)$	$\Large sin(120^\circ)$	0	$\Large cos(210^\circ)V$	$\Large sin(210^\circ)V$	0
恆星3	1	$\Large cos(240^\circ)$	$\Large sin(240^\circ)$	0	$\Large cos(330^\circ)V$	$\Large sin(330^\circ)V$	0

$\Large V=\frac {\sqrt 3} 3$

討論:我們從上面的幾張圖中不難發現,在三體運動中重力場是複雜的,所以彼此間的作用關係不是相互撞上,就是被彈射出去,幾乎無法穩定的存在三恆星系統,所以有則形成巢狀雙星,但是我們也發現巢狀雙星不久就會破壞了平衡,也無法長久的穩定存在。

三、行星系統中的小行星(asteroid)穩定度問題

	質量	X	Y	Z	Vx	Vy	Vz
恆星	1	-0.01	0	0	0	0.5	0
行星	0.01	1	0	0	0	-0.5	0
小行星	0.0001	$\Large cos(\theta)$	$\Large sin(\theta)$	0	$\Large cos(\theta+90^\circ)$	$\Large sin(\theta+90^\circ)$	0

在這個計算中,我們把取為,計算時間到
- 小行星穩定存在的位置
- Result of asteroid at $60^\circ$ and - $60^\circ$

。

Result of asteroid at $90^\circ$

Result of asteroid at $180^\circ$

Result of asteroid at $240^\circ$

四、拉格朗日點的探討

拉格朗日點:在兩個天體的軌道面上，所特有的一些穩定點。例如，與太陽和木星的連心線夾60度的地方，各有一個拉格朗日點，如果有小行星在這兩個拉格朗日點上，它會在此點附近振盪，但不會離開這些點，而特洛伊小行星 (Trojan asteroids) 就是位在這兩個區域。事實上,拉格朗日點不只這兩個,由於其他的拉格朗日點不很穩定，位在其他拉格朗日點上的小天體，稍受擾動就會離開它位置。所以我們只考慮在60度和-60度的這兩點。

	質量	X	Y	Z	Vx	Vy	Vz
太陽	1.99E30	0	0	0	0	0	0
木星	1.90e27	X	0	0	0	V	0
小行星	$\Large\longrightarrow 0$	$\Large cos(\theta)X$	$\Large sin(\theta)X$	0	$\Large cos(\theta+90^\circ)V$	$\Large sin(\theta+90^\circ)V$	0

$\Large X=5.2AU$ , $\Large V=\frac {2 \pi R}T=\frac {2\pi \time5.2AU} {3.74\time10^8}$

在這個計算中,我們把 $\Delta$ t取為86400秒,計算時間到500年

以下的四張圖是用複數坐標系的坐標轉換觀念來畫的,圖中是以太陽為轉動坐標系的原點,太陽-木星連線為X軸,這樣可以方便於觀察它們之間的相互運動關係。
拉格朗日點
Result of"拉格朗日點"

從下面的三張圖我們可以發現小行星會在轉動坐標系中,相對於木星做來回的運動。
Result of asteroid at $90^\circ$

Result of asteroid at $180^\circ$

Result of asteroid at $240^\circ$

四-(1)、行星質量對拉格朗日點穩定性的影響

固定變因:小行星在拉格朗日點上

	質量	X	Y	Z	Vx	Vy	Vz
太陽	1.99E30	0	0	0	0	0	0
木星	m	X	0	0	0	V	0
小行星	0.01m	$\Large cos({60^\circ})X$	$\Large sin({60^\circ})X$	0	$\Large cos({150^\circ})V$	$\Large sin({150^\circ})V$	0

$\Large X=5.2AU$ , $\Large V = \frac {2 \pi R}T=\frac {2\pi \time5.2AU} {3.74\time10^8}$

以下的八張圖仍然是以複數坐標系的坐標轉換觀念來畫的,是以太陽為轉動坐標系的原點,太陽-木星連線為X軸。

木星質量 $\Large m=\frac 1 {1000}$ 太陽質量

木星質量 $\Large m=\frac 1 {100}$ 太陽質量

木星質量 $\Large m=\frac 2 {100}$ 太陽質量

木星質量 $\Large m=\frac 3 {100}$ 太陽質量

木星質量 $\Large m=\frac 4 {100}$ 太陽質量

在下面這一張圖中,我們也可以看到所謂的重力彈弓,在圖中可以發現當小行星經過木星時,就被加速而甩了出去。

木星質量 $\Large m=\frac 5 {100}$ 太陽質量
由上面的例子可知當行星質量越大時,對於恆星系統的影響越大,相對的系統越不穩定,如果今天木星的質量太大的話,它對地球的擾動就越明顯,這會使得我們處在一個不穩定的重力場中,對我們來講當然不是一件好事。今天的木星質量大約是太陽的千分之一,所以還不致於會使地球和小行星一樣發生嚴重的擾動現象,所以我們至今仍能很穩定的存在於太陽系中。

五、太陽系小行星長期軌道演化:

在做太陽系小行星長期軌道演化的結果時,遇到了一個難題就是,在計算過後常會有資料行數到了gnuplot讀取(1500萬行)後會出現當機或者相類似情形,而變化常常就在於後方,所以現今的努力方向是以資料的分段繪出或者將大批資料精簡以求完整的呈現資料的原貌。

捌、結論:

1.雙星系統中行星的運行軌跡:在雙星系統中,行星的初始條件若沒有經過特殊設計的話,通常無法穩定存在,尤其在雙星附近時更是如此,所以如果有行星能穩定存在於雙星系統的話,有則形成內行星或外行星。

2.三恆星系統的運動:在我們的模擬中,三恆星似乎不可能存在,因為三體問體的重力場複雜度很高,任一顆恆星彼此都受到其它兩顆恆星的影響,但是我們又能發現三恆星初期可能形成巢狀雙星,但後來還是無法穩定的運行,所以三恆系統星是一個不穩定的系統。

3.行星系統中的小行星穩定度問題:在行星系統中,小行星可以穩定存在的位置我們不難發現是在和恆星-行星夾六十度的地方,其餘的地方小行星的軌跡變化明顯的比在拉格朗日點的小行星來得大。

4.拉格朗日點的探討:從上一個問題中我們已經知道拉格朗日點是在 $60^o$ 和- $60^o$ 的地方,我們從上面的圖可以發現,其實小行星會在它們的重力平衡點做運動。那如果把行星的質量變大的話,我們可以發現行星對小行星的擾動就越明顯,最後還是難逃被甩出去的命運。

玖、未來展望

在今天,科學家預測太陽大約還有五十億年的壽命,也就是說,我們大概能跟太陽一起共存亡五十億年之久,但是,又有誰能保證說那天不會有一顆像六千五百萬年前一樣使地球發生大滅絕的巨大隕石或其它的天體撞上地球,而使人類也發生大滅絕,假使我們能預測到幾千萬甚至於幾億年後會有那顆不穩定的天體撞上地球,而能事先去摧毀它或使它改向,那麼將能使地球免於再一次的大滅絕。那如果我們能在程式中深入探討小行星接近地球時的引力問題,那麼對於這一類問題會有很大的幫助。

拾、參考文獻
- 國立科學教育館競賽活動
- 恆星系統與三體運動穩定性的模擬探討, 高雄中學:陳威尹&謝俊駒, 第四十屆高中地科科展第一名
來源: http://www.ntsec.gov.tw/activity/race-1/40/hight/16.pdf
來源: http://www.khjh.kh.edu.tw/science40/高中/高中地科1/高中地科1.htm
陳威尹的科展作品內文(網頁型式)
- 逆向操作, 清水高中:黃姝琳&王衍翔, 第四十二屆高中地科科展
來源: http://www.ntsec.gov.tw/activity/race-1/42/pdf/e/5/040509.pdf
- 小行星對於地球原始海水的貢獻, 游大立, 碩士論文
- MimeTeX-Examples ..
- More examples of mimeTeX

彰中地科館

地球科學課程

網路資源

大家一起研究

大立老師的收藏

PmWiki

pmwiki.org

高中研發人才培育學生研究計畫

利用網路上免費軟體進行天體運行的計算

作者: 黃茂誠, 許鈞筌

指導老師: 游大立

1.雙星系統中行星的運行軌跡:在雙星系統中,行星的初始條件若沒有經過特殊設計的話,通常無法穩定存在,尤其在雙星附近時更是如此,所以如果有行星能穩定存在於雙星系統的話,有則形成內行星或外行星。

3.行星系統中的小行星穩定度問題:在行星系統中,小行星可以穩定存在的位置我們不難發現是在和恆星-行星夾六十度的地方,其餘的地方小行星的軌跡變化明顯的比在拉格朗日點的小行星來得大。