馬澤濤

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

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

[中圖分類號(hào)] G642 [文獻(xiàn)標(biāo)識(shí)碼] A [文章編號(hào)] 2095-3437（2017）12-0050-03

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

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

一、基于時(shí)間的比較

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

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

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

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

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

二、基于精度的比較

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

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

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

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

由于這三個(gè)函數(shù)都是顯格式，我們進(jìn)行其他算法上的綜合比較。函數(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)。這從另外一個(gè)方面也可以說(shuō)明，可能會(huì)有函數(shù)ode113所求精度比ode45高的情況出現(xiàn)。

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

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

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

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

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

四、結(jié)語(yǔ)

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

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

[ 參 考 文 獻(xiàn) ]

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

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

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

[4] 胡慶婉.常微分方程初值問(wèn)題的數(shù)值求解及MATLAB實(shí)現(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].北京：機(jī)械工業(yè)出版社，2014（2）：384-393.

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

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

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

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