第十五章 常微分方程解法

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

2、理論知，初值問(wèn)題(1)的解必定存在唯一。所謂數(shù)值解法，就是求問(wèn)題（1）的解在若干點(diǎn) 處的近似值的方法，稱(chēng)為問(wèn)題（1）的數(shù)值解，稱(chēng)為由到的步長(zhǎng)。今后如無(wú)特別說(shuō)明，我們總?cè)〔介L(zhǎng)為常量。建立數(shù)值解法，首先要將微分方程離散化，一般采用以下幾種方法：（i）用差商近似導(dǎo)數(shù)若用向前差商代替代入（1）中的微分方程，則得 化簡(jiǎn)得如果用的近似值代入上式右端，所得結(jié)果作為的近似值，記為，則有 （2）這樣，問(wèn)題（1）的近似解可通過(guò)求解下述問(wèn)題 （3）得到，按式（2）由初值可逐次算出。式（3）是個(gè)離散化的問(wèn)題，稱(chēng)為差分方程初值問(wèn)題。需要說(shuō)明的是，用不同的差商近似導(dǎo)數(shù)，將得到不同的計(jì)算公式。（ii）用數(shù)值積分方法將問(wèn)題（

3、1）的解表成積分形式，用數(shù)值積分方法離散化。例如，對(duì)微分方程兩端積分，得 (4)右邊的積分用矩形公式或梯形公式計(jì)算。（iii）Taylor多項(xiàng)式近似將函數(shù)在處展開(kāi)，取一次Taylor多項(xiàng)式近似，則得再將的近似值代入上式右端，所得結(jié)果作為的近似值，得到離散化的計(jì)算公式 以上三種方法都是將微分方程離散化的常用方法，每一類(lèi)方法又可導(dǎo)出不同形式的計(jì)算公式。其中的Taylor展開(kāi)法，不僅可以得到求數(shù)值解的公式，而且容易估計(jì)截?cái)嗾`差。2 歐拉（Euler）方法2.1 Euler方法Euler 方法就是用差分方程初值問(wèn)題 （5）的解來(lái)近似微分方程初值問(wèn)題（1）的解，即由公式（2）依次算出的近似值。這組公式求

4、問(wèn)題（1）的數(shù)值解稱(chēng)為向前Euler公式。如果在微分方程離散化時(shí)，用向后差商代替導(dǎo)數(shù)，即，則得計(jì)算公式 （6）用這組公式求問(wèn)題（1）的數(shù)值解稱(chēng)為向后Euler公式。向后Euler法與Euler法形式上相似，但實(shí)際計(jì)算時(shí)卻復(fù)雜得多。向前Euler公式是顯式的，可直接求解。向后Euler公式的右端含有，因此是隱式公式，一般要用迭代法求解，迭代公式通常為 2.2 Euler方法的誤差估計(jì) 對(duì)于向前Euler公式（5）我們看到，當(dāng)時(shí)公式右端的都是近似的，所以用它計(jì)算的會(huì)有累積誤差，分析累積誤差比較復(fù)雜，這里先討論比較簡(jiǎn)單的所謂局部截?cái)嗾`差。假定用（5）式時(shí)右端的沒(méi)有誤差，即，那么由此算出 （7）局部截

5、斷誤差指的是，按（7）式計(jì)算由到這一步的計(jì)算值與精確值之差。為了估計(jì)它，由Taylor展開(kāi)得到的精確值是 （8）（7）、（8）兩式相減（注意到）得 （9）即局部截?cái)嗾`差是階的，而數(shù)值算法的精度定義為：若一種算法的局部截?cái)嗾`差為，則稱(chēng)該算法具有階精度。顯然越大，方法的精度越高。式（5）說(shuō)明，向前Euler方法是一階方法，因此它的精度不高。3 改進(jìn)的Euler方法3.1 梯形公式利用數(shù)值積分方法將微分方程離散化時(shí)，若用梯形公式計(jì)算式(4)中之右端積分，即并用代替，則得計(jì)算公式 這就是求解初值問(wèn)題（1）的梯形公式。直觀上容易看出，用梯形公式計(jì)算數(shù)值積分要比矩形公式好。梯形公式為二階方法。梯形公式也是

6、隱式格式，一般需用迭代法求解，迭代公式為 (10) 由于函數(shù)關(guān)于滿(mǎn)足Lipschitz條件，容易看出其中為L(zhǎng)ipschitz常數(shù)。因此，當(dāng)時(shí)，迭代收斂。但這樣做計(jì)算量較大。如果實(shí)際計(jì)算時(shí)精度要求不太高，用公式（10）求解時(shí)，每步可以只迭代一次，由此導(dǎo)出一種新的方法改進(jìn)Euler法。3.2 改進(jìn)Euler法按式（6）計(jì)算問(wèn)題（1）的數(shù)值解時(shí)，如果每步只迭代一次，相當(dāng)于將Euler公式與梯形公式結(jié)合使用：先用Euler公式求的一個(gè)初步近似值，稱(chēng)為預(yù)測(cè)值，然后用梯形公式校正求得近似值，即 （11）式（11）稱(chēng)為由Euler公式和梯形公式得到的預(yù)測(cè)校正系統(tǒng)，也叫改進(jìn)Euler法。為便于編制程序上機(jī)，式

7、（11）常改寫(xiě)成 （12）改進(jìn)Euler法是二階方法。4 龍格庫(kù)塔（RungeKutta）方法回到Euler方法的基本思想用差商代替導(dǎo)數(shù)上來(lái)。實(shí)際上，按照微分中值定理應(yīng)有 注意到方程就有 （13）不妨記，稱(chēng)為區(qū)間上的平均斜率。可見(jiàn)給出一種斜率，（13）式就對(duì)應(yīng)地導(dǎo)出一種算法。向前Euler公式簡(jiǎn)單地取為，精度自然很低。改進(jìn)的Euler公式可理解為取，的平均值，其中，這種處理提高了精度。如上分析啟示我們，在區(qū)間內(nèi)多取幾個(gè)點(diǎn)，將它們的斜率加權(quán)平均作為，就有可能構(gòu)造出精度更高的計(jì)算公式。這就是龍格庫(kù)塔方法的基本思想。4.1 首先不妨在區(qū)間內(nèi)仍取2個(gè)點(diǎn)，仿照（13）式用以下形式試一下 （14）其中為待

8、定系數(shù)，看看如何確定它們使（14）式的精度盡量高。為此我們分析局部截?cái)嗾`差，因?yàn)?，所以?4）可以化為 （15）其中在點(diǎn)作了Taylor展開(kāi)。（15）式又可表為注意到中，可見(jiàn)為使誤差，只須令 ， （16）待定系數(shù)滿(mǎn)足（16）的（15）式稱(chēng)為2階龍格庫(kù)塔公式。由于（16）式有4個(gè)未知數(shù)而只有3個(gè)方程，所以解不唯一。不難發(fā)現(xiàn)，若令 ，即為改進(jìn)的Euler公式?？梢宰C明，在內(nèi)只取2點(diǎn)的龍格庫(kù)塔公式精度最高為2階。4.2 4階龍格庫(kù)塔公式要進(jìn)一步提高精度，必須取更多的點(diǎn)，如取4點(diǎn)構(gòu)造如下形式的公式： （17）其中待定系數(shù)共13個(gè)，經(jīng)過(guò)與推導(dǎo)2階龍格庫(kù)塔公式類(lèi)似、但更復(fù)雜的計(jì)算，得到使局部誤差的11個(gè)方

9、程。取既滿(mǎn)足這些方程、又較簡(jiǎn)單的一組，可得 （18）這就是常用的4階龍格庫(kù)塔方法（簡(jiǎn)稱(chēng)RK方法）。5 線性多步法以上所介紹的各種數(shù)值解法都是單步法，這是因?yàn)樗鼈冊(cè)谟?jì)算時(shí)，都只用到前一步的值，單步法的一般形式是 （19）其中稱(chēng)為增量函數(shù)，例如Euler方法的增量函數(shù)為，改進(jìn)Euler法的增量函數(shù)為如何通過(guò)較多地利用前面的已知信息，如，來(lái)構(gòu)造高精度的算法計(jì)算，這就是多步法的基本思想。經(jīng)常使用的是線性多步法。讓我們?cè)囼?yàn)一下，即利用計(jì)算的情況。從用數(shù)值積分方法離散化方程的（4）式 出發(fā)，記，式中被積函數(shù)用二節(jié)點(diǎn)，的插值公式得到（因，所以是外插。 （20）此式在區(qū)間上積分可得 于是得到 （21）注意到插

10、值公式（20）的誤差項(xiàng)含因子，在區(qū)間上積分后得出，故公式（21）的局部截?cái)嗾`差為，精度比向前Euler公式提高1階。若取可以用類(lèi)似的方法推導(dǎo)公式，如對(duì)于有 （22）其局部截?cái)嗾`差為。如果將上面代替被積函數(shù)用的插值公式由外插改為內(nèi)插，可進(jìn)一步減小誤差。內(nèi)插法用的是，取時(shí)得到的是梯形公式，取時(shí)可得 （23）與（22）式相比，雖然其局部截?cái)嗾`差仍為，但因它的各項(xiàng)系數(shù)（絕對(duì)值）大為減小，誤差還是小了。當(dāng)然，（22）式右端的未知，需要如同向后Euler公式一樣，用迭代或校正的辦法處理。6 一階微分方程組與高階微分方程的數(shù)值解法6.1 一階微分方程組的數(shù)值解法設(shè)有一階微分方程組的初值問(wèn)題 （24）若記，則

11、初值問(wèn)題（24）可寫(xiě)成如下向量形式 （25）如果向量函數(shù)在區(qū)域： 連續(xù)，且關(guān)于滿(mǎn)足Lipschitz條件，即存在，使得對(duì)，都有 那么問(wèn)題（25）在上存在唯一解。問(wèn)題（25）與（1）形式上完全相同，故對(duì)初值問(wèn)題（1）所建立的各種數(shù)值解法可全部用于求解問(wèn)題（25）。6.2 高階微分方程的數(shù)值解法高階微分方程的初值問(wèn)題可以通過(guò)變量代換化為一階微分方程組初值問(wèn)題。設(shè)有階常微分方程初值問(wèn)題 （26）引入新變量，問(wèn)題（26）就化為一階微分方程初值問(wèn)題 （27）然后用6.1中的數(shù)值方法求解問(wèn)題（27），就可以得到問(wèn)題（26）的數(shù)值解。最后需要指出的是，在化學(xué)工程及自動(dòng)控制等領(lǐng)域中，所涉及的常微分方程組初值問(wèn)

12、題常常是所謂的“剛性”問(wèn)題。具體地說(shuō)，對(duì)一階線性微分方程組 （28）其中為階方陣。若矩陣的特征值滿(mǎn)足關(guān)系 則稱(chēng)方程組（28）為剛性方程組或Stiff方程組，稱(chēng)數(shù)為剛性比。對(duì)剛性方程組，用前面所介紹的方法求解，都會(huì)遇到本質(zhì)上的困難，這是由數(shù)值方法本身的穩(wěn)定性限制所決定的。理論上的分析表明，求解剛性問(wèn)題所選用的數(shù)值方法最好是對(duì)步長(zhǎng)不作任何限制。7 Matlab解法7.1 Matlab數(shù)值解7.1.1 非剛性常微分方程的解法Matlab的工具箱提供了幾個(gè)解非剛性常微分方程的功能函數(shù)，如ode45，ode23，ode113，其中ode45采用四五階RK方法，是解非剛性常微分方程的首選方法，ode23采

13、用二三階RK方法，ode113采用的是多步法，效率一般比ode45高。Matlab的工具箱中沒(méi)有Euler方法的功能函數(shù)。（I）對(duì)簡(jiǎn)單的一階方程的初值問(wèn)題 改進(jìn)的Euler公式為 我們自己編寫(xiě)改進(jìn)的Euler 方法函數(shù)eulerpro.m如下：function x,y=eulerpro(fun,x0,xfinal,y0,n);if nargin5,n=50;endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0;for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fun,x(i),y(i); y2=y(i)+h*feval(fun,x(i+1),y

14、1); y(i+1)=(y1+y2)/2;end例1 用改進(jìn)的Euler方法求解 ，解 編寫(xiě)函數(shù)文件doty.m如下：function f=doty(x,y);f=-2*y+2*x2+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文件定義的微分方程右端的函數(shù)。tspan=t0,tfinal是求解區(qū)間，y0是初值。例2 用RK方法求解，解 同

15、樣地編寫(xiě)函數(shù)文件doty.m如下：function f=doty(x,y);f=-2*y+2*x2+2*x;在Matlab命令窗口輸入：x,y=ode45(doty,0,0.5,1)即可求得數(shù)值解。7.1.2 剛性常微分方程的解法Matlab的工具箱提供了幾個(gè)解剛性常微分方程的功能函數(shù)，如ode15s，ode23s，ode23t，ode23tb，這些函數(shù)的使用同上述非剛性微分方程的功能函數(shù)。7.1.3 高階微分方程解法 例3 考慮初值問(wèn)題 解 （i）如果設(shè)，那么 初值問(wèn)題可以寫(xiě)成的形式，其中。（ii）把一階方程組寫(xiě)成接受兩個(gè)參數(shù)和，返回一個(gè)列向量的M文件F.m： function dy=F(t

16、,y);dy=y(2);y(3);3*y(3)+y(2)*y(1);注意：盡管不一定用到參數(shù)和，M文件必須接受此兩參數(shù)。這里向量必須是列向量。（iii）用Matlab解決此問(wèn)題的函數(shù)形式為T(mén),Y=solver(F,tspan,y0)這里solver為ode45、ode23、ode113，輸入?yún)?shù)F是用M文件定義的常微分方程組，tspan=t0 tfinal是求解區(qū)間，y0是初值列向量。在Matlab命令窗口輸入T,Y=ode45(F,0 1,0;1;-1)就得到上述常微分方程的數(shù)值解。這里Y和時(shí)刻T是一一對(duì)應(yīng)的，Y(:,1)是初值問(wèn)題的解，Y(:,2)是解的導(dǎo)數(shù)，Y(:,3)是解的二階導(dǎo)數(shù)。例

17、4 求 van der Pol 方程 的數(shù)值解，這里是一參數(shù)。解 （i）化成常微分方程組。設(shè)，則有（ii）書(shū)寫(xiě)M文件（對(duì)于）vdp1.m:function dy=vdp1(t,y);dy=y(2);(1-y(1)2)*y(2)-y(1);（iii）調(diào)用Matlab函數(shù)。對(duì)于初值，解為 T,Y=ode45(vdp1,0 20,2;0);（iv）觀察結(jié)果。利用圖形輸出解的結(jié)果：plot(T,Y(:,1),-,T,Y(:,2),-)title(Solution of van der Pol Equation,mu=1);xlabel(time t);ylabel(solution y);legend

18、(y1,y2);例5 van der Pol 方程，（剛性）解 (i)書(shū)寫(xiě)M文件vdp1000.m:function dy=vdp1000(t,y);dy=y(2);1000*(1-y(1)2)*y(2)-y(1);（ii）觀察結(jié)果t,y=ode15s(vdp1000,0 3000,2;0);plot(t,y(:,1),o)title(Solution of van der Pol Equation,mu=1000);xlabel(time t);ylabel(solution y(:,1);7.2 常微分方程的解析解在Matlab中，符號(hào)運(yùn)算工具箱提供了功能強(qiáng)大的求解常微分方程的符號(hào)運(yùn)算命令

19、dsolve。常微分方程在Matlab中按如下規(guī)定重新表達(dá)：符號(hào)D表示對(duì)變量的求導(dǎo)。Dy表示對(duì)變量y求一階導(dǎo)數(shù)，當(dāng)需要求變量的n階導(dǎo)數(shù)時(shí)，用Dn表示，D4y表示對(duì)變量y求4階導(dǎo)數(shù)。由此，常微分方程在Matlab中，將寫(xiě)成D2y+2*Dy=y。7.2.1 求解常微分方程的通解無(wú)初邊值條件的常微分方程的解就是該方程的通解。其使用格式為： dsolve(diff_equation) dsolve( diff_equation，var)式中diff_equation 為待解的常微分，第1種格式將以變量t為自變量進(jìn)行求解，第2種格式則需定義自變量var。例6 試解常微分方程 解 編寫(xiě)程序如下：syms

20、x ydiff_equ=x2+y+(x-2*y)*Dy=0;dsolve(diff_equ,x)7.2.2 求解常微分方程的初邊值問(wèn)題求解帶有初邊值條件的常微分方程的使用格式為： dsolve(diff_equation，condition1，condition2，var)其中condition1，condition2， 即為微分方程的初邊值條件。例7 試求微分方程 ，的解。解 編寫(xiě)程序如下：y=dsolve(D3y-D2y=x,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種格式可以加上初邊值條件，用于具體求解