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

彭東海

（中山職業技術學院素質教育中心，中山528404）

0 引言

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

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

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

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

二階微分方程

求解代碼如下：

得到方程（1）的特解：

對于線性系統：

求解代碼如下：

得到系統（2）的通解：

2 微分方程數值解

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

調用ode45 求解器求解給定隱式微分方程：

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

輸入代碼：

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

圖1 隱式方程的數值解曲線

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

表1

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

下面介紹應用dde23 求解時滯微分方程，對方程：

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

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

圖2

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

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

編寫bvpexa.m 文件如下：

即可得問題（5）的數值解。

圖3

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

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

給定方程：

先編寫一個m 文件：

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

即可得到方向場如下:

圖4

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

的解。

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

下面是對微分系統零解的穩定性仿真的方法。對給定的系統：

分析構造Lyapunov 函數：

顯然這是正定的，同時V(x) 的全導數：

根據上述定理知系統（9）是零解是漸近穩定的。

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

然后通過調用MATLAB 主程序：

得到其穩定性圖像仿真。

圖5

通過對圖5 的觀察，可以很直觀的理解系統（9）的零解是漸近穩定的。

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

例如考慮下面微分系統（11）在λ取不同值時零解的穩定性

其中λ∈R。

通過調用phase.m 可以繪制團如圖6。

圖6

通過觀察可以發現λ＞0 時，原點是穩定的；λ＜0時，原點是不穩定的，但是當|λ|足夠小時，系統存在極限環。

4 結語

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

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