用Python做Lorenz attrator的動畫模擬

建議大家去我的blog看這篇文章，效果要更好：http://zhanglei.eu/

最近迷上了Lorenz attractor，於是在閑暇時間，用Python做了Lorenz attractor的動畫模擬，結果還是不錯的。

什麼是Lorenz attractor

簡單來說，Lorenz attractor 是一個混沌的動力系統，系統對初始值非常敏感，即使初始值非常小的不同，隨著時間演化，到最後能造成極大的差別。

這個動力系統首先由法國天氣學家Lorenz發現。他在研究天氣動力系統時，提出了一個非常簡單的模型：

這組方程描述了大氣里溫度梯度，各處壓強，流速等信息，但對真實情況做了太多的假設和簡化，在物理上並沒有參考意義，是一個「玩具模型」。但在數學上，卻帶了了非常有意思的發現。

Lorenz在一次計算時，在兩組初始值相差非常小的情況下，他驚奇的發現隨著時間演化，兩組結果竟然產生了翻天覆地的變化，甚至截然不同。由相同方程支配的系統竟然在初始值差別非常小的情況下產生巨大的不同，這是反直覺的，正如在他隨後的論文中詩意地寫道：

於是Lorenz attractor又被稱為蝴蝶效應。

這裡需要提一下，大家在電影里看到的蝴蝶效應，是藝術家對這一現象的詩意表達，並不是真實發生的現象。

Python求解

首先，我們將求解以下方程：

我們取p=10, r=30, b=3, 兩組初始值分別是（X,Y,Z）=（0，1，0）和（0，1，0.0000001），差別非常小！！！

Python里，我們利用scipy.integrate模塊下的odeint函數進行求解，具體的使用方法可以參看scipy的的documentation。

代碼如下：

from scipy.integrate import odeint import numpy as npdef lorenz(w,t,p,r,b): #定義求解的系統，返回函數組 x,y,z=w return np.array([p*(y-x),x*(r-z),x*y-b*z])t =np.arange(0,50,0.01)#定義求解區間，求解50s，時間間隔0.01strack1=odeint(lorenz, (0.0, 1.00, 0.000000001),t,args=(10.0,30.0,3.0))#得到結果track2=odeint(lorenz, (0.0, 1.00, 0.0),t,args=(10.0,30.0,3.0))

簡單的幾行代碼，就給我們帶來了求解結果，現在需要把結果可視化一下，我們首先利用matplotlib畫出兩組結果的運動軌跡：

代碼如下:

from mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as pltfig=plt.figure(1)ax=Axes3D(fig)line1=ax.plot(track1[:,0],track1[:,1],track1[:,2],"b", markersize=8)line2=ax.plot(track2[:,0],track2[:,1],track2[:,2],"r", markersize=8)plt.show()

我們可以看到，紅色軌跡和藍色軌跡出現了分離，而且軌跡的形狀像一隻蝴蝶，這也是Lorenz要拿蝴蝶舉例子的原因。

軌跡繞著兩個「眼」運動，而且不能脫離兩個「眼「，保持在其附近。這兩個」眼「好像把軌跡吸引過去一樣，這就是為什麼叫Lorenz attractor(吸引子)的原因。

我們來看看X坐標隨著時間的演化：

我們可以看出，兩組結果X坐標在23s左右出現分離，按照相反的路線開始演化，而且變化非常劇烈。還可以看出，運動軌跡其實有一定周期規律，但不是我們常見的周期函數，這也暗示了混沌雖然變化劇烈，難以預測，但其從整體上來說，又有一定的規律性。

Python動畫模擬

可視化的結果只給了我們整體的軌跡曲線，並不能反應每時每刻的運動狀態。於是我們緊接著用matplotlib.animation包進行3D的動畫模擬，主要刻畫出每一時刻的運動狀態。

代碼也很簡單，主要用了matplotlib.animation里的FuncAnimation函數，代碼如下：

from mpl_toolkits.mplot3d import Axes3Dimport matplotlib.pyplot as pltfrom matplotlib import animation #這個包用於動畫製作ax=Axes3D(fig)#生成3D的軸line1=ax.plot([],[],"b:")#初始化數據，line是估計，point是軌跡的頭部point1=ax.plot([],[],"bo",markersize=10)line2=ax.plot([],[],"r:")point2=ax.plot([],[],"ro",markersize=10)images=[]def init():#該函數用於初始化動畫 line1=ax.plot([],[],"b:",markersize=8) point1=ax.plot([],[],"bo",markersize=10) line2=ax.plot([],[],"r:",markersize=8) point2 = ax.plot([], [], "ro", markersize=10） return line1,point1,line2,point2 def anmi(i):#anmi函數用於每一幀的數據更新，i是幀數。 ax.clear() line1=ax.plot(track1[0:10*i,0],track1[0:10*i,1],track1[0:10*i,2],"b:", markersize=8) point1 = ax.plot(track1[10*i-1:10*i,0],track1[10*i-1:10*i,1],track1[10*i-1:10*i,2],"bo", markersize=10) line2 = ax.plot(track2[0:10*i, 0], track2[0:10*i, 1], track2[0:10*i, 2], "r:",markersize=8) point2 = ax.plot(track2[10*i-1:10 * i, 0], track2[10*i-1:10 * i, 1], track2[10*i-1:10 * i, 2], "ro", markersize=10) return line1,point1,line2,point2#animation.FuncAnimation用於動畫製作，輸入init函數和anmi函數，frams是幀數，一共畫500禎，interval是相鄰禎的間隔，單位是毫秒，此處為2毫秒anim = animation.FuncAnimation(fig, anmi, init_func=init, frames=500, interval=2, blit=False,repeat=False)

效果如下：

在這個動畫里，我們清楚的看到混沌現象，兩個小球初期值僅僅相差0.000001，結果在29s左右的時候竟然按照相反的軌跡演繹，之後軌跡也不會重合，非常有意思！

用VPython進行3D實物動畫模擬模擬

到目前為止，雖然動畫的結果基本滿意，也準確反應了Lorenz Attractor的特性。但對於追求精益求精的工程師，目前的3D動畫有點太過簡陋了，要是有更美、更漂亮的3D實物動畫，那就更好了。

於是上網查，果然發現了一個神器——VPython， 它是python環境下比較簡單的3D圖形庫，功能也很強大，對3D動畫的支持也非常好，提供的操作非常多，還能動態的畫出科學曲線，很適合科學工作者去進行數據3D可視化。

這個包的具體使用方法可以參考網站：http://vpython.org/

那麼說干就干！！

我們建立4個窗口，從左到右依次為：X坐標隨時間的變化曲線，3D運動軌跡，X坐標變化的3D呈現，3D曲線在Y-Z平面內的投影，運行效果如下：

Lorenz attractor3D動畫模擬—在線播放—優酷網，視頻高清在線觀看 http://v.youku.com/v_show/id_XMTY2MjM2NDAyNA==.html

效果還是很好看噠！！哈哈哈哈~

代碼可以參看我的blog: http://zhanglei.eu/2016/07/28/%E7%94%A8Pytho%E5%81%9ALorenz-attrator%E7%9A%84%E5%8A%A8%E7%94%BB%E4%BB%BF%E7%9C%9F/

大家有問題可以私信問我哈~

還請大家積極關注我個人創建的微信公眾平台，紅酒學一學，主要講一講我對葡萄酒文化的學習體會，還會總結我的學習筆記，寫幾篇適合新手的入門教程，也會偶爾寫寫其他的酒。請關注這個平台，讓我和大家一同進步，掃描下面二維碼關注。

或者微信搜索：紅酒學一學

http://weixin.qq.com/r/XEyHn4TEqdnWrcD89xlO (二維碼自動識別)