精通MATLAB科學計算（第3版）(王正林)12-3r

第章常微分方程求解

常微分方程作為微分方程的基本類型之一，在自然界與工程界有很廣泛的應用。很多

問題的數(shù)學表述都可以歸結(jié)為常微分方程的定解問題。很多偏微分方程問題，也可以化為

常微分方程問題來近似求解。MATLAB中提供了求解函數(shù)和求解器，可以實現(xiàn)常微分方

程的解析求解和數(shù)值求解，而且還可以通過編程實現(xiàn)一些常見的數(shù)值解法。

通過本章學習，讀者能夠熟練使用MATLAB的求解函數(shù)和求解器，以及通過編程進

行常微分方程的求解。

12.1MATLAB中的求解函數(shù)

一般微分方程式描述系統(tǒng)內(nèi)部變量的變化率如何受系統(tǒng)內(nèi)部變量和外部激勵，如輸入

的影響。當常微分方程式能夠解析求解時，可用MATLAB的符號工具箱中的功能找到精

確解;在常微分方程難以獲得解析解的情況下使用MATLAB的常微分方程求解器solver,

可以方便地在數(shù)值上求解。

12.1.1符號解函數(shù)dsolve

一階常微分方程式(first-orderOrdinaryDifferentialEquation,ODE)可寫為：

y=g(x,y)

其中x為獨立變量，而)，是x的函數(shù)。它的解是：

y=.f(x,y)

該解可以滿足)/=f=g(xj),在初始條件.i，(xo)=yo下，可以得到唯一解。

MATLAB常微分方程符號解的語法是：

dsolve('equation*,1condition1)

其中，equation代表常微分方程式即./=gQy),且須以Dy代表一階微分項y,D2y

代表二階微分項y",condition則為初始條件。

函數(shù)dsolve用來解符號常微分方程、方程組，如果沒有初始條件，則求出通解，如果

有初始條件，則求出特解。

第12章常微分方程求解

dsolve的調(diào)用格式如下：

(1)dsolve('equation')

給出微分方程的解析解，表示為，的函數(shù)；

(2)dsolve('equation'「condition')

給出微分方程初值問題的解，表示為/的函數(shù)；

(3)dsolveC'equation','vf)

給出微分方程的解析解，表示為v的函數(shù)；

(4)dsolve(*equation','condition*,V)

給出微分方程初值問題的解，表示為v的函數(shù)。

【例12-1】常微分方程符號解求解實例1。計算微分方程或+3孫=配-『的通解。

dx

>>dsolve('Dy+3*x*y=x*exp(-xA2)1)

輸出為：

ans=(l/3*exp(-x*(x-3*-))+C1)*exp(-3*x*t)

t))+C1)*?士―3***、?&)?

由于系統(tǒng)默認的自變量是/,顯然系統(tǒng)把x當作常數(shù)，把y當作/的函數(shù)求解.輸入命令:

A

>>dsolve('Dy+3*x*y=x*exp(-x2)'z*x')

輸出為：

ans=exp(-x、)+exp(-3/2*xA2)*C1

exp(-xA2)+oxp(3/2*M2)*C1

_32

可知通解為y=*G，其中G為常數(shù)。

【例12-2】常微分方程符號解求解實例2。計算微分方程中'+2y-/=0在初始條

件川”1=2e下的特解。

>>dsolve(*x*Dy+2*y-exp(x)=0',*y(1)=2*exp(1)',*x')

輸出為：

ans=(一xp(x)*x-exp(x)+2*xxp(1))/x八2

(oxp(x)*x-oxp(x)+2*€派斷(44-)-^^笠

可知特解為y=+2e。

X

【例12-3】常微分方程符號解求解實例3。求爐+2丁+產(chǎn)=0的通解。

>>dsolve(*D2y+2*Dy+exp(x)=0',1x*)

Y???281

精通MATLAB科學計算(第2忡-------------

輸出為：

ans=-l/3*exp(x)-l/2*exp(-2*x)*C1+C2

可知通解為"」,-LN*。。？，其中G,C2為常數(shù)。

32

282????

第12章常微分方程求解

12.1.2求解器solver

在求常微分方程數(shù)值解方面，MATLAB具有豐富的函數(shù)，將其統(tǒng)稱為solver,其一般

格式為：

[TfY]=solver(odefun,tspan,yO)

該函數(shù)表示在區(qū)間tspan=Do,5上,用初始條件川求解顯式常微分方程V=

odefun為顯式常微分方程V=f(t,y)中的/(f,y),tspan為求解區(qū)間，要獲得問題在其

他指定點2…上的解,則令他劭=汨/1,/2,…0(要求"單調(diào))，y0為初始條件。

solver為命令ode45、ode23、odel13,ode15s,ode23s、ode23t、ode23tb之一，其中

ode45、ode23、odel13屬于非剛性ODE類型,這些命令的特點如表12-1所示;ode15s,

ode23s、ode23t、ode23tb屬于剛性ODE類型，這些命令的特點如表12-2所示。

表12-1非剛性ODE求解命令

求解器solver功能說明

ode45一步算法；4,5階龍格■■庫塔方程；累計截斷誤差達(&03大部分場合的首選算法

ode23一步算法；2,3階龍格-庫塔方程；累計截斷誤差達(加個使用于精度較低的情形

odel13多步法；Adams算法；高低精度均可到IO。~I。/,計算時間比ode45短

表12-2剛性ODE求解命令

求解器solver功能說明

ode23t梯形算法適度剛性情形

ode15s多步法；Gear's反向數(shù)值微分；精度中等若ode45失效時，可嘗試使用

ode23s一步法；2階Rosebrock算法；低精度當精度較低時，計算時間比odel5s短

ode23tb梯形算法；低精度當精度較低時，計算時間比odel5s短

函數(shù)ode45與ode23是常使用的求解方法，函數(shù)ode45的使用與ode23完全一樣。兩

個函數(shù)的差別在于必須與所用的內(nèi)部算法相關。兩個函數(shù)都運用了基本的龍格-庫塔

(Runge-Kutta)數(shù)值積分法的變形。

ode23運用一^組合的2/3階龍格-庫塔-芬爾格(Runge-Kutta-Fehlerg)算法，而ode45運

用組合的4/5階龍格-庫塔-芬爾格算法。

通常，ode45可取較多的時間步。所以，要保持與ode23相同誤差時，在訪和夕之間

可取較少的時間步。然而，在同一時間內(nèi)，ode23每時間步至少調(diào)用3次，而ode45每時

間步至少調(diào)用6次。

正如使用高階多項式內(nèi)插常常得不到最好的結(jié)果一樣，ode45也不總是比。de23好。

如果。de45產(chǎn)生的結(jié)果，對作圖間隔太大，則必須在更細的時間區(qū)間內(nèi)對數(shù)據(jù)進行內(nèi)插，

比如用函數(shù)interpl。這個附加時間點會使ode23更有效。作為一條普遍規(guī)則，在所計算的

1???283

精通MATLAB科學計算（第2版i-----------------------------------------

導數(shù)中，如有重復的不連續(xù)點，為保持精度致使高階算法減少時間步長，這時低階算法更

有效。

采用求解器solver求解ODE的基本過程如下：

（1）根據(jù)問題所屬學科中的規(guī)律、定律、公式，用微分方程與初始條件進行描述。

F（y,y',..yn-'\t）=0

y（0）=Jo,y（O）=^i,..yn-l＞（0）=y?-\

寫為向量的形式為y=[y；MD；M2）,…,?。╩一1）],〃與加可以不等。

（2）運用數(shù)學中的變量替換：為=y（n-l）,為/可（n?2）,...,yi=y\=y，把高階（大于2

階）的方程（組）寫成一階微分方程組：

%力”,y)

1力(“)

y=

7?

相應的初始條件為：

（3）根據(jù)（1）與（2）的結(jié)果，編寫能計算導數(shù)的M-函數(shù)文件odefileo

（4）將文件odefile與初始條件傳遞給求解器solver中的一個，運行后就可得到ODE

在指定時間區(qū)間上的解列向量y（其中包含j，及不同階的導數(shù)\

【例12-4]求解器solver應用實例。采用ode45求解描述某非剛性體的運動方程的

”=y2y3伊（o）=o

微分方程：\yi=一川3，其初始條件為了2（。）=1,并畫出解的圖形。

y'3=0.51^1^2[乃（。）=1

解：編寫函數(shù)文件rigid.m:

functiondy=rigid(tzy)

dy=zeros(3,1);%行向量

dy(1)=y(2)*y(3);

dy(2)=-y(l)*y(3);

dy(3)=-0.51*y(l)*y(2);

建立文件exl204.m,如下所示：

284????

第12章常微分方程求解

%exl204.m

11

options=odeset(*RelTol*zle-4,AbsTol,[le-4le-4le-5]);

[T,Y]=ode45(@rigid,[012],[011],options);

1

plot(TZY(:Z1)，」,T,Y(:,2)J-?1T,Y(:,3),J)

legend(?yl\'y2\'y3')

運行程序，輸出的圖形結(jié)果為圖12-1所示。

圖12-1非剛性體運動的微分方程圖

12.2歐拉法

12.2.1簡單歐拉法

歐拉法(Euler)是簡單有效的常微分方程數(shù)值解法，歐拉法有多種形式的算法，其中

簡單歐拉法是一種單步遞推算法。

1.基本原理

對于一個簡單的一階微分方程：

y=f(x,y)

其中初始值為：

y(x0)=y0

歐拉方法等同于將函數(shù)微分轉(zhuǎn)換為數(shù)值差分，如下式,

y(x)J.+—

〃h

故原方程變形為：

yn+\=y?+hf(x,l,y?)

i???285

精通MATLAB科學計算(第2版i---------------

2.算法的程序?qū)崿F(xiàn)

在MATLAB中編程實現(xiàn)的歐拉法函數(shù)為：MyEulero

功能：以歐拉法求解常微分方程。

調(diào)用格式:[outx.outy]=MyEuler(fun,xO,xt,jX),PointNum)0

其中，fun為函數(shù)名；

xO為自變量區(qū)間初值；

“為自變量區(qū)間終值；

可為函數(shù)在xO處的值；

PointNum為自變量在[M),xf]上所取的點數(shù)；

outx為輸出自變量區(qū)間上所取點的x值；

outy為對應點上的函數(shù)值。

歐拉法的MATLAB程序代碼如下：

function[outx,outy]=MyEuler(fun,xOzxt,yO,PointNum)

金用前向差分的歐拉方法解微分方程y^=fun

率函數(shù)f(x,y)：fun

務自變量的初值和終值：xOzxt

%y0表示函數(shù)在xO處的值，輸入初值為列向量形式

?自變量在[xO,xt]上取的點數(shù):PointNum

%outx:所取的點的x值

%outy:對應點上的函數(shù)值

ifnargin<5|PointNum<=03如果函數(shù)僅輸入4個參數(shù)值，則PointNum默認值為100

PointNum=100;

end

ifnargin<4%y0默認值為0

y0=0;

end

h=(xt-xO)/PointNum;,計算步長h

x=xO+[0:PointNum]'*h;,自變量數(shù)組

yd,：)=yO(:)z；，輸入存為列向量

fork=1:PointNum

f=feval(fun,x(k),y(k,:))8計算f(x,y)在每個迭代點的值

f=f(:)z;

y(k+1,:)=y(k,:)+h*;金對于所取的點x迭代計算y值

end

outy=y;

outx=x;

plot(x,y),畫出方程解的函數(shù)圖

286????

-------------------------------------------------------第章常微分方程求解

3.實例分析

【例12-5】歐拉法求解常微分方程應用實例。用歐拉法求常微分方程：

y=sinx+y

y(xo)=l,xo=0

當〃=0.2和〃=0.4時的數(shù)值解，與精確解進行對比，并畫出解的圖形。

解：首先建立函數(shù)文件myfunOl.m:

functionf=myfun01(xzy)

f=sin(x)+y;

通過在相同的積分區(qū)間上設定不同的取點數(shù)，從而改變步長，可得到不同的歐拉解。

以下程序即可實現(xiàn)此功能，并給出了幾乎等于精確解的微分方程符號解，如文件exl201.m

所示。

%exl205.m

functionexl205()

[xl,yl]=MyEuler('myfunOl',0,2*pi,1,16);率歐拉法所得的解

hl=2*pi/15莒計算取步長

[xllAyll]=MyEuler(*myfun01'，0,2*pi,l,32);》歐拉法所得的解

h2=pi/15,計算步長

y=dsolve(1Dy=y+sin(t),,1y(0)=11);%該常微分方程的符號解

fork=l:33

t(k)=xll(k);

y2(k)=subs(y,t(k));務求其對應點的茗散解

end

plot(xl,yl,1+b1,xll,yll,,og,,xll,y2,**rf)

legend('h=0.4的歐拉法解I，h=0.2的歐拉解I，符號解，)

運行程序后，輸出如圖12-2所示的圖形。

+h=0.4的歐拉法解

800h=0.2的歐拉解

?符號解

Y???287

精通MATLAB科學計算(第2忡--------------------------------

圖12-2歐拉法所得方程解與符號解(精確解)對比

從圖12-2中可以看出，用歐拉法得到的解和用符號法得到的解之間存在一定的誤差，

且取的步長越小，歐拉解越接近精確解。

由于向前差商比較簡單，此處僅舉向前差商的例子。另外還有兩種差商，分別為向后

差商和中心差商。有興趣的讀者可以參考相關資料，編寫相應的MATLAB程序，其基本

思想是一致的0

12.2.2改進的歐拉法

改進的歐拉法實際上是通過預估-校正過程對梯形公式進行修正從而簡化迭代過程的

方法，又稱為Henu算法。

1.基本原理

改進的歐拉法的具體求解過程如下：

對方程方=f(x,y)兩邊由X”到x.+i積分,并利用梯形公式,可得：

yn+\=y?

其中初始值為：

y(xo)=yo

由上可知，所得到的迭代式為隱式格式。計算中為了保證一定的精確度，避免迭代過

程不菲的計算量，一般先用顯式公式算出初始值，再用隱式公式進行一次修正。此過程稱

為預估-校正過程。

具體迭代公式如下：

K+I=y?+hf(xn,y?)

M>+i=yn+5/⑸,y”)+/(x?+i,yn+\)]

2.算法的程序?qū)崿F(xiàn)

在MATLAB中編程實現(xiàn)的改進的歐拉法函數(shù)為：MyEulerProo

功能：用改進的歐拉法求解常微分方程。

調(diào)用格式:[Xout,Yout]=MyEulerPro(fun,xO,xt,PointNumber)o

其中，ftin為函數(shù)名；

M)為自變量區(qū)間初值；

X，為自變量區(qū)間終值；

川為函數(shù)在xO處的值；

288????

第12章常微分方程求解

PointNumber為自變量在［xO,詞上所取的點數(shù)；

Xout為輸出自變量區(qū)間上所取點的x值；

Yout為對應點上的函數(shù)值。

改進的歐拉法的MATLAB程序代碼如下：

function［Xout,Yout］=MyEulerPro(fun,xO,xtzyO,PointNumber)

%用改進的歐拉法解微分方程y，=fun

%函數(shù)f(x,y)：fun

當自變量的初值和終值：xO,xt

%y0表示函數(shù)在xO處的值，輸入初值為列向量形式

許自變量在［xO,xt］上取的點數(shù):PointNumber

%Xout:所取的點的x值

%Yout:對應點上的函數(shù)值

ifnargin<5IPointNumer<=0告如果函數(shù)僅輸入4個參數(shù)值，則PointNumer默認值

為100

PointNumer=100;

end

ifnargin<4%y0默認值為0

y0=0;

end

h=(xt-xO)/PointNumber;為計算所取的兩離散點之間的距離

x=x0+[0:PointNumber]'*h;%表示出離散的自變量x

y(lr:)=y0(:)^;

fori=l:PointNumber常迭代計算過程

fl=h*feval(fun,x(i),y(i,:))；

fl=fl(:)

f2=h*feval(fun,x(i+1)zy(i,:)+fl);

f2=f2(:)

y(i+l,:)=y(i,:)+1/2*(fl+f2);

end

Xout=x;

Yout=y;

3.實例分析

【例12-6】改進的歐拉法求解常微分方程應用實例。利用改進的歐拉法求常微分方

程：

y=sinx+y

y(xo)=l,xo=0

即對【例12-4］用改進的歐拉法求解，比較兩種解法的結(jié)果，并畫出解的圖形。

解：編寫MATLAB程序exl206.m：

functionexl206()

」???289

精通MATLAB科學計算(第2蝌

和exl206比較改進的歐拉法,簡單歐拉方法以及微分方程符號解

[x3,y3]=MyEulerPro(*myfun01',0,2*pi,1,128);

[x,yl]=MyEuler(*myfun01*,0,2*pi,1,128);超歐拉法所得的解

y=dsolve('Dy=y+sin(t),,,y(0)=11);超該常微分方程的符號解

fork=l:129

t(k)=x(k);

y2(k)=subs(y,t(k));

招求其對應點的離散解

end

plot(x,ylJ-blx3,y3,Pg',x,y2J*r，)

legend(，簡單歐拉法解，J改進歐拉解，J符號解

，)

運行程序后，輸出如圖12-3所示的圖形。

如圖12-3所示，改進歐拉解與精確解幾乎完

全吻合，而簡單歐拉解與精確解還有一定的誤差。

圖12-3改進的歐拉法解、歐拉解、

精確解的對比圖由此可見，改進的歐拉法較之簡單歐拉法更為

精確。

290????

第12章常微分方程求解

12.3龍格甚塔法

盡管改進的歐拉法相對于簡單歐拉法較為精確。但是對于很多實際的問題，運用這兩

種方法仍然得不到足夠精確的解。龍格-庫塔法(Runge-Kutta)是較之前兩種方法計算精

度更高的方法。

1.基本原理

在龍格-庫塔法中，四階龍格-庫塔法的局部截斷誤差約為。(『)，被廣泛應用于解微分

方程的初值問題。其算法公式為：

h

%+1—yn+二(ki+2k2+2k3+左4)

6

其中：

ki=f(x?,yn)

ki=J\xn+^h,y?+；/?肩)

左3=f(x?+；；萬左2)

3)

k4=f(x?+h,y?+府

2.四階龍格-庫塔算法編程實現(xiàn)

在MATLAB中編程實現(xiàn)的四階龍格-庫塔算法函數(shù)為：MyRunge_Kuttao

功能：用四階龍格-庫塔算法求解常微分方程。

調(diào)用格式:[x,y]=MyRunge_Kutta(ftin,xO,xl,yO,PointNum,varargin)o

其中，fbn為函數(shù)名；

xO為自變量區(qū)間初值；

口為自變量區(qū)間終值；

yO為函數(shù)在M)處的值；

PointNum為自變量在[xO,xf]上所取的點數(shù)；

Varargin為可選項參數(shù)；

x為輸出自變量區(qū)間上所取點的x值；

y為對應點上的函數(shù)值。

四階龍格-庫塔法的MATLAB程序代碼如下：

function[x,y]=MyRunge_Kutta(fun,xOzxtzyO,PointNum,varargin)

%Runge-Kutta方法解微分方程形為y"(t)=f(x,y(x))

考此程序可解高階的微分方程。只要將其形式寫為上述微分方程的向量形式

????291

精通MATLAB科學計算(第2蚪

殳函數(shù)f(x,y)：fun

,自變量的初值和終值：xO,xt

SyO表示函數(shù)在xO處的值，輸入初值為列向量形式

*自變量在［xO,xt］上取的點數(shù):PointNum

%varargin為可選輸入項可傳適當參數(shù)給函數(shù)f(x,y)

%x:所取的點的x值

務y:對應點上的函數(shù)值

ifnargin<4|PointNum<=0

PointNum=100;

end

ifnargin<3

yO=0;

end

y(1/：)=yO(：)1；e初值存為行向量形式

h=(xt-xO)/(PointNum-1);書計算步長

x=x0+[0:PointNum]**h;%得x向量值

fork=1:PointNum當?shù)嬎?/p>

fl=h*feval(fun,x(k),y(k,:),varargin{:});

fl=fl(:)%得公式中kl

f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2Zvarargin{1});

f2=f2(:)1;超得公式中k2

f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2zvarargin{:});

f3=f3(:)彩得公式中k3

f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});

￡4=f4⑴，告得公式中k4

y(k+1,:)=y(k,:)+(fl+2*(f2+f3)+f4)/6;y(n+1)

end

3.實例分析

【例12-7］龍格-庫塔法求解常微分方程應用實例。采用龍格-庫塔法求解微分方程，

y(x())=O,xo=0

比較與簡單歐拉法、改進的歐拉法解的求解結(jié)果以及各種方法的運行時間，并畫出解

的圖形。

解：容易得到該微分方程的真值解為：

y=i-exp(-x)

MATLAB程序如下：

exl207.m

292???>

第12章常微分方程求解

常用三種不同方法解微分方程

clear,elf告清內(nèi)存中的變量

xO=O;

xt=2;

Num=100;

h=(xt-xO)/(Num-1);

x=xO+[0:Num]*h;

a=1;

yt=1-exp(-a*x);告真值解

fun=inline('-y+l*z*x*z*y');招用inline構(gòu)造函數(shù)f(x,y)

yO=0;號設定函數(shù)初值

PointNum=4;%設定取點數(shù)

[xl,ye]=MyEuler(fun,xO,xt,y0fPointNum);

[x2zyh]=MyEulerPro(fun,xO,xt,y0zPointNum);

[x3fyr]=MyRunge_Kutta(fun,xO,xtzyO,PointNum);

11

plot(xzytzk'zxlzye,'b:,x2,yh,*g:',x3zyr,'r:')

legend(1真值1,1簡單歐拉法解1,,改進的歐拉法解r,*Runge_Kutta法解1)

holdon

plot(xl,ye,*bo*,x2,yh,*b+',x3,yr,1r*1)

PointNum=1000;%估計各算法運行PointNum步迭代的時間

tic電計時開始

[xl,ye]=MyEuler(fun,xO,xt,yO,PointNum);

time_Euler=toe帶得到歐拉法的運行時間

tic

[xlzyh]=MyEulerPro(fun,xO,xtzyO,PointNum);

time_EulerPro=toe專得到改進的歐拉法運行時間

tic

[xlzyr]=MyRunge_Kutta(fun,x0zxtzyO,PointNum);

time_RK4=toe宅得至！IRunge_Kutta法運彳亍時間

運行后得到：

>>exl207

time_Euler=0.2544

0.2544

time_EulerPro=0.4887

0.4887

time_RK4=1.0195

1.0195

運行的結(jié)果如圖12-4所示。

Y???293

精通MATLAB科學計算(第2蝌

由運算結(jié)果可知，龍格-庫塔法可得到與真值幾乎相同的微分方程數(shù)值解。但其相應的

程序運行時間較長，幾乎是簡單歐拉法的五倍(歐拉法運行時間是0.2544s,而龍格-庫塔

法運行時間是1.0195s卜

4.求解器solvei■中的龍格-庫塔法求解

前面講到，求解器solver中的ode23采用二階、三階龍格-庫塔法；ode45采用四階、

五階龍格-庫塔法。

【例12-8】求解器solver中的龍格-庫塔法求解應用實例1。分別采用二、三、四、

五階龍格-庫塔方法求解以下方程，

y'=-3y2+2x2+3x,H0?x?1

<

j(xo)=l,xo=0

解法一：用二階、三階龍格-庫塔函數(shù)求解。

編寫程序文件ex1208a.m:

%exl208a.m用ode23得到微分方程解并計算出該算法運行時間

fun=inline('-3*yA2+2*x.A2+3*x','x','y');%用inline構(gòu)造函數(shù)f(x,y)

[x,y]=ode23(fun,(0,1],1);3可得到x,y輸出向量值

ode23(fun,[0,1],1),holdon&可得到輸出的函數(shù)圖

結(jié)果如圖12-5a所不。

294????

-------------------------------------------------------第章常微分方程求解

解法二：用四階、五階龍格-庫塔函數(shù)求解。

編寫程序文件ex1208b.m:

%exl208b.m用ode45得到微分方程解，并計算出該算法運行時間

fun=inline('-3*yA2+2*x.A2+3*x\'x1,*yT);告用inline構(gòu)造函數(shù)f(x,y)

ode45(fun,[0,1],1),holdon%可得到輸出的函數(shù)圖

tic

[x,y]=ode45(fun,[0,1],1);

tl=toc

tic

[x,y]=ode23(fun,[0,1],1);

t2=toc

命令窗口顯不tl、t2值。

?exl208b

tl=0.0256

0.0256

t2=0.0166

0.0166

如圖12-5b所示顯示函數(shù)的解結(jié)果。

????295

精通MATLAB科學計算(第2蝌

由圖12-5a和圖12-5b所示可知，用兩個函數(shù)求解得到的結(jié)果相同。而ode23比ode45

的運行速度要快一些。

【例12-9】求解器solver中的龍格-庫塔法求解應用實例2。采用ode45求解如下二

階方程：

V'(2)=X2)[-0.5+0.02XD]

初始條件為：x=0時，XI)=25,y(2)=2,并畫出解的圖形。

解：用ode45求解。首先建立函數(shù)文件myfun02.m:

myfun02.m

functiondx=myfun02(xzy)

dx=zeros(2,1);

dx(l)=y(l)*(l-0.1*y(2));

dx(2)=y(2)*(-0.5+0.02*y(1));

編寫exl209.m對微分方程求解

%exl209.m

[x,y]=ode45('myfun02'z[015],[252]);

plot(x,y(:,l),'-',x,y(:,2),'*'),畫出y(1),y(2)的函數(shù)圖

運行程序，輸出結(jié)果如圖12-6所示。

296????

常微分方程求解

12.4預估-校正法

所謂預估-校正(predictor-corrector)算法是指算法中主要包括的兩個過程，預估過程

和校正過程。前面講到的改進的歐拉法是一種簡單的預估-校正算法。此外還有兩種比較常

用的形式，它們分別是ABM(Adams-Bashfbrth-Moulton)法和Hamming法。

12.4.1ABM法

ABM法是最常用的傳統(tǒng)的預估-校正法。

1.基本原理

ABM法與上面幾種方法的主要不同是：它屬于多步算法，%的值是由，

,(X._2,y”-2),(x“_3,%-3)這幾個點通過預估-校正過程求得的，要分以下兩步

進行。

第一步:用四階的拉格朗日多項式來近似得到/(XJ)的積分，并用前四個點的/值表

示。從而可得ya+i的預估表達式p“+i：

pn+\=yn+&(-9加3+37加2-59/t+55,/?)

第二步：對p向的預估值進行校正，得到為+i的校正值金+|。

c“+i=%+一5￡T+19人+9/m)

通常使用的是改進的ABM算法。即在上述迭代式的基礎上，再加入更正公式，因此

可得如下的迭代運算公式：

????297

精通MATLAB科學計算(第2忡-----------------------------

P"+1=y.+&(-9加3+37小-59小+556)

24

251、

加〃+1—Pn+\+z~Pn)

?！?1=yn+~^^(/〃一2一5力】一1+19/〃+9/(X〃+],〃[”+]))

19,、

y/i+i=c〃+i—270^C，，+1-Pz)

由于算法采用多步近似，因而其具有較小的截斷誤差。①5),其計算精度優(yōu)于前面介

紹的眾多方法。

2.算法的程序?qū)崿F(xiàn)

MATLAB中的odel13函數(shù)就是該算法的程序?qū)崿F(xiàn)，由于迭代過程較復雜，此處不給

出具體程序。有興趣的讀者可以參考odel13的實現(xiàn)程序。其具體命令格式與上一節(jié)介紹

的ode23和ode45相同。

12.4.2Hamming法

Hamming法是另一種形式的多步預估一校正法，其截斷誤差也是。(二),運算效率與

ABM相當。下面對其進行具體介紹。

1.基本原理

Hamming預估-校正公式為：

4h“

Pn+\=y?-3+—(2.AI-2-fn-\+2力,)

112、

叫i+l=Pn+\+z—P")

cn+i=1[9y?-y?-2+3A(-/；,_i+2fn+f(xn+i,mn+\))]

8

9,、

yn+l=Q+l--P〃+l)

2.算法程序?qū)崿F(xiàn)

在MATLAB中編程實現(xiàn)的Hamming算法函數(shù)為：MyHamming。

功能：以Hamming算法求解常微分方程。

調(diào)用格式：[x,y]=MyHammingGx0,xt,j'O,N,KC)O

其中，/為函數(shù)名；

xO為自變量區(qū)間初值；

W為自變量區(qū)間終值；

298????

第12章常微分方程求解

yO為函數(shù)在xO處的值；

N為自變量在［xO，M上所取的點數(shù)；

KC為確定是否使用改正公式；

x為輸出自變量區(qū)間上所取點的x值；

y為對應點上的函數(shù)值。

Hamming算法的MATLAB程序代碼如下：

function［x,y］=MyHamming(f,xO,xt,yO,NzKC)

%用hamming方法解向量微分方程yz(t)=f(t,y(t))

軟函數(shù)f(x,y)：f

告自變量的初值和終值：xO,xt

%函數(shù)在xO處的值：yO

出自變量在［xO,xt］上取的點數(shù)：N

為選擇是否使用改正公式：KC

生x:所取的點的x值

%y:對應點上的函數(shù)值

ifnargin<6,

KC=1;帶默認使用更正公式

end

ifnargin<5IN<=0

N=100;與最大迭代數(shù)默認為100

end

ifnargin<4

yO=0;考初值默認為0

end

yO=yO(:)1;%使yO為行向量

h=(xt-xO)/(N-l);超步長

xtO=x0+2*h;

[x,y]=MyRunge_Kutta(fzxO,xtO,yOz3)務用Runge-Kutta得初始三點值

x=[x(1:3)*x(4):h:xt]1;務重新調(diào)整x向量，使其前三點與

Runge-Kutta法計算的三點相同

fork=2:4

F(k-1,:)=feval(f,x(k)zy(k,:));告計算f2f3f4

end

P=y(4,:);%預估量初值

c=y(4,:);超校正量初值

h34=h/3*4;稱預估公式中系數(shù)

KC11=KC*112/121;務更正公式中系數(shù)

KC91=KC*9/121;與最后計算y值的公式中更正項的系數(shù)

h312=3*h*[-121];當校正公式中各f項系數(shù)

詞<<<299

精通MATLAB科學計算(第2版i-----------------------------------------

fork=4:N-1

pl=y(k-3,:)+h34*(2*(F(l,:)+F(3Z:))-F(2,:));%p(n+l)計算公式

ml=pl+KC11*(c-p);當m(n+l)計算公式

cl=(-y(k-2,:)+9*y(kz:)+...

h312*[F(2:3,f