數(shù)學實驗課件 第9章9.2

9.2微分方程的數(shù)值解

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

上的近似值yk

，

．稱

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

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

這樣導出計算公式（稱為Euler格式）

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

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

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

歐拉公式易于計算，但精度不高，收斂速度慢．在實際應用中，我們采用龍格-庫塔（Runge-Kutta）方法，基本思想為進一步可得在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中返回的值相對應．

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

求解微分方程先求解析解，再求數(shù)值解，并進行比較.解>>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)畫好的圖形，如果下面再畫圖，兩個圖形和并在一起>>[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已知方程當

時，上面方程可化為求上面方程的解析解和數(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向量必須為一列向量運行MATLAB代碼>>clear;close;>>[t,y]=ode45(@fun93,[0,10],[15,0]);>>plot(t,y(:,1));

%畫

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

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

由圖9-2可見，

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

例9.4

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

（9-1）

變量

x

和

y

分別計算獵物和捕食者的數(shù)量.當沒有捕食者時，獵物數(shù)量將增加，當獵物匱乏時，捕食者數(shù)量將減少.使用初始條件x(0)=y(0)=20，使捕食者和獵物的數(shù)量相等．求當α=0.01和β=0.02時方程的解．解在MATLAB中，兩個變量x和y可以表示為向量y中的兩個值．同樣，導數(shù)是向量yp中的兩個值.當α=0.01和β=0.02時，方程組（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);繪制兩個種群數(shù)量對時間的圖，見圖9-3a．plot(t,y(:,1),'--',t,y(:,2))title('捕食者-獵物模型')xlabel('時間')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ù)所描繪的曲線.可以看出，獵物數(shù)量減少時，捕食者的數(shù)量也在減少，捕食者的種群數(shù)量會隨著獵物種群數(shù)量的增加而不斷增加．獵物的種群數(shù)量增加時，捕食者數(shù)量也在增加，但是當捕食者達到一定的程度后，獵物又在不斷減少.即獵物種群數(shù)量增加→捕食者種群數(shù)量增加→獵物種群數(shù)量減少→捕食者種群數(shù)量減少→獵物種群數(shù)量增加→再次循環(huán)．這種變化趨勢反映了生態(tài)系統(tǒng)中普遍存在的負反饋調(diào)節(jié)．

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

例9.5（羅倫茲吸引子與“蝴蝶效應”）吸引子在1963年由美國麻省理工學院的氣象學家羅倫茲（E.N.Lorenz）發(fā)現(xiàn).羅倫茲教授在研究天氣的不可預測性時，通過簡化方程，獲得了具有三個自由度的系統(tǒng)．在計算機上用他所建立的微分方程模擬氣候變化，意外地發(fā)現(xiàn)，初始條件的極微小差別可以引起模擬結(jié)果的巨大變化，這表明天氣過程以及描述它們的非線性方程時如此的不穩(wěn)定，以至巴西熱帶雨林的一只蝴蝶偶然拍動一下翅膀，幾星期后可以在美國德克薩斯州引起一場龍卷風，這就是“蝴蝶效應”．

羅倫茲根據(jù)牛頓定律建立的溫度、壓強和風速之間的微分方程組為給定初值條件求當

時，微分方程組的解．解當

時，得到微分方程組：將三個方程的右端函數(shù)寫成向量形式：新建函數(shù)文件flo.m，命令如下：functionz=fl