數(shù)學(xué)實(shí)驗(yàn)課件：常微分方程實(shí)驗(yàn)

常微分方程實(shí)驗(yàn)微分方程差不多是和微積分同時(shí)先后產(chǎn)生的，蘇格蘭數(shù)學(xué)家耐普爾創(chuàng)立對(duì)數(shù)的時(shí)候，就討論過微分方程的近似解.

牛頓在建立微積分的同時(shí)，對(duì)簡單的微分方程用級(jí)數(shù)來求解．瑞士數(shù)學(xué)家雅各布·貝努利、歐拉、法國數(shù)學(xué)家克雷洛、達(dá)朗貝爾、拉格朗日等人又不斷地研究和豐富了微分方程的理論．常微分方程的形成與發(fā)展是和力學(xué)、天文學(xué)、物理學(xué)，以及其他科學(xué)技術(shù)的發(fā)展密切相關(guān)的.

數(shù)學(xué)的其他分支的新發(fā)展，如復(fù)變函數(shù)、李群、組合拓?fù)鋵W(xué)等，都對(duì)常微分方程的發(fā)展產(chǎn)生了深刻的影響，當(dāng)前計(jì)算機(jī)的發(fā)展更是為常微分方程的應(yīng)用及理論研究提供了非常有力的工具．

牛頓研究天體力學(xué)和機(jī)械動(dòng)力學(xué)的時(shí)候，利用了微分方程這個(gè)工具，從理論上得到了行星運(yùn)動(dòng)規(guī)律.法國天文學(xué)家勒維烈和英國天文學(xué)家亞當(dāng)斯使用微分方程各自計(jì)算出那時(shí)尚未發(fā)現(xiàn)的海王星的位置.這些都使數(shù)學(xué)家更加深信微分方程在認(rèn)識(shí)自然、改造自然方面的巨大力量．本章介紹利用MATLAB求解常微分方程的解析解和數(shù)值解.常微分方程的解析解9.1CONTENTS9.2微分方程的數(shù)值解9.2.1歐拉法9.2.2龍格-庫塔法9.1常微分方程的解析解

定義9.1

含有未知函數(shù)的導(dǎo)數(shù)的方程稱為微分方程.

如果未知函數(shù)是一元函數(shù)，稱為常微分方程.

常微分方程的一般形式為.

由幾個(gè)微分方程聯(lián)立而成的方程組稱為微分方程組．

定義9.2

微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階．若方程中未知函數(shù)及其各階導(dǎo)數(shù)都是一次的，稱為線性常微分方程，一般表示為.

若上式中的系數(shù)

均與

x無關(guān)，稱之為常系數(shù).

二階及二階以上的微分方程稱為高階微分方程．

一階常微分方程與高階微分方程可以互化，已給一個(gè)n階方程設(shè)

，

可將上式化為一階方程組

反過來，在許多情況下，一階微分方程組也可化為高階方程.

所以一階微分方程組與高階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解．

MATLAB中主要用dsolve求符號(hào)解析解，調(diào)用格式如下：

S=dsolve(eqn,cond)

求帶初始條件的微分方程特解，若缺省cond，則求微分方程的通解.利用“diff”和“==”表示微分方程．例如diff(y,x)==y表示dy/dx=y．例9.1

求下列微分方程的解析解（1）（2）（3）解（1）>>symsy(t)ab

>>eqn=diff(y,t)==a*y+b;

>>dsolve(eqn)

ans=

-(b-C1*exp(a*t))/a可知微分方程的通解為

，其中

，C1是任意常數(shù)．（2）>>symsy(x)>>eqn=diff(y,x,2)==sin(2*x)-y;>>Dy=diff(y,x);>>cond=[y(0)==0,Dy(0)==1];>>dsolve(eqn,cond)ans=(5*sin(x))/3-sin(2*x)/3可知微分方程的特解為.（3）>>symsf(t)g(t)>>eqns=[diff(f,t)==f+g,diff(g,t)==g-f];>>cond=[f(0)==1,g(0)==1];>>S=dsolve(eqns,cond)S=包含以下字段的struct:g:[1×1sym]f:[1×1sym]>>S.fans=exp(t)*cos(t)+exp(t)*sin(t)>>S.gans=exp(t)*cos(t)-exp(t)*sin(t)可知微分方程組的特解為.9.2微分方程的數(shù)值解

除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無解析解，應(yīng)用中主要依靠數(shù)值解法．考慮一階常微分方程初值問題其中所謂數(shù)值解法，就是尋求y(x)在一系列離散節(jié)點(diǎn)

上的近似值yk

，

．稱

為步長，通常取為常量h．最簡單的數(shù)值解法是Euler法．9.2.1歐拉法

Euler法的思路極其簡單：在節(jié)點(diǎn)處用差商近似代替導(dǎo)數(shù)

這樣導(dǎo)出計(jì)算公式（稱為Euler格式）

它能求解各種形式的微分方程.Euler法也稱折線法．

Euler方法只有一階精度，改進(jìn)方法有二階Runge-Kutta法、四階Runge-Kutta法、五階Runge-Kutta-Felhberg法和先行多步法等，這些方法可用于解高階常微分方程（組）初值問題.邊值問題采用不同方法，如差分法、有限元法等.

數(shù)值算法的主要缺點(diǎn)是它缺乏物理理解．9.2.2龍格-庫塔法

歐拉公式易于計(jì)算，但精度不高，收斂速度慢．在實(shí)際應(yīng)用中，我們采用龍格-庫塔（Runge-Kutta）方法，基本思想為進(jìn)一步可得在MATLAB中，利用ode23，ode45求微分方程數(shù)值解．

[t,y]=ode23(odefun,tspan,y0)

[t,y]=ode45(odefun,tspan,y0)

求微分方程組y′=f(t,y)從t0到tf的積分，初始條件為y0.其中tspan=[t0tf].解數(shù)組y中的每一行都與列向量t中返回的值相對(duì)應(yīng)．

ode45是最常用的求解微分方程數(shù)值解的命令，采用四階和五階Runge-Kutta算法，是一種自適應(yīng)步長（變步長）的常微分方程數(shù)值解法，對(duì)于剛性方程組不宜采用.ode23與ode45類似，采用二階和三階Runge-Kutta算法，只是精度低一些．ode45是解決數(shù)值解問題的首選方法．例9.2

求解微分方程先求解析解，再求數(shù)值解，并進(jìn)行比較.解>>clear;>>symsy(t);>>eqn=diff(y,t)==-y+t+1;>>cond=y(0)==1;>>dsolve(eqn,cond)ans=t+exp(-t)可得解析解為

．下面再求其數(shù)值解，先編寫M文件fun92.m．%M函數(shù)fun92.mfunctionf=fun92(t,y)f=-y+t+1;再用命令>>clear;t=0:0.1:1;>>y=t+exp(-t);plot(t,y);

%繪制解析解的圖形>>holdon;

%保留已經(jīng)畫好的圖形，如果下面再畫圖，兩個(gè)圖形和并在一起>>[t,y]=ode45(@fun92,[0,1],1);>>plot(t,y,'ro');

%繪制數(shù)值解圖形>>xlabel('t'),ylabel('y')結(jié)果見圖9-1．

由圖9-1所示，解析解和數(shù)值解吻合得很好．圖9-1解析解與數(shù)值解例9.3已知方程當(dāng)

時(shí)，上面方程可化為求上面方程的解析解和數(shù)值解．解先求解析解，>>symsy(t)>>eqn=diff(y,t,2)==9.8*sin(t);>>Dy=diff(y,t);>>cond=[y(0)==15,Dy(0)==0];>>dsolve(eqn,cond)ans=(49*t)/5-(49*sin(t))/5+15可知方程的解析解為

.下面求數(shù)值解．令

可將原方程化為如下方程組建立函數(shù)文件fun93.m如下%M文件fun93.mfunctionf=fun93(t,y)f=[y(2),9.8*sin(y(1))]';

%f向量必須為一列向量運(yùn)行MATLAB代碼>>clear;close;>>[t,y]=ode45(@fun93,[0,10],[15,0]);>>plot(t,y(:,1));

%畫

隨時(shí)間變化圖,y(:2)則表示

的值>>xlabel('t'),ylabel('y1')結(jié)果見圖9-2．

由圖9-2可見，

隨時(shí)間t周期變化．圖9-2數(shù)值解

例9.4

Lotka-Volterra方程，也即捕食者-獵物模型的一對(duì)一階常微分方程

（9-1）

變量

x

和

y

分別計(jì)算獵物和捕食者的數(shù)量.當(dāng)沒有捕食者時(shí)，獵物數(shù)量將增加，當(dāng)獵物匱乏時(shí)，捕食者數(shù)量將減少.使用初始條件x(0)=y(0)=20，使捕食者和獵物的數(shù)量相等．求當(dāng)α=0.01和β=0.02時(shí)方程的解．解在MATLAB中，兩個(gè)變量x和y可以表示為向量y中的兩個(gè)值．同樣，導(dǎo)數(shù)是向量yp中的兩個(gè)值.當(dāng)α=0.01和β=0.02時(shí)，方程組（9-1）可以表示為：

yp(1)=(1-alpha*y(2))*y(1)

yp(2)=(-1+beta*y(1))*y(2)MATLAB中的函數(shù)文件lotka.m來表示Lotka-Volterra方程，typelotkafunctionyp=lotka(t,y)%LOTKALotka-Volterrapredator-preymodel.%Copyright1984-2014TheMathWorks,Inc.yp=diag([1-.01*y(2),-1+.02*y(1)])*y;首先使用ode23在區(qū)間0<t<15中求解lotka中定義的微分方程．t0=0;tfinal=15;y0=[20;20];[t,y]=ode23(@lotka,[t0tfinal],y0);繪制兩個(gè)種群數(shù)量對(duì)時(shí)間的圖，見圖9-3a．plot(t,y(:,1),'--',t,y(:,2))title('捕食者-獵物模型')xlabel('時(shí)間')ylabel('數(shù)量')legend('獵物','捕食者','Location','North')再使用

ode45

求解該方程組，得到相軌線圖9-3b.[T,Y]=ode45(@lotka,[t0tfinal],y0);plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'--');title('相軌線')legend('ode23','ode45')

圖9-3捕食者-獵物模型a)捕食者—獵物模型b）相軌線

圖9.3a中實(shí)線表示的是微分方程的解，虛線表示的是解曲線的導(dǎo)數(shù)所描繪的曲線.可以看出，獵物數(shù)量減少時(shí)，捕食者的數(shù)量也在減少，捕食者的種群數(shù)量會(huì)隨著獵物種群數(shù)量的增加而不斷增加．獵物的種群數(shù)量增加時(shí)，捕食者數(shù)量也在增加，但是當(dāng)捕食者達(dá)到一定的程度后，獵物又在不斷減少.即獵物種群數(shù)量增加→捕食者種群數(shù)量增加→獵物種群數(shù)量減少→捕食者種群數(shù)量減少→獵物種群數(shù)量增加→再次循環(huán)．這種變化趨勢(shì)反映了生態(tài)系統(tǒng)中普遍存在的負(fù)反饋調(diào)節(jié)．

圖9.3b中，使用不同的數(shù)值方法求解微分方程會(huì)產(chǎn)生略微不同的答案.可以看出利用ode45函數(shù)得到的圖形是平滑的，要比ode23函數(shù)得到的結(jié)果精度更高．

例9.5（羅倫茲吸引子與“蝴蝶效應(yīng)”）吸引子在1963年由美國麻省理工學(xué)院的氣象學(xué)