付志粉，馬建立，李 洋，王 兵

(安徽理工大學(xué)力學(xué)與光電物理學(xué)院，安徽淮南 232001)

0 引言

微分方程的初值問題求解是自然科學(xué)、工程技術(shù)和經(jīng)濟生活等領(lǐng)域經(jīng)常遇到的問題，此類問題的求解有些可以得到解析解，有些則不能，此時利用數(shù)值法求解顯得尤為重要.數(shù)值法求解微分方程的基本思想是離散化，求解方法分兩類，即單步法(如Euler法及Runge-Kutta法等)和多步法(如Adams法)[1-2]，具體的求解過程頗為繁雜，常常需要用到MATLAB或C語言編寫代碼，這給此類方程的求解帶來諸多不便.Simulink是MathsWorks公司基于MATLAB軟件平臺開發(fā)的可用于對動態(tài)系統(tǒng)進行建模、仿真和分析的工具包.在MATLAB環(huán)境中利用Simulink提供的可視化編程界面，用戶利用其內(nèi)置的各種功能化模塊編寫程序，無需進行繁雜的代碼編寫，就可搭建出復(fù)雜的數(shù)學(xué)模型并對其進行仿真求解.鑒于此，本文使用MATLAB/Simulink工具箱對一階微分方程(組)、二階微分方程及高階微分方程的初值問題進行建模、仿真和圖解，以期為該類問題的數(shù)值求解提供參考.

1 微分方程初值問題的MATLAB/Simulink求解過程

使用MATLAB/Simulink工具箱求解微分方程的初值問題可分為以下過程：

(1)啟動Simulink工具箱，選擇相關(guān)模塊.在MATLAB命令窗口中輸入“Simulink”或者單擊MATLAB工具欄中的Simulink圖標(biāo)，即可打開Simulink窗口.單擊該窗口“File”菜單中“new”選項下的“Blank Model”命令，新建一空白模型文件.根據(jù)所求解微分方程左右兩邊各項的特點從模塊庫中選擇相應(yīng)的模塊添加至新建的空白模型文件中.微分方程求解時常用的模塊有積分器Integrator模塊、加法器Add模塊、增益Gain模塊、乘法器Divide/Dot Product模塊及用來顯示求解結(jié)果的示波器Scope模塊等.

(2)搭建微分方程模型.將(1)中所選定的模塊放置好后，按照一定的順序連接起來，構(gòu)成Simulink仿真模型.

(3)設(shè)置仿真參數(shù).在積分器模塊屬性中設(shè)置微分方程的初始條件，單擊工具欄的“Model Configuration Parameters”選項后在打開的頁面中設(shè)置仿真起始和終止時間及具體的方程求解算法類型等.

(4)運行程序，求解方程.單擊工具欄中的程序運行按鈕“Run”，待仿真結(jié)束后單擊示波器模塊，即可得到微分方程解的圖示.

2 MATLAB/Simulink數(shù)值求解微分方程算例

2.1 一階微分方程初值問題的MATLAB/Simulink求解

一階微分方程的求解方法有分離變量法、全微分法、常數(shù)變易法、特征方程法及待定系數(shù)法等[3].有關(guān)一階微分方程初值問題的求解通常是先采用前述方法求出方程的通解，然后將初值代入通解中確定常數(shù)C，最終得到滿足初值條件的一階微分方程的解.這樣的方程求解過程步驟多、計算量大且易出錯.采用MATLAB/Simulink工具箱求解則可以將復(fù)雜、費時的方程求解過程變得簡單、直觀、省時.如求解方程：

y′=y-2x/y

(1)

其中y(0)=1，0≤x≤10的解.

所搭建的Simulink模型框圖如圖1(a)所示.其中積分器Integrator的輸入信號為y′，輸出信號為y.將時鐘信號源Clock模塊提供的信號x和積分器Integrator模塊的輸出信號y分別接入除法器Divide模塊的輸入端，經(jīng)過除法器Divide模塊后的信號即為x/y，此信號經(jīng)過增益Gain模塊擴大-2倍后接入加法器Add模塊的下輸入端口.加法器Add的上輸入端口引入信號y，下輸入端口為信號-2x/y，加法器Add模塊的輸入信號即為y-2x/y，其輸出信號為y′(積分器Integrator模塊的輸入信號).雙擊積分器Integrator模塊在彈出的對話框中設(shè)置初值條件，在工具欄中設(shè)置仿真時間為10 s(此即設(shè)置x的取值范圍).示波器Scope模塊是用來圖示信號y隨x的變化情況的.設(shè)置好參數(shù)后點擊工具欄中的程序運行按鈕“Run”后，點擊示波器即可得到方程的求解結(jié)果，如圖1(b)所示.

圖1 方程(1)的Simulink求解模型框圖及解的結(jié)果Fig.1 Simulink solution model block diagram and solution results of equation (1)

3.2 一階微分方程組初值問題的MATLAB/Simulink求解

1963年，美國氣象學(xué)家Lorenz在研究大氣對流時建立了著名的Lorenz方程，該方程在混沌和可預(yù)報性方面有廣泛的研究.Lorenz方程為一個三元一階微分方程組，其形式為：

(2)

其中α，β，γ為固定值，t為時間，對于給定的初值x1(0)、x2(0)及x3(0)，可利用MATLAB/Simulink工具箱求解該方程.為便于與文獻[4]中的結(jié)果比較，此處取α=8/3，β=10，γ=28，初值x1(0)=x2(0)=0，x3(0)=10-10.所搭建的仿真模型如圖2(a)所示，采用三接口示波器Scope模塊，自上而下依次接入變量x1(t),x2(t)及x3(t),方程的求解結(jié)果即變量的時間響應(yīng)見圖2(b).需要注意的是加法器Add3模塊采用的是3輸入端口的，各輸入端口的“+”“-”可根據(jù)需要設(shè)置，自上而下該加法器的輸入量依次為x3(t)，28x2(t)及x1(t)x2(t).三個To Workspace模塊將積分器輸出的量x1(t),x2(t)及x3(t)分別寫入MATLAB的工作空間中.圖2(c)示出的是在命令窗口中使用plot3函數(shù)繪制的變量x1,x2及x3的三維相空間圖.可以看到采用Simulink求解該方程所得的結(jié)果與文獻[4]中采用代碼編程求解的結(jié)果一致，但采用Simulink求解要簡單、直觀一些.

圖2 Lorenz模型的Simulink求解模型框圖及解的結(jié)果Fig.2 Lorenz model Simulink solution model block diagram and solution results

3.3 二階微分方程初值問題的MATLAB/Simulink求解

1927年荷蘭科學(xué)家Balthasarvan der Pol在研究LC回路電子管振蕩器時建立了一個二階微分方程(簡稱van der Pol方程)，該微分方程在電學(xué)和力學(xué)中有著廣泛的應(yīng)用.van der Pol方程的表達式如下：

y″+μ(y2-1)y′+y=0

(3)

其中μ>0為標(biāo)量參數(shù)，此方程在數(shù)學(xué)上無法精確求解.下面采用MATLAB/Simulink工具箱數(shù)值求解該方程，令方程初始條件為y(0)=-0.2，y′(0)=-0.7，其中μ取2.搭建的Simulink模型框圖如圖3(a)所示，框圖中積分器integrator1模塊的輸入信號為y″，輸出信號(積分器integrator2模塊的輸入信號)為y′，示波器Scope模塊用來顯示y及y′隨時間t的變化，雙綜示波器XY Graph模塊顯示的是y′(t)及y(t)在相平面的運動情況.初值條件及μ的取值可在積分器integrator1、integrator2及增益Gain模塊中修改.仿真時間設(shè)為20 s，仿真結(jié)果如圖3(b)及3(c)所示.與文獻[5]中采用MATLAB語言編程所得結(jié)果比較后發(fā)現(xiàn)，兩者所得結(jié)果一致，但Simulink采用模塊化的編程方式編寫的程序易于理解，可顯示信息的流動，且程序調(diào)試、修改方便，具有直觀、高效的特點.

圖3 van der Pol方程的Simulink求解模型框圖及解的結(jié)果Fig.3 Simulink solution model block diagram of Van der Pol equation and solution results

3.4 高階微分方程初值問題的MATLAB/Simulink求解

對于高階微分方程的求解，通常采用降階法，即利用變換將高階方程化為較低階的方程后進行求解，此過程往往較為復(fù)雜.采用MATLAB/Simulink工具箱數(shù)值求解，可化繁為簡，化抽象為具體.如求解四階非剛性微分方程：

y(4)+2costy″+sinty=7e-t

(4)

其中y(0)=y′(0)=y″(0)=1,y?(0)=0.搭建的Simulink模型框圖如圖4(a)所示，其中積分器integrator1到integrator4模塊的左端分別為輸入信號y(4)，y?，y″及y′，2cost和sint分別由正弦波信號發(fā)生器Sine Wave1及Sine Wave2模塊產(chǎn)生，加法器模塊自上而下輸入的信號分別為7e-t，sinty及2costy″.在各積分器模塊中設(shè)置初值條件，仿真時間設(shè)為100 s，運行后的結(jié)果見圖4(b).可以看到此結(jié)果與文獻[6]采用MATLAB語言編寫代碼求解的結(jié)果一致.相比較而言，采用Simulink求解，數(shù)學(xué)思想清晰，比較符合人們的思維習(xí)慣.

圖4 四階非線性微分方程的Simulink求解模型框圖及解的結(jié)果Fig.4 Fourth-order nonlinear differential equations in Simulink solution model block diagram andsolution results

3 結(jié)論

本文結(jié)合具體算例，利用MATLAB/Simulink對帶有初值條件的一階微分方程(組)、二階微分方程及高階微分方程進行了求解.結(jié)果表明：采用Simulink提供的可視化編程環(huán)境能夠方便地對帶有初值問題的微分方程進行求解和結(jié)果可視化，該方法程序直觀、簡潔，求解速度快，不失為一種求解微分方程初值問題數(shù)值解的有效方法.