基于Matlab常微分方程數(shù)值解的分析與比較

馬澤濤

[摘 要]本文借助Matlab常微分方程求解工具箱，從時間與精度兩個方面對剛性和非剛性方程的數(shù)值求解進行分析與比較，進而對常微分方程的求解給出一般的建議。

[關(guān)鍵詞]常微分方程；時間；精度； 剛性方程；非剛性方程

[中圖分類號] G642 [文獻標識碼] A [文章編號] 2095-3437（2017）12-0050-03

常微分方程現(xiàn)代數(shù)學(xué)的一個重要組成部分，是研究自然科學(xué)、社會科學(xué)中事物運動變化規(guī)律最基本的數(shù)學(xué)理論和方法，是各種應(yīng)用型學(xué)科和數(shù)學(xué)理論研究都不可缺少的工具。但絕大多數(shù)常微分方程都不能求得解析解。因此，要分析與比較不同數(shù)值解法[1][2][3][4]，面對用戶不同層次的要求，找到最優(yōu)解法是急需解決的問題。根據(jù)常微分方程理論，一階常微分方程初值問題

上述的方法在求解剛性方程和非剛性方程時，有各自的優(yōu)點，又都存在著一定的局限性。例如：一般情況下，若采用函數(shù)ode45求解剛性方程時，結(jié)果會很糟糕甚至是死機。即使求得了相應(yīng)的數(shù)值解，所得的結(jié)果也會和精確解相差較大。因此，在將常微分方程分為剛性方程和非剛性方程后，針對它們各自的特殊性，從求得數(shù)值解所用的時間和精度這兩個方面，來對比總結(jié)出較好的方法。

一、基于時間的比較

上述ode系列函數(shù)所基于的數(shù)學(xué)思想方法不同，主要是步長取法上的差異，而步長的選取對數(shù)值計算有很大的影響。步長選取較大雖然會減少方程求解在計算機中的運行時間，但保證不了求解精度；步長選取過小，雖然保證了求解的精度，卻會引起計算機中的多次舍入誤差的積累及運行時間的增加。因此，產(chǎn)生了變步長的算法，它能很好地平衡這兩方面。然而，有的用戶僅需要了解方程解的大概分布情況，對精度的要求不太高，因此針對時間方面，我們做了如下的比較。

對于剛性方程（2）來說，首先，基于ode15s，ode23s，ode23t和ode23tb這四個函數(shù)都只能求得低精度的數(shù)值解情況下，又從上表的各函數(shù)運行所需時間可以看出，函數(shù)ode15s和ode23t所需時間明顯少于ode23tb和ode23s，并且函數(shù)ode23tb又明顯優(yōu)于ode23s。

其次，在這4個函數(shù)都是隱格式情況下，函數(shù)ode15s采用的是多步長法，而ode23s，ode23t和ode23tb采用的都是單步長法。從上表可以看出，多步長法在運行時間方面具有較大的優(yōu)勢。在單步法中，函數(shù)ode23t的運行時間又少于其他函數(shù)。因此，在求解這類剛性方程且求解精度相差不大的四種函數(shù)中，函數(shù)ode15s是最佳的求解辦法，ode23t次之。

而對于非剛性方程（3），從上表中可以看到，函數(shù)ode113運行時間最少，ode23次之，ode45最多。

這是因為，這三個函數(shù)中，函數(shù)ode45與ode23采用的都是變步長的單步法，而ode113采用的是多步法。變步長法即在解變化大的地方采用小步長，解變化小的地方采用大步長。但許多理論和研究都表明，在解變化小的地方仍然不能采用大步長，否則會導(dǎo)致誤差的急劇增加。于是，為了保證函數(shù)ode45和ode23的求解精度，并沒有在解變化小的地方采用大步長，從而導(dǎo)致了其運算次數(shù)的增加。相應(yīng)的，和多步法的函數(shù)ode113相比，運算時間上也不占優(yōu)勢。因此，求解非剛性方程時，函數(shù)ode113運行時間最少，是最佳的求解辦法，ode23次之，ode45最慢。

二、基于精度的比較

其次，單純地要求運行時間少并不利于更好的研究解的性質(zhì)。在精度方面，我們可以發(fā)現(xiàn)，多步法提出了一種新的逼近，但是許多必需的計算將被求解過程中已經(jīng)算得的值的插值所代替。且采用多步長算法的函數(shù)的求解精度往往比采用單步長算法的函數(shù)的求解精度要高。于是，在不考慮時間因素的情況下，作了如下比較。

對于剛性方程（4）來說，首先，可以知道函數(shù)ode23t和ode23tb在求解該方程時失效。因此，在求解剛性方程時，在四個求解剛性方程的函數(shù)中，函數(shù)ode23t和ode23tb不是首要選擇。

再討論函數(shù)ode15s和ode23s的精度。函數(shù)ode15s采用的是多步法，而ode23s采用的是單步法。多步法利用前面若干個節(jié)點得到下一點的近似值，相比只利用了前一個節(jié)點得到下一點的單步法得到的數(shù)據(jù)更具說服力與科學(xué)性。再比較二者的二范數(shù)，由于0.2058<0.3885，所以函數(shù)ode15s是求解這個微分方程的最佳方法，ode23s次之。

而對于非剛性方程（5），首先可以發(fā)現(xiàn)這三個函數(shù)都能求得相應(yīng)的二范數(shù)，故函數(shù)ode45，ode23，ode113都可用來求解此方程。其次，根據(jù)二范數(shù)的比較可知，函數(shù)ode45所求精度最高，ode113其次，ode23最差。

由于這三個函數(shù)都是顯格式，我們進行其他算法上的綜合比較。函數(shù)ode45和ode23都是單步法，其中函數(shù)ode45采用的是Runge-Kutta法的四、五階算法，函數(shù)ode23是Runge-Kutta法的二、三階算法，函數(shù)ode45階數(shù)更高，則所求精度也越高；而函數(shù)ode113是多步法，故函數(shù)ode113優(yōu)于ode23。然而，函數(shù)ode113是變精度變階次的算法，在這里它所求階次比函數(shù)ode45要低，故此例中函數(shù)ode45更優(yōu)。這從另外一個方面也可以說明，可能會有函數(shù)ode113所求精度比ode45高的情況出現(xiàn)。

三、基于實際應(yīng)用的分析

應(yīng)用1：導(dǎo)彈追蹤問題

從上表的運行時間上可以看出，函數(shù)ode45，ode23與ode113所需時間和我們之前所得結(jié)論相吻合。然而，從二范數(shù)的比較中可以發(fā)現(xiàn)函數(shù)ode23所求精度高于函數(shù)ode113，與之前所得結(jié)果有一點出入，這是由于函數(shù)ode113是變精度變階次的算法，在求解該方程時，比Runge-Kutta法的二、三階算法所求精度更低。因此，在面對非剛性方程時，在精度方面要對函數(shù)ode113格外注意，在條件允許的情況下，可以嘗試其他函數(shù)取更優(yōu)法。endprint

應(yīng)用2：化學(xué)反應(yīng)問題

從上表的所需時間看出，這四個函數(shù)均與方程（2）所得結(jié)論相吻合。然而，在精度方面由于方程（4）我們所選方程的特殊性，函數(shù)ode23t和ode23tb不能用于求解該方程，而在此化學(xué)反應(yīng)中，這四個函數(shù)都可用于求解此方程，并且函數(shù)ode15s所求精度最高，其次分別是函數(shù)ode23t，ode23s及ode23tb。究其原因，也是由于函數(shù)ode15s采用的是多步長算法，而函數(shù)ode23t，ode23s及ode23tb采用的是單步長算法。所以，在選用函數(shù)ode23t和ode23tb時，也應(yīng)該格外注意其使用其局限性，實際上這就需要我們關(guān)注方程的求解區(qū)間。

四、結(jié)語

通過以上幾個例子及應(yīng)用，我們可以發(fā)現(xiàn)：在已知方程剛性與否的情況下，在運用Matlab中ode系列函數(shù)求常微分方程數(shù)值解時，對于大多數(shù)剛性方程，若只要求運行時間少但對精度要求不高，多步法與單步法相比占有很大的優(yōu)勢，因此函數(shù)ode15s是首選，ode23t次之；而在對精度有一定要求但不考慮時間的情況下，多步法同樣優(yōu)于單步法，從而函數(shù)ode15s是首選，ode23s次之。對于一般的非剛性方程，若只考慮時間，采用變精度變階次的Adams？鄄Bashforth？鄄Moulton校正法的函數(shù)ode113是首選；若只考慮精度，采用自動變步長的R？鄄K方法的函數(shù)ode45是首選，ode113次之。

因此，綜合考慮求解精度和運行時間這兩方面，對于非剛性方程來說，函數(shù)ode113在某些方程中都是最佳選擇，其次考慮嘗試用函數(shù)ode45；而對于剛性方程來說，函數(shù)ode15s是最佳選擇。當然，考慮到還有眾多常微分方程在未能準確判斷其剛性與否的情況下，可以優(yōu)先嘗試函數(shù)ode45，而ode45失效時再嘗試ode15s。

[ 參 考 文 獻 ]

[1] 李慶揚，王能超，易大義.數(shù)值分析（第五版）[M].北京：清華大學(xué)出版社，2008.

[2] 魏明強.一階常微分方程數(shù)值解中四種算法的實例比較[J].中國傳媒大學(xué)報（自然科學(xué)版），2016（2）：41-44.

[3] 孫美玲.常微分方程數(shù)值解法的Matlab計算與可視比較[J].高教學(xué)刊，2016（19）：60-61.

[4] 胡慶婉.常微分方程初值問題的數(shù)值求解及MATLAB實現(xiàn)[J].SCIENCE & TECHNOLOGY INFORMATION，2012（7）：34-35.

[5] 王高雄，等.常微分方程[M].北京：高等教育出版社，2006（7）：120-121.

[6] Richard L.Burden，J.Douglas Faires.Numberical Analysis[J].Higher Education Press Thomson Learning，Inc，2001（3）：348-353.

[7] 宋葉志，等.MATLAB數(shù)值分析與應(yīng)用[M].北京：機械工業(yè)出版社，2014（2）：384-393.

[8] 張曉.Matlab微分方程高效解法：譜方法原理與實現(xiàn)[M].北京：機械工業(yè)出版社，2015（9）：9-12.

[9] 袁東錦.數(shù)值分析=Numberical Analysis[M].南京：東南大學(xué)出版社，2005（8）：158-166.

[10] 李慶揚，等.數(shù)值計算原理[M].北京：清華大學(xué)出版社，2000（9）：372-376.

[責(zé)任編輯：林志恒]endprint