屈小妹

(湖北師范大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，湖北 黃石 435002)

0 引言

微分方程的研究來源廣泛，歷史久遠(yuǎn)。含有未知函數(shù)導(dǎo)數(shù)的方程，稱之為微分方程，其一般形式為：

f(x,y,y′,y″,…,y(n))=0

其中，未知函數(shù)是一元函數(shù)的，叫常微分方程；未知函數(shù)是多元函數(shù)的叫做偏微分方程。微分方程在幾何學(xué)、力學(xué)、經(jīng)濟(jì)學(xué)、物理學(xué)、生物學(xué)、化學(xué)等領(lǐng)域有非常重要的應(yīng)用。大部分微分方程很難求出十分精確的解，而只能得到近似解。因而，采用數(shù)值方法獲取微分方程的近似解受到廣泛關(guān)注，常見的求數(shù)值解的方法有歐拉折線法、亞當(dāng)姆斯法、龍格-庫塔法等[1]。

含有未知函數(shù)的導(dǎo)數(shù)高于一階的微分方程叫高階微分方程。n階線性微分方程的一般形式為：

x(n)+a1(t)x(n-1)+…+a(n-1)(t)x′+an(t)x=f(t)

若f(t)=0,稱為齊次方程;若f(t)≠0,稱為非齊次方程。特別的，x″+p(t)x′+q(t)=f(t)稱為二階線性微分方程。一般來說，高階微分方程的求解比較復(fù)雜，高階線性微分方程的解析解研究已積累了一定成果[2,3]，其求解方法有：常數(shù)變異法、特征根法、比較系數(shù)法、拉普拉斯變換法等[4]。針對(duì)高階非線性微分方程，其正解的適定性研究得到較多關(guān)注[5,6]。對(duì)于更加廣泛的、非線性的一般高階微分方程，從理論方面進(jìn)行研究還存在相當(dāng)大的困難。人們常采用降階法，即利用變換將高階方程化為較低階的方程。已知一個(gè)n階微分方程

y(n)=f(x,y,y′,y″,…,y(n-1))

(1)

設(shè)y1=y，y2=y′,…,yn=y(n-1),則解(1)式相當(dāng)于求解下面n個(gè)一階微分方程組

(2)

常微分方程的所有數(shù)值算法都是以一階微分方程組為求解對(duì)象，對(duì)(2)可應(yīng)用大量數(shù)值方法求解。

MATLAB軟件具有強(qiáng)大的數(shù)值矩陣計(jì)算能力、大量計(jì)算算法和繪圖能力，可以幫助人們脫離復(fù)雜的計(jì)算困境，也能用強(qiáng)大的圖形功能對(duì)數(shù)值計(jì)算結(jié)果進(jìn)行顯示[7,8]?，F(xiàn)今對(duì)于高階微分方程應(yīng)用MATLAB軟件求解的文章還比較少，對(duì)于高階延遲微分方程如何應(yīng)用MATLAB求解和程序?qū)崿F(xiàn)的研究文章更是少之又少。本文研究應(yīng)用MATLAB求解高階微分方程初值問題，通過多個(gè)算例分析高階線性微分方程、非剛性、剛性van der Pol方程及高階延遲微分方程的算法求解及程序?qū)崿F(xiàn)。

1 應(yīng)用MATLAB求解高階常微分方程

1.1 二階線性常微分方程初值問題

對(duì)于帶初邊值條件的二階微分方程，若方程存在解析解，MATLAB命令格式為：

dsolve(‘equation’, ‘con1, con2,…’,’variable’),其中equation是待解的微分方程，con1,con2,…為其初邊值條件，variable為自定義的變量名。

解：1)求解析解:

輸入命令：z=dsolve(‘D2z+2*Dz+z=0’, ‘z(0)=4,Dz(0)=-2’,’t’)

得到結(jié)果z=4*exp(-t)+6*exp(-t)*t.

2)計(jì)算數(shù)值解，并繪制圖形。

function dz=fun1(t,z)

dz=[z(2);-2*z(2)-z(1)];

end

輸入命令：[t,z]=ode45(‘fun1’,[0,10],[4;-2]);

plot(t,z(:,1),'-',t,z(:,2),’*’);

legend(‘z_1’,’Dz_1’,’Location’,’NorthEast’);

運(yùn)行后求得解函數(shù)z=z1(t)圖形如圖1所示。

圖1 例1運(yùn)行結(jié)果

1.2 van der Pol方程

van der Pol方程[9]是荷蘭科學(xué)家Balthasarvan der Pol于1927年研究極限環(huán)獲得的，它在電學(xué)和力學(xué)中有廣泛的應(yīng)用，此方程為二階微分方程。當(dāng)μ為小量時(shí)，研究人員可采用頻數(shù)展開法、諧波平衡法及伽遼金法等求出定常周期解。μ不為小量時(shí)，許多方法費(fèi)時(shí)、費(fèi)力，甚至難以反映解的所有特征，因而采用數(shù)值方法模擬是一種行之有效的方法[9～11]。

① 當(dāng)μ=0.8時(shí)，生成的方程組為非剛性方程組，可以使用非剛性求解器ode45(中階方法)、ode23(低階方法)、ode113(變階方法)。求解過程如下：

1)建立名為van1.m的M函數(shù)文件

function dydt=van1(t,y)

dydt=[y(2); 0.8*(1-y(1)^2)*y(2)-y(1)];

2)調(diào)用van1.m求解，在命令窗口輸入

tspan=[0,20]; %設(shè)置仿真時(shí)間20秒

x0=[2;0]; %設(shè)置仿真初始值x(0)=2, x’(0)=0

[t,y]=ode45(@van1,tspan,x0);

plot(t,y(:,1),'-o');

xlabel(' Time t');

ylabel(' Solution y');

運(yùn)行后求得解函數(shù)y1(t)及其圖形如圖2所示。

圖2 例2運(yùn)行結(jié)果(μ=0.8)

②當(dāng)μ=3000時(shí)，生成的方程組為剛性方程組，解會(huì)在較長(zhǎng)的時(shí)間段中顯示振蕩，可以使用ode15s(變階方法)、ode23s(低階方法)、ode23t(梯形法則)、ode23tb(梯形法則+向后差分公式)等剛性求解器。

求解過程如下：

1)建立名為van2.m的M函數(shù)文件

function dydt=van2(t,y)

dydt=[y(2); 3000*(1-y(1)^2)*y(2)-y(1)];

2)調(diào)用van2.m求解，在命令窗口輸入

tspan=[0,9000]; %設(shè)置仿真時(shí)間9000秒

x0=[2;0]; %設(shè)置仿真初始值

[t,y]=ode23s(@van2,tspan,x0);

plot(t,y(:,1),'-o');

xlabel(' Time t');

ylabel(' Solution y');

運(yùn)行后求得解函數(shù)y1(t)及其圖形如圖3所示。

圖3 例2運(yùn)行結(jié)果(μ=3000)

1.3 四階非剛性微分方程

解：1)令y=y1,y′=y2,y(2)=y3,y(3)=y4,

2)建立名為fun2.m的M函數(shù)文件

function dydt=fun2(t,y)

dydt=[y(2);y(3);y(4);7*exp(-t)-sin(t)*y(1)-2*cos(t)*y(3)];

end

3)調(diào)用fun2.m求解，在命令窗口輸入

tspan=[0,10];

y0=[1;1;1;0];

[t,y]=ode45(@fun2,tspan,y0);

plot(t,y(:,1),'-o');

運(yùn)行后求得解函數(shù)y=y1(t)及其圖形如圖4所示。

圖4 例3運(yùn)行結(jié)果

2 MATLAB求解高階延遲微分方程

延遲微分方程(DDE)是當(dāng)前時(shí)間的解與過去時(shí)間的解相關(guān)的微分方程，其理論解較之常微分方程更難獲得，因而對(duì)其數(shù)值解的研究具有重要的科學(xué)意義[1]。具有常延遲的微分方程形式如下：

y′(t)=f(t,y(t),y(t-τ1)，…，y(t-τk))

這里，t為自變量，y為因變量，而y′表示y關(guān)于t的一階導(dǎo)數(shù)。延遲τ1,…,τk是正常數(shù)。MATLAB常使用dde23函數(shù)求解常延遲微分方程(組)。

解：1)創(chuàng)建向量定義兩個(gè)延遲：

Lags=[1,0.2]; %τ1=1,τ2=0.2

2)建立名為dde2.m的M函數(shù)文件，定義該方程組：

function dydt=dde2(t,y,z)

ylag1=z(:,1); %逼近y1(t-1)

ylag2=z(:,2); %逼近y2(t-0.2)

dydt=[-0.8*ylag1(1)+0.1*ylag2(2);-0.5*ylag1(1);0.05*y(2)-ylag1(1)];

end

3)編寫歷史解代碼：

function s=his2(t)

s=ones(3,1);

end

4)求解方程及繪圖：

tspan=[0,100];

sol=dde23(@dde2,Lags,@his2,tspan);

plot(sol.x,sol.y,′-o′)

xlabel(′t′);

ylabel(′y′);

legend(′y_1′,′y_2′,′y_3′,′Location′,′NorthEast′);

運(yùn)行后繪制三個(gè)解分量對(duì)時(shí)間的圖見圖5.

圖5 例4運(yùn)行結(jié)果

3 結(jié)語

高階常微分方程和延遲微分方程數(shù)值求解往往需要先轉(zhuǎn)化為一階微分方程組，再選擇合適的算法調(diào)用恰當(dāng)?shù)腗ATLAB求解器進(jìn)行求解。對(duì)于非線性系統(tǒng)仿真利用M函數(shù)可以直接給出動(dòng)態(tài)系統(tǒng)的導(dǎo)數(shù)描述，較為方便。應(yīng)用MATLAB進(jìn)行編程求解高階微分方程，程序直觀、簡(jiǎn)潔，求解時(shí)間快，易學(xué)易懂。