基于MATLAB 的常微分方程的求解及數(shù)值仿真

彭東海

（中山職業(yè)技術(shù)學(xué)院素質(zhì)教育中心，中山528404）

0 引言

常微分方程是描述客觀世界運動變化的有效工具，在自然科學(xué)、工程技術(shù)、社會學(xué)等領(lǐng)域有廣泛的應(yīng)用，是高等數(shù)學(xué)課程的重要內(nèi)容之一。對于一些比較特殊的常微分方程，如線性方程、可分離變量的方程、存在積分因子的方程等，人們可以求出其解析解；對于稍微復(fù)雜一些的方程，可以借助微分方程定性理論分析解的性態(tài)，如漸近性、穩(wěn)定性、周期性等，而這部分內(nèi)容，我國現(xiàn)行高等數(shù)學(xué)教材多不涉及；對于更一般的常微分方程問題的求解常常是困擾人們的一個難題。在工程技術(shù)和社會科學(xué)等應(yīng)用領(lǐng)域，人們往往借助于現(xiàn)代計算機技術(shù)，對常微分方程進行數(shù)值求解或者數(shù)值仿真，從而得到相應(yīng)的微分方程解的基本性態(tài)。MATLAB 是用于高性能數(shù)值計算和可視化的商業(yè)數(shù)學(xué)軟件，其將數(shù)值分析、矩陣計算、信號處理、符號計算和圖形圖像集成到一個易于使用的交互式環(huán)境中。在此環(huán)境中，通過調(diào)用不同模塊或工具箱，無需太多傳統(tǒng)編程即可表達出來，使得其在工程計算、信號處理、金融建模設(shè)計與分析等領(lǐng)域有廣泛的使用，此外，MATLAB 還是美國數(shù)學(xué)建模競賽最常用的參賽使用軟件。

在常微分方程求解方面，針對不同類型的方程與不同問題，MATLAB 提供了有效的求解器，比較典型的如 ODE45、DDE23、BVP4C、BVPINIT 等，分別可以對常見的常微分方程、時滯微分方程、邊值問題、初值問題進行相應(yīng)的求解。結(jié)合我們多年的教學(xué)經(jīng)驗和學(xué)生使用MATLAB 的情況，本文通過實例介紹用MATLAB 相關(guān)求解器求解微分方程及對微分系統(tǒng)作數(shù)值仿真的基本方法及需要注意的問題。

1 求解常微分方程的解析解

MATLAB 有一個豐富的函數(shù)庫可用于求解常微分方程，不僅數(shù)值計算功能強大，由于MathWorks 公司收購了Maple 的內(nèi)核，其符號運算的功能也很強大。

二階微分方程

求解代碼如下：

得到方程（1）的特解：

對于線性系統(tǒng)：

求解代碼如下：

得到系統(tǒng)（2）的通解：

2 微分方程數(shù)值解

MATLAB 有許多用于數(shù)值求解常微分方程的工具。這里主要介紹常用的內(nèi)置函數(shù)ode45 和dde23，它們是基于Dormand-Price 提出的顯式Runge-Kutta（4,5）對，屬于單步解算命令。

調(diào)用ode45 求解器求解給定隱式微分方程：

將（3）改寫成矩陣形式，令：

輸入代碼：

得到方程（3）的解的圖像。

圖1 隱式方程的數(shù)值解曲線

當然，對于一般的微分方程還有其他的求解器，如表1。

表1

時滯微分方程相對于不含時滯的傳統(tǒng)微分方程能夠更為準確地描述客觀事物的變化規(guī)律，數(shù)理生態(tài)學(xué)等研究領(lǐng)域也出現(xiàn)越來越多含有時滯的微分方程模型。但是，時滯微分方程的求解往往更加困難。MATLAB 中提供了 dde23、ddesd 等求解器。

下面介紹應(yīng)用dde23 求解時滯微分方程，對方程：

給定y(0)=1，建立一個 MATLAB 的m 文件：

得到解的圖像及相位圖（圖2）。

圖2

由此可以判斷上述方程（4）的解是振蕩的。

通常，一個微分方程往往具有許多解，而我們感興趣的是滿足某些特定條件的解。對于邊值問題，可以用bvp4c 求解，只須先將微分方程寫成一階微分系統(tǒng)。例如，給定微分方程：

編寫bvpexa.m 文件如下：

即可得問題（5）的數(shù)值解。

圖3

3 常微分方程方向場的繪制

微分方程方向場可以給出微分方程解的幾何解釋，從微分方程本身的特性了解到它的任一解所應(yīng)具有的某些幾何特征，或者近似的描繪出微分方程解的積分曲線，這就是微分方程方向場的應(yīng)用所在。

給定方程：

先編寫一個m 文件：

然后通過下面主程序代碼：

即可得到方向場如下:

圖4

圖中4 所示的圓點是方程（8）的初始條件下的數(shù)值解，而實線所示為方程的解析解。

的解。

Lyapunov 穩(wěn)定性定理若有原點的鄰域U和一個正定（負定）函數(shù)V(x)，使得V˙(x)是半負定（半正定）的，則系統(tǒng)（9）的零解是穩(wěn)定的，且使得V˙(x)是負定（正定）時，系統(tǒng)（8）的零解是漸近穩(wěn)定的。

下面是對微分系統(tǒng)零解的穩(wěn)定性仿真的方法。對給定的系統(tǒng)：

分析構(gòu)造Lyapunov 函數(shù)：

顯然這是正定的，同時V(x) 的全導(dǎo)數(shù)：

根據(jù)上述定理知系統(tǒng)（9）是零解是漸近穩(wěn)定的。

仿真先編寫一個phase.m 程序：

然后通過調(diào)用MATLAB 主程序：

得到其穩(wěn)定性圖像仿真。

圖5

通過對圖5 的觀察，可以很直觀的理解系統(tǒng)（9）的零解是漸近穩(wěn)定的。

借助MATLAB 的仿真繪圖，對考察二維微分系統(tǒng)參數(shù)的分支及極限環(huán)的存在性時，可以作為理論分析的有效補充，也是很有幫助的。

例如考慮下面微分系統(tǒng)（11）在λ取不同值時零解的穩(wěn)定性

其中λ∈R。

通過調(diào)用phase.m 可以繪制團如圖6。

圖6

通過觀察可以發(fā)現(xiàn)λ＞0 時，原點是穩(wěn)定的；λ＜0時，原點是不穩(wěn)定的，但是當|λ|足夠小時，系統(tǒng)存在極限環(huán)。

4 結(jié)語

借助MATLAB 軟件不僅能進行符號計算求解特定類型的常微分方程的解析解，如常系數(shù)線性方程，變量分離型方程，而且能在解析解無法求出時，得到數(shù)值解，如非線性初值問題、邊值問題、時滯方程等。更重要的是，MATLAB 軟件能通過圖形圖像，直觀地展示常微分方程解的性態(tài)，如穩(wěn)定性、漸近性、極限環(huán)的存在性等。這些，我們通過MATLAB 的幾個典型求解器都做了一一展示。此外，針對實際應(yīng)用中可能遇到的誤差問題，可以通過optimset 對求解器中參數(shù)option 進行設(shè)置，能夠?qū)⒄`差控制在我們要求的范圍之內(nèi)。對于剛性的微分方程，也要注意求解器的選擇，否則會造成計算機處理時間過長等問題。另外，在作微分系統(tǒng)定性分析時，數(shù)值仿真雖然可以起到直觀、快速判斷的效果，但是，要得到嚴格、完備的定性分析結(jié)果，我們還是需要回到理論計算上來。

通過MATLAB 的軟件我們可以比較方便地實現(xiàn)微分方程的求解和相關(guān)仿真，為學(xué)生在學(xué)習(xí)掌握常微分方程的解法和理論時，提供必要的幾何直觀結(jié)果，能激發(fā)學(xué)生學(xué)習(xí)興趣，對提高學(xué)生在微分方程的學(xué)習(xí)效果具有明顯的促進作用。