常微分方程的解法

1、-175-第十五章常微分方程的解法建立微分方程只是解決問題的第一步，通常需要求出方程的解來說明實際現(xiàn)象，并加以檢驗。如果能得到解析形式的解固然是便于分析和應用的，但是我們知道，只有線dx分重要的手段。性常系數(shù)微分方程，并且自由項是某些特殊類型的函數(shù)時，才可以肯定得到這樣的解， 而絕大多數(shù)變系數(shù)方程、非線性方程都是所謂“解不出來”的，即使看起來非常簡單的 方程如dy . y2 x2，于是對于用微分方程解決實際問題來說，數(shù)值解法就是一個十§1常微分方程的離散化F面主要討論一階常微分方程的初值問題，其一般形式是y(a)在下面的討論中，我們總假定函數(shù)f (x, y)連續(xù)，且關(guān)于y滿足李普希茲

2、(Lipschitz)條件，即存在常數(shù)L ,使得I f(x,y) f(x,y)匚 L |y y |這樣，由常微分方程理論知，初值問題(1)的解必定存在唯一。所謂數(shù)值解法，就是求問題(1)的解y(x)在若干點a = x0 ：捲：x2: xN = b處的近似值yn( n=1,2,，N)的方法，y.( n=1,2,，N)稱為問題(1)的數(shù)值解， hn二Xn 1 -Xn稱為由X.到人1的步長。今后如無特別說明，我們總?cè)〔介L為常量h。建立數(shù)值解法，首先要將微分方程離散化，一般采用以下幾種方法：(i )用差商近似導數(shù)若用向前差商y(xn 1 蟲"代替y'(xn)代入(1)中的微分方程，則

3、得y(Xn AWXJf(xn,y(xn)(n =0,1, )h化簡得y(Xn 1)： y(Xn) hf (Xn, y(Xn)如果用y(Xn)的近似值yn代入上式右端，所得結(jié)果作為y(Xn1)的近似值，記為yn .1，則有yn 十二 yn +hf(Xn,yn) (n =0,1,)這樣，問題(1)的近似解可通過求解下述問題丫時=yn +hf(Xnn)(門=0,1，) 0= y(a)得到，按式(2)由初值y0可逐次算出y1, y2/。式(3)是個離散化的問題，稱為差 分方程初值問題。需要說明的是，用不同的差商近似導數(shù)，將得到不同的計算公式。(ii )用數(shù)值積分方法將問題(1)的解表成積分形式，用數(shù)值

4、積分方法離散化。例如，對微分方程兩端 積分，得xny(Xn 1)-y(Xn)二.f(x, y(x)dx (n =0,1;)xn右邊的積分用矩形公式或梯形公式計算。(iii) Taylor多項式近似將函數(shù)y(x)在Xn處展開，取一次 Taylor多項式近似，則得y(Xn i) ： y(Xn) hy'(Xn)二 y(Xn) hf(Xn,y(Xn)再將y(Xn)的近似值yn代入上式右端，所得結(jié)果作為y(Xn.J的近似值yn 1，得到離散化的計算公式y(tǒng)n 1 =yn hf (Xn,yn)以上三種方法都是將微分方程離散化的常用方法，每一類方法又可導出不同形式的計算公式。其中的Taylor展開法，

5、不僅可以得到求數(shù)值解的公式，而且容易估計截斷誤£。§歐拉(Euler)方法2.1 Euler 方法Euler方法就是用差分方程初值問題(5)'yn+ = yn +hf(Xn,yn) (n = 0,1，)<y。二 y(a)的解來近似微分方程初值問題(1 )的解，即由公式(2)依次算出y(xn)的近似值yn(n =1,2，)。這組公式求問題(1)的數(shù)值解稱為向前Euler公式。如果在微分方程離散化時，用向后差商代替導數(shù)，即y'(xny(Xn d)，h則得計算公式>n+ = Yn +hf (Xn*ynG 5 = 0,1，) (6)yo= y(a)用這組

6、公式求問題(1)的數(shù)值解稱為向后 Euler公式。向后Euler法與Euler法形式上相似，但實際計算時卻復雜得多。向前Euler公式是顯式的，可直接求解。向后Euler公式的右端含有yn 1，因此是隱式公式，一般要用迭代法求解，迭代公式通常為yn =yn+hf(Xn,yn)y yn hf(Xn1,ynk)1)(k =0,1,2,)2.2 Euler方法的誤差估計對于向前Euler公式(5)我們看到，當n =1,2/時公式右端的yn都是近似的, 所以用它計算的 yn d會有累積誤差，分析累積誤差比較復雜，這里先討論比較簡單的 所謂局部截斷誤差。假定用(5)式時右端的yn沒有誤差，即yn =：y

7、(xn)，那么由此算出yn 1 =y(Xn) hf (Xn, 丫區(qū))局部截斷誤差指的是，按(7)式計算由Xn到Xn .1這一步的計算值yn，1與精確值之差y(Xn 1) - yn 1。為了估計它，由Taylor展開得到的精確值 y(Xn 1)是h23y(xn 1)=y(xn) hy'(xn)y''(Xn) O(h )2(7 )、( 8 )兩式相減(注意到 y'=f(x,y)得h232y''(Xn) O(h ) ： O(h )2y(Xn 1)- yn 1(7)y(Xn 1)(8)(9)即局部截斷誤差是h2階的，而數(shù)值算法的精度定義為：若一種算法的局

8、部截斷誤差為o(hp 1)，則稱該算法具有p階精度。顯然p越大，方法的精度越高。式(5)說明，向前 Euler方法是一階方法，因此 它的精度不高。§3 改進的Euler方法3.1 梯形公式利用數(shù)值積分方法將微分方程離散化時，若用梯形公式計算式(4)中之右端積分，即Xn 1hx f(x, y(x)dx f (Xn, y(Xn) f (Xn 1, y(Xn 1)Xn并用yn, yn 1代替y(Xn), y(Xn 1)，則得計算公式hyn 1 二 yn f (Xn, yn) f (Xn 1, Yn 1)2這就是求解初值問題(1)的梯形公式。直觀上容易看出，用梯形公式計算數(shù)值積分要比矩形公式

9、好。梯形公式為二階方法。梯形公式也是隱式格式，一般需用迭代法求解，迭代公式為yn +hf (Xn, yn)WnV0 = yn +£ f(Xn, yn) + f (Xn 卅,yn'i)(10)2(k =0,1,2,)由于函數(shù)f (x, y)關(guān)于y滿足Lipschitz條件，容易看出I (k 1)(k) 1 . hL (k)(k -4).| yn 1- yn 1 | yn 1 一 yn 1 |2其中L為Lipschitz常數(shù)。因此，當0 ：比:1時，迭代收斂。但這樣做計算量較大。2如果實際計算時精度要求不太高，用公式(10)求解時，每步可以只迭代一次，由此導出一種新的方法一改進

10、Euler法。3.2 改進Euler法按式(6)計算問題(1)的數(shù)值解時，如果每步只迭代一次，相當于將Euler公式與梯形公式結(jié)合使用：先用Euler公式求yn 1的一個初步近似值 yn ,，稱為預測值，然后用梯形公式校正求得近似值yn ,，即yn yn hf(Xn,yn)預測th_ 亠(11)j yn 1 = yn ' 2 f ( Xn , yn ) f (xn 1, yn 1 )校正式(11)稱為由Euler公式和梯形公式得到的預測一校正系統(tǒng)，也叫改進Euler法。為便于編制程序上機，式(11)常改寫成yp = Yn +hf(Xn,yn)« yq = yn +hf (Xn

11、 +h, yp)(12)1 %卅石仏+yq)改進Euler法是二階方法。§4 龍格一庫塔(RungeKutta)方法回到Euler方法的基本思想一用差商代替導數(shù)一上來。實際上，按照微分中值定理 應有心卅)-畑)*仏+餉),0卡1h注意到方程y'=f(x, y)就有y(Xnd =y(Xn) +hf (Xn +Th,y(Xn + 日h)( 13)不妨記K = f(xn rh,y(xn rh),稱為區(qū)間人公“計上的平均斜率?？梢娊o出一種 斜率K , (13)式就對應地導出一種算法。向前Euler公式簡單地取f(Xn,yn)為K，精度自然很低。改進的 Euler公式可理 解為 K 取

12、 f(Xn，yn) , f(Xn 1，Vn 1)的平均值，其中 1 一 Yn ' hf(Xn,yn)，這種處 理提高了精度。如上分析啟示我們，在區(qū)間Xn,Xn內(nèi)多取幾個點，將它們的斜率加權(quán)平均作為 K，就有可能構(gòu)造出精度更高的計算公式。這就是龍格一庫塔方法的基本思想。4.1首先不妨在區(qū)間洙人1內(nèi)仍取2個點，仿照(13)式用以下形式試一下yn 1 =yn h(、k1& 二 f(Xn, yn)(14)阮= f(Xn +呦,0£°少£1其中-2/ /:為待定系數(shù)，看看如何確定它們使(14)式的精度盡量高。為此我們分析局部截斷誤差y(Xn 1) - yn

13、1，因為yn =y(Xn)，所以(14 )可以化為yn = y(Xn) +h(人+bk2)(15)& = f (Xn, y(Xn) = y'(Xn) 也=f(Xn +ah, 丫佃)+ 陽3=f (Xn, y(Xn) +ahfx(Xn, y(Xn)L+ Phki fy(Xn, y(Xn) +O(h2)其中k2在點(Xn,y(Xn)作了 Taylor展開。(15)式又可表為203yn 1 二 y(Xn) ( 12)hy'(Xn)2： h (fx ff y ) O(h )CL注意到y(tǒng)(Xn 1)= y(Xn) hy'(Xn)h223y''(Xn) O(

14、h )y'' = fx - ffy，可見為使誤差 y(Xn .J - yn 1=O(h3)，只須令(16)16）式有4個未知數(shù)待定系數(shù)滿足（16）的（15）式稱為2階龍格一庫塔公式。由于（1而只有3個方程，所以解不唯一。不難發(fā)現(xiàn)，若令、=，2 =丄，=：=1，即為改2進的Euler公式?？梢宰C明，在Xn，Xn 1內(nèi)只取2點的龍格一庫塔公式精度最高為2階。4.24階龍格一庫塔公式要進一步提高精度，必須取更多的點，如取4點構(gòu)造如下形式的公式：An+ = Yn +h（箱& 十婦k2 +%k3 +k4）k1 = f （Xn, yn）“ k2 = f （xn 中h, yn +（h

15、k1）（ 17）k3 = f （Xn +口2h, yn + hk jhk2）k4 二 f （Xn ： 3h, yn ：4hk1：5hk2：6hk3）其中待定系數(shù)'i/'i, '-i共13個，經(jīng)過與推導2階龍格一庫塔公式類似、但更復雜的計算，得到使局部誤差 y（xn J -yn 4 =o（h5）的11個方程。取既滿足這些方程、 又較簡 單的一組id,，可得(18)hyn = (ki +2k2 +2k3 +k4)6ki = f (Xn,yn)hhk1杯2 = f (Xn +二，y)2 2, 一 h hk2k3 = f(Xn +：, yn +片)k f (Xn h, yn h

16、k3)這就是常用的4階龍格一庫塔方法(簡稱 RK方法)。§5線性多步法以上所介紹的各種數(shù)值解法都是單步法，這是因為它們在計算yn “時，都只用到前一步的值,單步法的一般形式是yn* = yn +h®(Xn,yn,h) (n = 0,1,，N -1)( 19)其中(x, y,h)稱為增量函數(shù)，例如 Euler方法的增量函數(shù)為 f (x, y)，改進Euler法的 增量函數(shù)為1(x,y,h)f(x,y) f(x h,y hf(x,y)2如何通過較多地利用前面的已知信息，如yn, yn,，yn _r，來構(gòu)造高精度的算法計算yn 1，這就是多步法的基本思想。經(jīng)常使用的是線性多步法。

17、讓我們試驗一下r二1，即利用Yn, Yn j計算Yn 1的情況。 從用數(shù)值積分方法離散化方程的(4)式Xny(Xn 1)-y(Xn)二.f (X, y(x)dxXn出發(fā)，記f (Xn, yn)二fn， f &n/,n)二仁4，式中被積函數(shù)f ( X, y(x)用二節(jié)點 (Xn j, fn 4)， (Xn, fn)的插值公式得到(因 X Xj，所以是外插。f(x,y(x)十亠 fn4_(20)Xn_ XnXn 二 一Xn1(x Xn4)fn (XXn)fn4 h此式在區(qū)間Xn,人J上積分可得冷 13hhX f(X, y(x)dX fnfn4Xn22于是得到.hyn 1 二 yn(3fn

18、- fn J( 21)2注意到插值公式(20)的誤差項含因子(X - Xn 4)(X - Xn )，在區(qū)間xXn上積分后-#-181-得出h3，故公式（21）的局部截斷誤差為 0（h3），精度比向前Euler公式提高1階。若取r =2,3，可以用類似的方法推導公式，如對于r = 3有h24剛一5嘰37f嘰）(22)其局部截斷誤差為 0（h5）。如果將上面代替被積函數(shù)f （x, y（x）用的插值公式由外插改為內(nèi)插，可進一步減小誤差。內(nèi)插法用的是 yn 1, ，yn_r 1，取r =1時得到的是梯形公式，取 r = 3時 可得hyn 1 = yn （9 fn I 19 fn -5fnfn ）（ 2

19、3）24與（22）式相比，雖然其局部截斷誤差仍為o（h5），但因它的各項系數(shù)（絕對值）大為減小，誤差還是小了。當然，（22）式右端的fn彳未知，需要如同向后 Euler公式一樣，用迭代或校正的辦法處理。§6 一階微分方程組與高階微分方程的數(shù)值解法6.1一階微分方程組的數(shù)值解法設有一階微分方程組的初值問題y'i = fi(x, y1,y2,ym) y(a) = yi°(i =1,2,m)若記 y =（%,丫2, , Ym）T，y。=（%0,丫20, ,ym0）T，問題（24）可寫成如下向量形式= (f1, f2,（24）fm）T，則初值如果向量函數(shù)；yf(x, y)、

20、y(a) = yof (x, y)在區(qū)域D :(25)連續(xù)，且關(guān)于y滿足Lipschitz條件，即存在L 0，使得對-a,b， y1, y Rm， 都有| f（x, yj - f （x,y2）咗 L 屮-光那么問題（25）在a,b上存在唯一解y = y（x）。問題（25）與（1）形式上完全相同，故對初值問題（1）所建立的各種數(shù)值解法可 全部用于求解問題（25）。6.2高階微分方程的數(shù)值解法高階微分方程的初值問題可以通過變量代換化為一階微分方程組初值問題。 設有m階常微分方程初值問題y(m)= f (x, y, y',y(m') a 蘭 x 蘭 b畀(a) = y。, y'

21、;(a) = y0°,，y"z(a) = y0m)引入新變量yy, yy' ,yy（m，），問題（26）就化為一階微分方程初值問題y'i = y2yi(a) = y°y； = y3y2(a y01)=(27),/ 、(m_2)y m = ymym_j(a) = y°(m 1) y m = f(X, yi,ym)ym(a)二yo然后用6.1中的數(shù)值方法求解問題(27),就可以得到問題(26)的數(shù)值解。最后需要指出的是，在化學工程及自動控制等領域中，所涉及的常微分方程組初值問題常常是所謂的“剛性”問題。具體地說，對一階線性微分方程組dy Ay

22、 T(x)( 28)dx其中y： Rm,A為m階方陣。若矩陣A的特征值i(i =1,2,-,m)滿足關(guān)系Re i ：0 (i =1,2,m)max | Re、min | Re i |1豈竺1豈劉則稱方程組(28)為剛性方程組或 Stiff方程組，稱數(shù)s = max | Re i |/min | Re i |1蘭哲| | 1生四| i 1為剛性比。對剛性方程組用前面所介紹的方法求解，都會遇到本質(zhì)上的困難，這是由數(shù)值方法本身的穩(wěn)定性限制所決定的。理論上的分析表明，求解剛性問題所選用的數(shù)值方法最好是對步長h不作任何限制。§ Matlab 解法7.1 Matlab數(shù)值解7.1.1非剛性常微分

23、方程的解法Matlab的工具箱提供了幾個解非剛性常微分方程的功能函數(shù)，如ode45, ode23,ode113，其中ode45采用四五階RK方法，是解非剛性常微分方程的首選方法，ode23采用二三階RK方法，ode113采用的是多步法，效率一般比 ode45高。Matlab的工具箱中沒有 Euler方法的功能函數(shù)。(I)對簡單的一階方程的初值問題；y = f(x,y)(x。)= y。改進的Euler公式為yp = yn +hf (Xn ,yn)Tq =% +hf (Xn +hp)1Yn+ =-(yYq)L2我們自己編寫改進的 Euler方法函數(shù)eulerpro.m如下：functionx,y=

24、eulerpro (fun, x0,xfi nal,y 0,n);ifnargin<5,n=50;e ndh=(xfi nal-xO)/n; x(1)=xO;y(1)=yO; for i=1: nx(i+1)=x(i)+h;y1= y(i)+h*feval(fu n, x(i),y(i); y2=y(i)+h*feval(fu n,x(i+1),y1); y(i+1)=(y1+y2)/2;end例1用改進的Euler方法求解2y' - -2y 2x 2x , (0 _ x _ 0.5), y(0) = 1 解編寫函數(shù)文件doty.m如下：fun ctio n f=doty(x,y

25、); f=-2*y+2*xA2+2*x;在Matlab命令窗口輸入：x,y=eulerpro('doty',0,0.5,1,10)即可求得數(shù)值解。(II)ode23，ode45, ode113的使用Matlab的函數(shù)形式是t,y=solver('F',tspan， y0)這里solver為ode45，ode23, ode113，輸入?yún)?shù)F是用M文件定義的微分方程 y' = f(x,y)右端的函數(shù)。tspan=t0,tfinal是求解區(qū)間，y0是初值。例2用R防法求解2y'=2y 2x 2x，(0 乞 x0.5)，y(0)=1解 同樣地編寫函數(shù)文件

26、 doty.m如下：fun ctio n f=doty(x,y); f=-2*y+2*xA2+2*x;在Matlab命令窗口輸入：x,y=ode45('doty',0,0.5,1)即可求得數(shù)值解。7.1.2剛性常微分方程的解法Matlab的工具箱提供了幾個解剛性常微分方程的功能函數(shù)，如ode15s，ode23s,ode23t，ode23tb，這些函數(shù)的使用同上述非剛性微分方程的功能函數(shù)。7.1.3高階微分方程y(n f (t, y, y',y(n)解法例3考慮初值問題y'''_3y''_y'y =0y(0)=0y'

27、(0)=1y''(0) = _1解(i)如果設y1y,y2y',y y''，那么Y1 (0) = 0y 二 y3y'3 =3y3y2y1V2(0) =1y3(0)一 1初值問題可以寫成 Y'=F(t,Y),Y(0) =丫。的形式，其中丫 二浙;丫2畀3。(ii)把一階方程組寫成接受兩個參數(shù)t和y，返回一個列向量的M文件F.m：fun ctio ndy=F(t,y);dy=y(2);y(3);3*y(3)+y(2)*y(1);注意：盡管不一定用到參數(shù)t和y , M文件必須接受此兩參數(shù)。這里向量dy必須是列向量。(iii )用Matlab解決

28、此問題的函數(shù)形式為T,Y=solver('F',tspa n, y0)這里solver為ode45、ode23、ode113，輸入?yún)?shù)F是用M文件定義的常微分方程組， tspan=tO tfinal是求解區(qū)間，yO是初值列向量。在Matlab命令窗口輸入T,Y=ode45('F',0 1,0;1;-1)就得到上述常微分方程的數(shù)值解。這里Y和時刻T是一一對應的，Y(:,1)是初值問題的解，Y(:,2)是解的導數(shù)，Y(:,3)是解的二階導數(shù)。例4 求van der Pol 方程y''-叩-y2)y'+y=o的數(shù)值解，這里 0是一參數(shù)。解(i

29、)化成常微分方程組。設y1 = y, y2 = y'，則有M = y2 丿 2 譏二W -力M沖(ii)書寫 M 文件(對于- 1)vdp1.m:function dy=vdp1(t,y); dy=y(2);(1-y(1)A2)*y(2)-y(1);(iii )調(diào)用Matlab函數(shù)。對于初值y(0) = 2, y'(0) =0，解為T,Y=ode45('vdp1' ,0 20,2;0);(iv )觀察結(jié)果。利用圖形輸出解的結(jié)果：plot(T,Y(:,1),'-',T,Y(:,2),'-')title('Soluti on

30、of van der Pol Equati on,m u=1');xlabel( 'time t' );ylabel('solutio n y');legend('y1', 'y2');例 5 van der Pol 方程，"=1000 （剛性）解(i) 書寫M文件vdplOOO.m:fun ction dy=vdp1OOO(t,y);dy=y(2);1000*(1-yA2)*y(2)-y(1);(ii )觀察結(jié)果t,y=ode15s( 'vdplOOO' ,0 3000,2;0);plot(t,

31、y(:,1),'o')title( 'Solutio n of van der Pol Equatio n,mu=1000');xlabel( 'time t' );ylabel( 'solution y(:,1)');7.2常微分方程的解析解在Matlab中，符號運算工具箱提供了功能強大的求解常微分方程的符號運算命令 dsolve。常微分方程在 Matlab中按如下規(guī)定重新表達：符號D表示對變量的求導。Dy表示對變量y求一階導數(shù)，當需要求變量的n階導數(shù)時，用Dn表示，D4y表示對變量y求4階導數(shù)。由此，常微分方程y'&#

32、39;，2y' = y在Matlab中，將寫成'D2y+2*Dy=y'。7.2.1求解常微分方程的通解無初邊值條件的常微分方程的解就是該方程的通解。其使用格式為：dsolve('diff_equati on')dsolve(' diff_equation'， 'var')式中diff_equation為待解的常微分，第1種格式將以變量t為自變量進行求解，第2種格式則需定義自變量var。例6試解常微分方程2X2y (x -2y)y'=0解編寫程序如下：syms x ydiff_equ='xA2+y+(x-2*

33、y)*Dy=0'dsolve(diff_equ,'x')7.2.2 -求解常微分方程的初邊值問題求解帶有初邊值條件的常微分方程的使用格式為：dsolve('diff_equation' , 'condition1 , condition2 ,''var') 其中condition1 , condition2 ,即為微分方程的初邊值條件。例7試求微分方程, y(1)=8, y'(1)=7,y''(2) = 4的解。解編寫程序如下：y=dsolve( 'D3y-D2y=x', '

34、y(1)=8,Dy(1)=7,D2y(2)=4', 'x')7.2.3求解常微分方程組求解常微分方程組的命令格式為：dsolve('diff_equ1 , diff_equ2 ,', 'condition1 , condition2 ,')dsolve('diff_equ1 , diff_equ2 ,', 'condition1 , condition2 ,', 'var')第1種格式用于求解方程組的通解軍，第2種格式可以加上初邊值條件，用于具體求解。例8試求常微分方程組：,'x')"f''+3g =s in