//最近試著自學(xué)了python對初值問題的數(shù)值模擬。

//原理也非常簡單，把微分方程近似為差分方程，然后用循環(huán)語句計算質(zhì)點位置的變化。

//本程序基于python的xlrd，numpy，matplotlib，mpl_toolkits.mplot3D庫實現(xiàn)。

//xlrd可以從excel導(dǎo)入數(shù)據(jù)，numpy處理向量運算，matplotlib和mpl_toolkits繪制圖像。

0. 理論基礎(chǔ)

三體系統(tǒng)是著名的混沌系統(tǒng)，三個質(zhì)點僅受萬有引力作用而產(chǎn)生復(fù)雜的運動，且決定系統(tǒng)運動的微分方程無法求得解析解，只能數(shù)值求解。同時，系統(tǒng)的后續(xù)運動對初值高度敏感，隨著時間推移數(shù)值近似的解誤差也越來越大。雖然對混沌系統(tǒng)的詳細研究相當(dāng)復(fù)雜，有包括龐加萊截面法、功率譜分析法、分形維數(shù)、李雅普諾夫指數(shù)等刻畫混沌現(xiàn)象的方法，但不影響我們基于基本的微分方程對該系統(tǒng)進行簡單的數(shù)值模擬。

設(shè)三個星體質(zhì)量分別為，位置用矢量表示。三體系統(tǒng)的微分方程如下：

在現(xiàn)實中萬有引力屬于相對微弱的作用力(相對其他基本力)，因此引力的作用也通常是質(zhì)量極大的天體在相當(dāng)長的時間里相互作用得以體現(xiàn)(相對地球上的日常而言)。所以這個模型中我們簡單起見直接令，所有量綱相關(guān)的討論都不考慮。

作為數(shù)值模擬的需要，我們將微分方程化為差分方程，以1質(zhì)點為例：

受力

變速

位移

其中是程序中設(shè)定的一個很小的時間間隔。越小則誤差越小，但計算量越大。

I. 準備工作

開始之前你需要確認你的python裝有前面提到的xlrd，numpy，matplotlib，mpl_toolkits幾個庫。如果你使用的IDE是Anaconda，這些包都是自帶的，不需要額外安裝。關(guān)于其他如Pycharm, Visual Studio等環(huán)境我也不是很了解，但總之如果你缺少某個包(例如numpy)，只要在命令行中運行以下命令：

系統(tǒng)就會自動從網(wǎng)絡(luò)上搜索相關(guān)庫并安裝。確認所有需要的庫都安裝后就可以開始了！

II. 開始編程

①包與數(shù)據(jù)導(dǎo)入

我們首先在程序中導(dǎo)入所有需要的庫：

由于本系統(tǒng)需要的初始條件高達21個(每個質(zhì)點的質(zhì)量、初始位置三個分量、初速度三個分量)，我們?nèi)绻胕nput()命令在程序運行時逐個輸入麻煩且不直觀，我們選擇把初值寫在excel文件中并利用xlrd包導(dǎo)入。首先找一個地方創(chuàng)建一個Excel文件DataInput.xlsx，例如在D盤的文檔中。Python文件中利用以下命令：

就可以將數(shù)據(jù)來源指定為DataInput.xlsx的sheet1。在excel中，按如圖方式將三個質(zhì)點的質(zhì)量、初位置、初速度寫入，保存。(你可以根據(jù)自己的喜好指定初始值，圖中的數(shù)據(jù)僅供參考)

接下來我們將這些數(shù)據(jù)存儲在變量m, 數(shù)組r,v中，作為后續(xù)計算的準備工作：

注意雖然我們的數(shù)據(jù)是從第2行第2列開始的，但代碼中是從(1,1)開始導(dǎo)入數(shù)據(jù)的，這是因為大部分編程語言中數(shù)組或類似結(jié)構(gòu)下標都默認從0開始，下標要注意減1.

② 計算準備

我們需要準備九個空列表存放各個質(zhì)點坐標隨時間的變化：

設(shè)置和總步數(shù)：

創(chuàng)建一個三維圖，用于可視化結(jié)果：

我們定義幾個函數(shù)，求點間的施力：

四個函數(shù)的功能分別是：輸入兩個位置返回距離；輸入距離和質(zhì)量返回引力大小；輸入兩個位置返回連接它們的單位矢量(以數(shù)組形式)；調(diào)用引力大小和單位矢量函數(shù)返回引力矢量(數(shù)組形式)。其中，調(diào)用了numpy中的數(shù)學(xué)、數(shù)組運算函數(shù)。

接下來，定義兩個void函數(shù)，用于進行質(zhì)點的移動和受力改變速度：

③ 計算

準備工作都做好了，現(xiàn)在只需要一次循環(huán)語句完成所有計算即可。

④ 繪圖

所有數(shù)據(jù)的計算都完成了，接下來就簡單了，以下命令將幾個質(zhì)點的軌跡繪制在前面定義的三維圖中：

最后輸入預(yù)覽命令

我們就會成功看到前面的封面圖片了！

可以看到，確實很混沌。

如果你希望保存這個圖片，即可以直接在預(yù)覽窗口上點擊保存，也可以用命令

來保存這個結(jié)果。

那么以上就是這樣一個新手向的Python數(shù)值計算筆記了。當(dāng)然本人也只是Python初學(xué)者，偶爾用它來數(shù)據(jù)分析，這個程序可能沒有寫出最簡潔高效的代碼，如果你有更好的方法也歡迎指出。