XPPAUT在Lorenz方程中的應用

張 鑫,陳欣怡,2,魏章志

1.南京財經大學應用數學學院,江蘇南京,210023;2.蘇州美萊美容醫院,江蘇蘇州,215000;3.宿州學院動力系統與控制研究中心,安徽宿州,234000

隨著科學技術的飛速發展,在研究數學問題時,計算機的輔助功能變得越來越重要。目前可以使用的軟件數不勝數,比如:Matlab、Maple和Mathematica。XPPAUT軟件是由匹茲堡大學的Bard Ermentrout[1]開發的,并將該軟件的使用指南記錄在其著作中,而孝鵬程等[2]將該著作翻譯成了中文。XPPAUT是一款免費的動力系統分析工具,主要功能分為兩個部分:方程數值求解部分(XPP)和分岔分析部分(AUTO)。 雖然XPPAUT使用方便且功能強大,可以對常微分方程、時滯方程和隨機方程進行數值模擬、分岔分析,自帶的動畫功能可以更加生動形象地體現出問題的特點,但是其在國內普及不是很廣泛,國內介紹相關使用方法的教材也很少。 下面主要研究XPPAUT在常微分方程中的應用。

1 XPPAUT簡介

1.1 XPPAUT界面介紹

ODE文件是XPPAUT所使用的最基礎的文件,這些文件的拓展名都是.ode,這里通稱ODE文件。 只要在啟動Xming的前提下,將ODE文件拖到XPPAUT上方就可以打開該文件,這是使用XPPAUT最簡單的方法。 事實上,文檔也可以代替ODE文件。只要創建一個文本文檔,在其中輸入符合XPPAUT運行規則的代碼,就可以被XPPAUT正常讀取和運行。

在ODE文件中,通常使用par定義參數,用init定義初始值,且參數和初始值的默認值均為0。 代碼不用區分大小寫字母的,且必須以done作為結尾。 此外,ODE文件還有一些常用的符號,以其中兩個為例,代碼每一行至多有1 024個字符,用符號“”可以換行繼續編輯;“#”可以加入一些注釋,但不是文件必須有的部分。

XPPAUT界面有兩部分,第一部分稱為XPP界面。將ODE文件拖到XPPAUT上方后,打開的界面為XPP界面,主要完成數值求解任務,如圖1所示。

圖1 XPPAUT的XPP界面

(1)在界面最上方是快捷菜單欄,包含:ICs、BCs、Delay、Param、Eqns和Data六個部分。 ICs可以設置變量初值,還可以改變繪圖中橫縱坐標代表的變量;BCs用于邊界值問題;Delay可以將時滯方程初始化;Param用于改變參數;Eqns可以查看微分方程;Data可以查看數據。

(2)在其下方是輸入窗口,必要時可輸入內容。

(3)左側一列是主菜單欄,是XPPAUT最重要的部分。幾個常用的菜單如下: Initialconds用來設置初始條件及求解方程;Dir.Field/flow可繪制方向場;Window/Zoom可調整繪圖顯示區域;Graphic stuff用于圖像設置;File可以操作文件;Erase可以清除圖像;Xivst用來繪制指定量的時間序列圖;Restore可以重新繪制圖像。

(4)右側最大空白的地方是繪制窗口,可以顯示所繪制的圖像。

(5)窗口底部為信息欄,會有一些文字提示。

第二部分是AUTO界面。點擊File-Auto進入AUTO界面,如圖2所示。 AUTO可用于分岔和延續,具有非常強的技巧性,運行過程中可能會遇到各種各樣的問題。

圖2 XPPAUT的AUTO界面

(1) 左側一列為菜單欄,其中Axes可以對繪圖進行設置,Numerics用來設置數值參數以及分岔方向,Run是畫出分岔圖的按鈕,reDraw用來重新繪制圖形,ABORT可以中斷計算。

(2)左側圓形為狀態窗口,可以看出當前點的穩定性。當加號“+”在圓外時,表示點是不穩定的,當“+”在圓內時,表示點是穩定的。

(3)右側最大的空白區域為繪圖窗口,主要是顯示分岔圖。

(4)最下方為提示信息窗口,有簡單的提示操作文字。

(5)在提示信息窗口上方為信息顯示窗口,可以顯示當前點的信息。

第三部分是IPAD上的XPPAUT界面。IPAD也有XPPAUT,界面比電腦端簡單,打開XPP軟件,如圖3所示。 同樣有快捷菜單欄、主菜單欄、繪制窗口和信息欄。 通過點擊右上角“+”建立ODE文件,如圖4所示。

圖3 IPAD的XPPAUT界面

圖4 創建ODE文件界面

1.2 XPPAUT的研究現狀

近30年來,國內外與XPPAUT軟件有關的學術成果不是很多,而且從所占比例上來看以期刊論文為主,如圖5所示。 而研究XPPAUT在常微分方程中的應用的相關內容更是少之又少。

圖5 XPPAUT在各頻道檢索量統計

從研究發展的趨勢上來看,近幾年XPPAUT的研究量較20年前相比有明顯的增加,如圖6和圖7所示。 在2004年以前,并沒有與XPPAUT相關的學位論文,與其相關的期刊論文也是屈指可數。 2004年后的幾年至現在,學位論文和期刊論文都呈現出折線形發展趨勢,整體呈上升狀態。

圖6 與XPPAUT相關的學位論文發展趨勢曲線

圖7 期刊發表XPPAUT相關論文趨勢曲線

在國內,主要應用此軟件到數理學科和化學學科中,其次是生物科學和工業技術。

2006年,李義龍[3]在其博士論文《幾類生物數學模型的定性與分支研究》中用XPPAUT對幾類生物數學模型進行數值模擬; 2013年,劉宣亮[4]發表了《XPPAUT在常微分方程教學中的應用》,主要將XPPAUT中的部分功能用于常微分方程的課堂教學中,例如,“方向場”“垂直等傾斜線”和“軌線”等概念,可以用圖來幫助理解,還可以用XPPAUT的動畫功能演示振動問題。 由此,不僅可以讓學生對所學內容理解更加深刻,還可以調動學生學習的積極性。 2015年,孔維旺等[5]在《XPPAUT在系統生物學中的應用》寫到XPPAUT是系統生物學必備工具之一,有著十分可觀的前景,并用該軟件研究了系統生物學中的四類典型問題,分別是求系統的數值解、相平面分析、分岔分析和模擬網絡模體的信息處理功能。由此可見,XPPAUT可以靈活地運用在這門數學學科中。外文文獻XPPAUT軟件主要應用在理科研究中,其次是計算機科學。

1987年,XPPAUT誕生,創造者將其使用指南記錄在Simulating,AnalyzingandAnimatingDynamicalSystems:AGuidetoXppautforResearchersandStudents中,才有了之后XPPAUT在各種領域中的應用。書中不僅有該軟件詳細的安裝步驟,還有界面介紹、一些功能的使用方法以及在具體學科和問題中的應用。 Frank W?rnle等[6]在2005年利用XPPAUT及延拓技術、分岔理論對鐵磁諧振電路的動力學進行了分析。2015年,Alexandra Teslya等[7]在其博士論文中用XPPAUT對捕食模型進行數值分析。此類相關研究很多,將不一一贅述。

很多軟件可以用來解決與微分方程有關的問題,例如利用Matlab中的MatCont安裝包分析非線性動力學方程組的周期解以及分支圖等、DDE-BIFTOOL安裝包分析時滯微分方程的動力學行為等,還可利用Maple以及Mathematica等軟件, 但是這些軟件在解決數值問題時,運行十分緩慢。 與上述軟件相比,XPPAUT有著非常明顯的優勢。首先,其數值積分的速度較快,語法簡單,還有一些方便地快捷鍵。其次,XPPAUT是一款免費的軟件,不需要每年更新證書,并且XPPAUT還有手機端和平板電腦端的程序,對于移動時代的今天來說是巨大的便利,是一個不可忽略的優勢。 此外,XPPAUT還可以自定義字體、顏色以及背景等,增強了趣味性。 最后,該軟件的動畫功能生動形象,有利于在教學領域的應用,調動學生學習的積極性。

但是,XPPAUT的缺陷在于其圖像后處理功能較弱,無法滿足一些對圖像質量要求較高的情況,此時往往需要借助第三方軟件的幫助。雖然XPPAUT存在著一些缺點,但是在其超高的計算能力面前,這些缺陷是可以被忽視的,無法改變它是一個卓越的軟件這個事實。

2 XPPAUT在常微分方程中的應用

下面以Lorenz方程為例說明XPPAUT在常微分方程中的應用

2.1 Lorenz方程

1963年,美國氣象學家E.N.Lorenz在研究氣象預報模擬氣候變化的時候,第一次發現了如下的方程

(1)

其中,a,b,c為系統參數;x,y,z∈Rn為系統變量[8]。

系統(1)是利用三維截斷這一方式處理無窮維動力系統的 Rayleigh-Benard 熱對流問題而得到的,該方程具有“蝴蝶效應”。 該效應說明,事物發展的結果往往對初始條件極為敏感,初始條件的微小變化,可能會導致結果產生巨大的差異,這是世界上的第一個混沌模型,一般稱為Lorenz方程。 Lorenz方程的提出,奠定了混沌系統發展的基礎,并且推動了混沌學的發展和實際應用。

自此以后,人們開始了對于混沌系統的研究,混沌系統因對初值敏感、復雜度高、內在隨機性和遍歷性等特征,因此被廣泛應用于保密通訊、神經網絡等各個研究領域。 隨后很多具有混沌性質的三維系統被提出,比如Chen系統、 Lu系統、 Chua電路系統以及Rossler系統等[9]。 近年來,四維超混沌系統相繼被發現并得以研究,使得混沌系統的混沌特性得以深入研究,大力促進了混沌系統在相關領域的應用[9-11]。

主要是將XPPAUT應用到一般的Lorenz方程中,對模型進行數值模擬并畫出部分圖像,然后求出平衡點。

2.2 Lorenz方程的性質

性質1方程(1)在變換(x,y,z)→(-x,-y,-z)下仍然保持不變,即方程關于變量z對稱。

性質2方程(1)在c<1時關于原點是全局一致漸進穩定的,同時有

(2)

2.3 數值模擬

考慮模型(1),令a(y-x)=0,-xz+cx-y=0,xy-bz=0, 可得方程平衡點[8]為:

O(0,0,0),

令a=10,b=2.5,c=40,得出模型(1)的平衡點為:

O(0,0,0),

P+(9.874 2,9.874 2,39),

P-(-9.874 2,-9.874 2,39)。

再用XPPAUT對該模型進行數值模擬,畫出部分圖像并用XPPAUT求出其平衡點,步驟如下[12]:

(1)編寫代碼。

先在文件夾中創建一個文本文檔,該文檔是ODE文件,然后在文檔中輸入以下內容并保存:

#Lorenz1

x′=a×(y-x)

y′=-x×z+c×x-y

z′=x×y-b×z

para=10,b=2.5,c=40

initx=-7.5,y=-3.6,z=30

done

以上代碼中,定義參數a=10,b=2.5,c=40,并且定義了初始值x=-7.5,y=-3.6,z=30。

(2)在XPPAUT中運行文本文檔。

按照1.1的步驟運行上述文本文檔。

(3)計算平衡點。

用XPPAUT演示求解平衡點,以P-平衡點為例。 在運行狀態下,點擊Sing pts-(G)o,會彈出窗口,如圖8(a)。點擊YES,會彈出新的窗口,如圖8(b)。點擊YES,得到Equilibria窗口,如圖8(c)。 其中下方的X=-9.874 2,Y=-9.874 2,Z=39表示模型的負平衡點P-(-9.874 2,-9.874 2,39)。

圖8 利用XPPAUT計算平衡點

(4) 在XPPUT中進行如下操作。

啟動XPPAUT后,點擊Initialconds-(G)o,得到X vs T的圖像,點擊Window/ zoom-(F)it,自動調整繪圖窗口的顯示區域,可以顯示完整圖像,如圖9所示。

圖9 X vs T圖像

點擊ICs,選擇X,Z,再點擊xvsy,得到相平面Z vs X,其中A、B為平衡點,如圖10所示。 同理得到相平面Z vs Y,其中C、D為平衡點,如圖11所示。

圖10 相平面Z vs X

圖11 相平面Z vs Y

點擊Dir.field/ flow-(S)caled Dir.Fld,敲擊Enter接受默認值16為網格值,可得方向場,如圖12所示。用IPAD也可以畫出Lorenz方程的相平面,如圖13所示。

圖12 方向場

圖13 IPAD上Lorenz方程的相平面

Lorenz方程是三維的,點擊ICs,選擇X、Y、Z,點擊xvsy,可得三維相面,如圖14。 從三維相圖可以看出該軌跡幾乎要撐滿整個立方體區域,這是不穩定的平衡點導致的,使軌道在一定程度上拉伸、折疊,這是Lorenz方程的動力學特性的奇特之處。

圖14 Lorenz系統的三維相圖

3 總 結

對XPPAUT軟件的基礎介紹,同時進行了應用舉例選取了常微分方程中的經典Lorenz方程,列舉了一些相關的性質和定理,通過理論計算分別求出該模型的平衡點并用XPPAUT對這兩個模型進行數值模擬,畫出相平面、方向場和軌跡圖等。

常微分方程的學習及教學,離不開軟件的幫助,而這個專門用于微分方程的軟件XPPAUT,有一般的數值模擬,可以繪制圖像。 只選取了兩個較為簡單的模型,實際上該軟件還可以用于更加復雜的其他模型中,如分岔分析和動畫功能。對于其應用還可以進一步研究,之后可以將該軟件應用到時滯微分方程以及其他類型的方程中。 希望一些已經解決或者待解決的問題,通過XPPAUT的使用,會變得更加簡單,可以更快地解決。

XPPAUT在國內的普及率不高,遠低于我們所熟知的Matlab等軟件。期待在未來的日子里,越來越多的國內學者來學習和運用該軟件,讓其在生物、物理以及其他領域中發光發熱。