MATLAB在常微分方程上的應(yīng)用

趙亞琪 范進軍

( 山東師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院, 250358, 濟南 )

1 引 言

微分方程奇解的存在性、奇點的類型、極限環(huán)的存在性和零解的穩(wěn)定性的判斷是常微分方程理論中的重要內(nèi)容[1-3], MATLAB具有強大的繪圖功能且在常微分方程的抽象問題具體化方面有極大優(yōu)勢[4,5]. 本文從理論角度對微分方程的上述四類問題進行分析, 給出了微分方程奇點、極限環(huán)的定義及奇解的存在性、奇點類型、極限環(huán)的存在性與零解穩(wěn)定性的判定方法, 鑒于理論判斷較抽象, 考慮應(yīng)用MATLAB的繪圖功能探討解決四類問題的新方法, 通過結(jié)合圖形與理論判斷達(dá)到由繁到簡和形象直觀解決問題的目的.

2 預(yù)備知識

2.1奇解定義設(shè)一階微分方程

(1)

有一特解Λ:y=f(x)(x∈I).

如果對每一點P∈A, 在P點的某個鄰域內(nèi)方程(1)有一個不同于Λ的解在P點與Λ相切, 則稱Λ為方程(1)的奇解.

注1當(dāng)常微分方程有兩組積分曲線族, 且有一特解為其中一個積分曲線族的包絡(luò), 則此包絡(luò)是常微分方程的奇解.

2.3一階線性微分方程組奇點分類的直觀判斷依據(jù)給定二維常系數(shù)線性自治系統(tǒng)：

令p=-(a+d)=-tr(A),q=ad-bc=det(A).當(dāng)q≠0時, 一階線性微分方程組奇點分類的直觀判斷如圖1所示.

圖1 奇點類型圖(圖形中坐標(biāo)系是變換為標(biāo)準(zhǔn)方程后的坐標(biāo)系x′O y′)

2.4極限環(huán)概念設(shè)系統(tǒng)

(2)

具有閉軌線C.假如在C充分小鄰域中, 除C之外, 軌線全不是閉軌線, 且這些非閉軌線當(dāng)t→+或t→-時趨近于閉軌線C, 則稱閉軌線C是 (2)的一個極限環(huán).

2.5零解穩(wěn)定性的概念對于動力系統(tǒng)

(3)

其中:

如果對?ε>0,t0≥0,?δ(ε,t0)>0,?x0∈G,只要式(3)滿足初值x(t0)=x0的解x(t;t0,x0)都在t≥t0時有定義且‖x(t;t0,x0)‖<ε, 則稱方程(3)的零解是穩(wěn)定的.

3 MATLAB在常微分方程中的具體應(yīng)用

3.1運用MATLAB判斷奇解的存在性

例1考慮微分方程

(4)

易知方程通解為

(5)

y=x-C2,

(6)

顯然

y=0

(7)

是方程的一個解.

運用MATLAB繪制方程(4)的第一個積分曲線族(5)(圖2), 程序如下：

clear; clc; hold on;

x=linspace(-400,400,6 000);y=x-x;plot(x,y,'-k');

y1=0;

fori=-100:17:100

y1=1/27*(x-i).^3;

plot(x,y1,'-r');

end

hold off

由圖2可知曲線族有包絡(luò), 故方程(4)有奇解. 同時, 觀察可知方程(7)的曲線為方程(4)第一個積分曲線族(5)的包絡(luò), 故y=0是奇解.

運用MATLAB繪制方程(4) 第二個積分曲線族(6)(圖3), 程序如下：

clear; clc; hold on;

x=linspace(-400,400,6 000);y2=0;

fori=-100:17:100

y2=x-i;

plot(x,y2,'-r');

end

hold off

圖2 積分曲線族(5)的局部圖像

圖3 積分曲線族(6)的局部圖像

由圖3可知方程(4)第二個積分曲線族(6)無奇解.

1) 理論分析驗證.由(5)與(6)得到方程(4)的通積分

它的C-判別式為

由此得到

Λ:x=C,y=0(-∞

2) 在判斷奇解中MATLAB優(yōu)缺點的分析.

優(yōu)點：運用MATLAB繪制常微分方程的積分曲線族能夠快速并直觀判斷出奇解的存在性.

缺點：運用MATLAB繪圖是由x與y遍歷區(qū)間數(shù)值得到的, 繪制出的圖形有誤差, 對于一些復(fù)雜的曲線族會容易出現(xiàn)判斷失誤的情況.

3.2運用MATLAB繪圖判斷奇點類型

例2討論如下系統(tǒng)奇點類型

(8)

易知(0,0)和(-1,0)為動力系統(tǒng)(8)的奇點, 運用MATLAB繪制奇點附近的相圖(圖4), 程序如下：

clear; clc;x0=-1.3:0.01:0.3;y0=-1:0.01:1;

[x,y]=meshgrid(x0,y0);

dx=y;dy=-x-y-x.^2-y.^2;d=sqrt(x.^2+y.^2);

u=dx./d;v=dy./d; figure; streamslice(x,y,u,v);

hold on; plot([-1.3,0.3],[0,0],'k','linewidth',0.5);

hold on; plot([0,0],[-1,1],'k','linewidth',0.5);

hold on; plot([-1,-1],[-1,1],'k','linewidth',0.5)

由圖4直觀觀察, 易看出(-1,0)和(0,0)分別為系統(tǒng)(8)的鞍點和穩(wěn)定焦點.

例3判斷如下系統(tǒng)的奇點類型

(9)

運用MATLAB繪制圖形(圖5), 程序如下：

clear; clc; hold on;

x0=-0.3:0.01:0.3;y0=-0.3:0.01:0.3;

[x,y]=meshgrid(x0,y0);

dx=2*x+3*y;dy=2*x-3*y; d=sqrt(x.^2+y.^2);

u=dx./d;v=dy./d; figure; streamslice(x,y,u,v);

hold on; plot([-0.3,0.3],[0,0],'k','linewidth',0.5);

hold on; plot([0,0],[-0.3,0.3],'k','linewidth',0.5)

觀察圖5, 易看出(0,0)為系統(tǒng)(9)的鞍點.

圖4 動力系統(tǒng)(8)在兩奇點附近的相圖

圖5 動力系統(tǒng)(9)在奇點(0,0)附近的相圖

1) 理論分析驗證.在例2中, 由

得系統(tǒng)(8)有兩個奇點A(-1,0)和B(0,0).

對于奇點A(-1,0), 作可逆變量代換x'=x+1,y'=y, 則系統(tǒng)(8)可化為

(10)

系統(tǒng)(10)對應(yīng)的線性近似系統(tǒng)的系數(shù)矩陣為

高次項為

φ(x',y')=0，ψ(x',y')=-(x')2-(y')2.

計算:

p=-trD=-(0-1)=1>0,q=detD=-1<0,p2-4q=1+4=5>0，

從而系統(tǒng)(10)的奇點(0,0)為鞍點.

由于可逆線性變換不改變系統(tǒng)奇點的類型, 所以奇點A(-1,0)也是系統(tǒng)(8)的鞍點.

對于奇點B(0,0),系統(tǒng)(8)的線性近似系統(tǒng)為

(11)

相應(yīng)線性近似系統(tǒng)(11)的系數(shù)矩陣為

高次項為

φ(x,y)=0，ψ(x,y)=-x2-y2.

計算:

p=-trD=-(0-1)=1>0,q=detD=1>0,p2-4q=1-4=-3<0.

所以奇點B(0,0)是系統(tǒng)(8)也是系統(tǒng)(8)的穩(wěn)定焦點.

綜上所述, 運用MATLAB對例2中兩奇點類型的判斷正確. 類似地, 可以驗證例3對奇點類型的判斷正確.

2) 在判斷奇點類型中MATLAB優(yōu)缺點的分析.在判斷奇點類型中, 運用MATLAB可以通過軌線分布的特點而快速分辨出系統(tǒng)中奇點的類型, 但是在通過圖形直接確定奇點坐標(biāo)方面容易出現(xiàn)偏差, 難度也較大.所以, 運用MATLAB軟件繪圖對于判斷奇點的存在性有很大優(yōu)勢, 而對于奇點坐標(biāo)的確定運用理論判斷更方便快捷.

3.3運用MATLAB繪圖判斷動力系統(tǒng)是否存在極限環(huán)

例4判斷方程組(動力系統(tǒng))

(12)

是否有極限環(huán)并確定極限環(huán)的位置.

運用MATLAB繪制方程組(12)的軌線分布(圖6), 程序如下：

clear;x0=-3:0.1:3;y0=-3:0.1:3; [x,y]=meshgrid(x0,y0);

dx=y+x.*sqrt(x.^2+y.^2).*(1-(x.^2+y.^2)).^2;

dy=-x+y.*sqrt(x.^2+y.^2).*(1-(x.^2+y.^2)).^2;

d=sqrt(x.^2+y.^2);u=dx./d;v=dy./d; figure;

streamslice(x,y,u,v); hold on;t=0:0.01:2*pi;

x=@(t)cos(t);y=@(t)sin(t); fplot(x,y,'-r','linewidth',1);

hold off

圖6 動力系統(tǒng)(12)的局部相圖

1) 運用理論知識驗證分析.對于例4, 考察函數(shù)族V(x,y)=x2+y2=C(其中C為參數(shù)). 求V沿著方程組的全導(dǎo)數(shù)

G為由Γ1和Γ2所圍成的環(huán)域, 則G滿足Poincare-Bendixon環(huán)域定理的條件.由于r1、r2可與1任意接近, 所以單位圓x2+y2=1就是(12)的一個極限環(huán). 可見運用MATLAB繪圖對(12)極限環(huán)存在性的判斷正確.

2) 在動力系統(tǒng)是否存在極限環(huán)的判斷中MATLAB的優(yōu)缺點分析.在動力系統(tǒng)是否存在極限環(huán)的判斷中, 運用MATLAB繪圖可以清晰地觀察到方程組(動力系統(tǒng))中x與y隨時間t的增加而形成的軌線分布情況, 結(jié)合極限環(huán)的定義容易辨別極限環(huán)的存在性.但鑒于直觀辨別的誤差性較大, 不易得到極限環(huán)的方程式. 所以MATLAB在動力系統(tǒng)極限環(huán)存在性的判斷中有較大優(yōu)勢, 可以簡化問題, 但對于極限環(huán)表達(dá)式的確定需要運用理論分析得出.

3.4運用MATLAB判斷動力系統(tǒng)零解是否具有穩(wěn)定性

例5考慮如下自治系統(tǒng)零解的穩(wěn)定性

(13)

運用MATLAB繪制圖形(圖7、圖8)：

首先, 定義函數(shù)

function [x,y,wx,wy] =H(gxy,hx,hy,nx)

if nargin<4,nx=20;

elsenx=nx*20;

end

ny=floor(nx*(hy(2)-hy(1))/(hx(2)-hx(1)));

bx=(hx(2)-hx(1))/nx;by=(hy(2)-hy(1))/ny;x=hx(1):bx:hx(2);

y=hy(1):by:hy(2); [x,y]=meshgrid(x,y);kk=gxy(x,y);

wx=1./(1+kk.^2);wy=wx.*kk;

end

繪制方程組(13)的方向場(圖7), 程序如下：

clc; clear; clf;gxy=@(x,y)(-y./x);

hx=[-0.003,0.003];hy=[-0.003,0.003];

[x,y,wx,wy]=H(gxy,hx,hy); quiver(x,y,wx,wy,'color','b')

繪制細(xì)節(jié)圖(圖8), 程序如下：

clear; clc;x0=-0.000 004:0.000 000 1:0.000 004;y0=-0.000 004:0.000 000 1:0.000 004;

[x,y]=meshgrid(x0,y0);dx=-y;dy=x; d=sqrt(x.^2+y.^2);

u=dx./d;v=dy./d;

figure; streamslice(x,y,u,v)

由圖8及零解穩(wěn)定性定義可知, 方程組 (13)零解是一致穩(wěn)定的, 但不是漸進穩(wěn)定的.

圖7 方程組(13)在(0,0)附近的方向場

圖8 方程組(13)在(0,0)附近方向場細(xì)節(jié)圖

例6判斷如下系統(tǒng)零解的穩(wěn)定性

(14)

運用MATLAB繪制(14)的軌線圖(圖9), 程序如下：

clear; clc;x0=-0.004:0.000 1:0.004;y0=-0.004:0.000 1:0.004;

[x,y]=meshgrid(x0,y0); dx=-y+x.*(x.^2+y.^2-1); dy=x+y.*(x.^2+y.^2-1);

d=sqrt(x.^2+y.^2); u=dx./d;v=dy./d; figure; streamslice(x,y,u,v)

圖9 動力系統(tǒng)(14)在(0,0)點附近的軌線圖

由圖9及零解穩(wěn)定性定義可知, 方程組 (13)的零解是一致漸進穩(wěn)定的.

2) 在動力系統(tǒng)零解是否具有穩(wěn)定性的判斷中MATLAB的優(yōu)缺點分析.MATLAB在動力系統(tǒng)是否具有穩(wěn)定性的判斷中不僅可以判斷出其是否具有穩(wěn)定性, 還可以判斷出其一致穩(wěn)定性和漸進性.所以, 在動力系統(tǒng)是否具有穩(wěn)定性的判斷中, 相較于理論判斷運用MATLAB繪圖判斷具有很大的優(yōu)越性, 能夠快速并準(zhǔn)確得出結(jié)論.

4 結(jié) 語

MATLAB在常微分方程上的運用是多方面的, 也是快速解決問題的一種方式. 運用MATLAB繪圖功能在解決判斷微分方程奇解的存在性、奇點的類型、極限環(huán)的存在性、零解的穩(wěn)定性這四類問題上有較大作用, 不僅快捷、簡便, 而且正確度高, 在不需要太過考慮誤差的情況下也能夠判斷無誤. 未來, 我們對常微分方程方面的研究可以多借助MATLAB軟件, 通過繪制方程的軌線圖或相圖, 有助于使抽象的方程(組)或動力系統(tǒng)更加形象、直觀, 易于理解.