第十五章 常微分方程的解法

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

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

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

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

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

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

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

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

9、式類似、但更復(fù)雜的計算，得到使局部誤差的11個方程。取既滿足這些方程、又較簡單的一組，可得 （18）這就是常用的4階龍格庫塔方法（簡稱RK方法）。§5 線性多步法以上所介紹的各種數(shù)值解法都是單步法，這是因為它們在計算時，都只用到前一步的值，單步法的一般形式是 （19）其中稱為增量函數(shù)，例如Euler方法的增量函數(shù)為，改進Euler法的增量函數(shù)為如何通過較多地利用前面的已知信息，如，來構(gòu)造高精度的算法計算，這就是多步法的基本思想。經(jīng)常使用的是線性多步法。讓我們試驗一下，即利用計算的情況。從用數(shù)值積分方法離散化方程的（4）式 出發(fā)，記，式中被積函數(shù)用二節(jié)點，的插值公式得到（因，所以是外插

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

11、 一階微分方程組的數(shù)值解法設(shè)有一階微分方程組的初值問題 （24）若記，則初值問題（24）可寫成如下向量形式 （25）如果向量函數(shù)在區(qū)域： 連續(xù)，且關(guān)于滿足Lipschitz條件，即存在，使得對，都有 那么問題（25）在上存在唯一解。問題（25）與（1）形式上完全相同，故對初值問題（1）所建立的各種數(shù)值解法可全部用于求解問題（25）。6.2 高階微分方程的數(shù)值解法高階微分方程的初值問題可以通過變量代換化為一階微分方程組初值問題。設(shè)有階常微分方程初值問題 （26）引入新變量，問題（26）就化為一階微分方程初值問題 （27）然后用6.1中的數(shù)值方法求解問題（27），就可以得到問題（26）的數(shù)值解。最

12、后需要指出的是，在化學(xué)工程及自動控制等領(lǐng)域中，所涉及的常微分方程組初值問題常常是所謂的“剛性”問題。具體地說，對一階線性微分方程組 （28）其中為階方陣。若矩陣的特征值滿足關(guān)系 則稱方程組（28）為剛性方程組或Stiff方程組，稱數(shù)為剛性比。對剛性方程組，用前面所介紹的方法求解，都會遇到本質(zhì)上的困難，這是由數(shù)值方法本身的穩(wěn)定性限制所決定的。理論上的分析表明，求解剛性問題所選用的數(shù)值方法最好是對步長不作任何限制。§7 Matlab解法7.1 Matlab數(shù)值解7.1.1 非剛性常微分方程的解法Matlab的工具箱提供了幾個解非剛性常微分方程的功能函數(shù)，如ode45，ode23，ode1

13、13，其中ode45采用四五階RK方法，是解非剛性常微分方程的首選方法，ode23采用二三階RK方法，ode113采用的是多步法，效率一般比ode45高。Matlab的工具箱中沒有Euler方法的功能函數(shù)。（I）對簡單的一階方程的初值問題 改進的Euler公式為 我們自己編寫改進的Euler 方法函數(shù)eulerpro.m如下：function x,y=eulerpro(fun,x0,xfinal,y0,n);if nargin<5,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*fev

14、al(fun,x(i),y(i); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2;end例1 用改進的Euler方法求解 ，解 編寫函數(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，

15、輸入?yún)?shù) F 是用M文件定義的微分方程右端的函數(shù)。tspan=t0,tfinal是求解區(qū)間，y0是初值。例2 用RK方法求解，解 同樣地編寫函數(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的工具箱提供了幾個解剛性常微分方程的功能函數(shù)，如ode15s，ode23s，ode23t，ode23tb，這些函數(shù)的使用同上述非剛性微分方程的功能函數(shù)。7.1.3 高階微分方程解法 例3 考慮初值問題 解

16、（i）如果設(shè)，那么 初值問題可以寫成的形式，其中。（ii）把一階方程組寫成接受兩個參數(shù)和，返回一個列向量的M文件F.m： function dy=F(t,y);dy=y(2);y(3);3*y(3)+y(2)*y(1);注意：盡管不一定用到參數(shù)和，M文件必須接受此兩參數(shù)。這里向量必須是列向量。（iii）用Matlab解決此問題的函數(shù)形式為T,Y=solver('F',tspan,y0)這里solver為ode45、ode23、ode113，輸入?yún)?shù)F是用M文件定義的常微分方程組，tspan=t0 tfinal是求解區(qū)間，y0是初值列向量。在Matlab命令窗口輸入T,Y=ode

17、45('F',0 1,0;1;-1)就得到上述常微分方程的數(shù)值解。這里Y和時刻T是一一對應(yīng)的，Y(:,1)是初值問題的解，Y(:,2)是解的導(dǎo)數(shù)，Y(:,3)是解的二階導(dǎo)數(shù)。例4 求 van der Pol 方程 的數(shù)值解，這里是一參數(shù)。解 （i）化成常微分方程組。設(shè)，則有（ii）書寫M文件（對于）vdp1.m:function dy=vdp1(t,y);dy=y(2);(1-y(1)2)*y(2)-y(1);（iii）調(diào)用Matlab函數(shù)。對于初值，解為 T,Y=ode45('vdp1',0 20,2;0);（iv）觀察結(jié)果。利用圖形輸出解的結(jié)果：plot(T

18、,Y(:,1),'-',T,Y(:,2),'-')title('Solution of van der Pol Equation,mu=1');xlabel('time t');ylabel('solution y');legend('y1','y2');例5 van der Pol 方程，（剛性）解 (i)書寫M文件vdp1000.m:function dy=vdp1000(t,y);dy=y(2);1000*(1-y(1)2)*y(2)-y(1);（ii）觀察結(jié)果t,y=ode1

19、5s('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中，符號運算工具箱提供了功能強大的求解常微分方程的符號運算命令dsolve。常微分方程在Matlab中按如下規(guī)定重新表達：符號D表示對變量的求導(dǎo)。Dy表示對變量y求一階導(dǎo)數(shù)，當(dāng)需要求變量的n階導(dǎo)數(shù)時，用Dn表

20、示，D4y表示對變量y求4階導(dǎo)數(shù)。由此，常微分方程在Matlab中，將寫成'D2y+2*Dy=y'。7.2.1 求解常微分方程的通解無初邊值條件的常微分方程的解就是該方程的通解。其使用格式為： dsolve('diff_equation') dsolve(' diff_equation'，'var')式中diff_equation 為待解的常微分，第1種格式將以變量t為自變量進行求解，第2種格式則需定義自變量var。例6 試解常微分方程 解 編寫程序如下：syms x ydiff_equ='x2+y+(x-2*y)*Dy=

21、0'dsolve(diff_equ,'x')7.2.2 求解常微分方程的初邊值問題求解帶有初邊值條件的常微分方程的使用格式為： dsolve('diff_equation'，'condition1，condition2，'，'var')其中condition1，condition2， 即為微分方程的初邊值條件。例7 試求微分方程 ，的解。解 編寫程序如下：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種格式可以加上初邊值條件，用于具體求解。例8