微分方程數(shù)值解

1、微分方程數(shù)值解及其應(yīng)用緒論自然界中的許多事物的運(yùn)動(dòng)和變化規(guī)律都可以用微分方程來(lái)描述，因此對(duì)工程和科學(xué)技術(shù)中的實(shí)際問(wèn)題的研究中, 常常需要求解微分方程但往往只有少數(shù)較簡(jiǎn)單和典型的微分方程可求出其解析解，在大多數(shù)情況下,只能用近似法求解，數(shù)值解法是一類重要的近似方法本文主要討論一階常微分方程的初值問(wèn)題的數(shù)值解法，探討這些算法在處理來(lái)自生活實(shí)際問(wèn)題中的應(yīng)用，并結(jié)合MATLAB軟件，動(dòng)手編程予以解決微分方程的初值問(wèn)題1.1 預(yù)備知識(shí)在對(duì)生活實(shí)際問(wèn)題的研究中，通常需要考慮一階微分方程的初值問(wèn)題 （1）這里是矩形區(qū)域：上的連續(xù)函數(shù)對(duì)初值問(wèn)題（1）需要考慮以下問(wèn)題：方程是否一定有解呢?若有解，有多少個(gè)解呢?

2、下面給出相關(guān)的概念與定理定義1 條件：矩形區(qū)域：上的連續(xù)函數(shù)若滿足：存在常數(shù)，使得不等式對(duì)所有都成立，則稱在上關(guān)于滿足條件.定理 1 解的存在唯一性定理：設(shè)在區(qū)域上連續(xù)，關(guān)于滿足條件，則對(duì)任意的，常微分方程初值問(wèn)題（1）當(dāng)時(shí)存在唯一的連續(xù)解該定理保證若一個(gè)函數(shù)關(guān)于滿足條件，它所對(duì)應(yīng)的微分方程的初值問(wèn)題就有唯一解.在解的存在唯一性得到保證的前提下，自然要考慮方程的求16解問(wèn)題求解微分方程雖然有多種解析方法，但根據(jù)工程和科學(xué)實(shí)踐問(wèn)題所得到的微分方程往往很復(fù)雜，在很多情況下不能或很難給出解析解，有時(shí)即使能求出形式解，也往往因形式過(guò)于復(fù)雜或計(jì)算量太大而不實(shí)用，因此從實(shí)際問(wèn)題中歸結(jié)出來(lái)的微分方程主要依靠

3、數(shù)值解法 定義 2 微分方程數(shù)值解：對(duì)初值問(wèn)題（1）尋求數(shù)值解就是尋求解在一系列離散節(jié)點(diǎn)上的近似解，相鄰兩個(gè)節(jié)點(diǎn)的間距稱為步長(zhǎng).在一般情況下假定為常數(shù)，這時(shí)節(jié)點(diǎn)為要求微分方程數(shù)值解，首先要建立數(shù)值算法，即對(duì)初值問(wèn)題（1）中的方程離散化，建立求解數(shù)值解法的遞推公式一類是計(jì)算時(shí)只用到前一點(diǎn)的值，稱為單步法；另一類是用到前面點(diǎn)的值稱為步法.對(duì)初值問(wèn)題（1）式的單步法可用一般形式表示為，其中多元函數(shù)與有關(guān)，當(dāng)含有時(shí)，方法是隱式的；若中不含，則為顯式方法，所以顯式單步法可表示為. （2）設(shè)是初值問(wèn)題（1）的準(zhǔn)確解，稱為顯式單步法（2）的局部截?cái)嗾`差. 若存在最大正整數(shù),使顯式單步法（2）式的局部截?cái)嗾`差

4、滿足，則稱（2）式有階精度.1.2幾種常用的數(shù)值解法及其分析、比較1.2.1歐拉法與后退歐拉法1）歐拉法：歐拉曾簡(jiǎn)單地用差分代替微分，即利用公式將初值問(wèn)題（1）離散化，則問(wèn)題（1）可化為， （3）此方法稱為歐拉法.歐拉方法的幾何意義在數(shù)值計(jì)算公式中體現(xiàn)了出來(lái).在平面上，一階微分方程的解稱作它的積分曲線.積分曲線上一點(diǎn)的切線斜率等于函數(shù)，按函數(shù)在平面上建立一個(gè)方向場(chǎng)，那么，積分曲線上每一點(diǎn)的切線方向均與方向場(chǎng)在該點(diǎn)的方向相一致.基于上述幾何解釋，從初始點(diǎn)出發(fā)，先依方向場(chǎng)在該點(diǎn)的方向上推進(jìn)到上一點(diǎn)，再?gòu)囊婪较驁?chǎng)的方向推進(jìn)到上一點(diǎn)，循環(huán)前進(jìn)便作出一條折線，因此歐拉方法又稱為折線法.若初值已知，則由（

5、3）式可逐步算出 為了分析計(jì)算公式的精確度，通?？捎锰├照归_(kāi)將在處展開(kāi)，則有在的前提下，可得歐拉法（3）的誤差為容易看出，歐拉法（3）式具有一階精度.2）向后歐拉方法：如果對(duì)微分方程（1）從到積分，得 ， （4）如果（4）式右端積分用右矩形公式近似，則得到另一個(gè)公式 ， （5）稱為后退歐拉法.值得一提的是:后退歐拉法與歐拉公式有著本質(zhì)的區(qū)別，后者是關(guān)于的直接計(jì)算公式，它是顯式的，而（5）式的右端含有關(guān)于的表達(dá)式，它是隱式的.在利用后退歐拉法時(shí)，我們通常利用迭代法求解，實(shí)質(zhì)就是逐步顯示化.具體迭代過(guò)程如下：首先利用歐拉公式給出迭代初值，把它代入（5）式的右端，使之轉(zhuǎn)化為顯式，直接計(jì)算得.如此反復(fù)

6、進(jìn)行，得 ，則得到后退歐拉法的迭代公式，可以看出，后退歐拉法具有一階精度，且計(jì)算比較麻煩.1.2.2梯形方法為得到比歐拉法精確度高的計(jì)算公式，在等式（4）式右端積分中若用梯形求積公式近似，并用代替，代替，則得， （6）稱其為梯形方法.梯形方法與后退歐拉法一樣，都是隱式單步法，可用迭代法求解，其迭代公式為 . （7）為了分析梯形公式的收斂性，將（6）與（7）式相減，得，因?yàn)闈M足條件，于是有，其中為關(guān)于的常數(shù).如果選取充分小，使得，則當(dāng)時(shí)有，這說(shuō)明迭代過(guò)程（7）式是收斂的4.容易推導(dǎo)得出梯形法（7）式是二階方法.經(jīng)分析，梯形方法雖然提高了精度，但是以增加計(jì)算量為代價(jià)的.從上述的迭代公式可以看出，每

7、迭代一次都要重新計(jì)算的值，而且迭代又要進(jìn)行若干次，計(jì)算相當(dāng)?shù)膹?fù)雜.為此，有沒(méi)有比較簡(jiǎn)便的計(jì)算方法呢？下面給出改進(jìn)的歐拉方法.1.2.3改進(jìn)的歐拉方法由前面的討論可知，梯形法計(jì)算相對(duì)復(fù)雜，現(xiàn)對(duì)上面的梯形法進(jìn)行簡(jiǎn)化，具體方法是只計(jì)算一兩次就轉(zhuǎn)入下一步的計(jì)算，先用歐拉公式（3）求得一個(gè)初步的近似解，稱為預(yù)測(cè)值，再利用公式（6）把它校正一次，這樣建立的預(yù)測(cè)-校正系統(tǒng)通常稱為改進(jìn)的歐拉公式.具體公式如下 （8）改進(jìn)的歐拉法與梯形法一樣，是二階方法.1.2.4 Runge-Kutta方法由前面討論可知，從（4）式可以看出，若要使公式階數(shù)提高，就必須使右端積分的數(shù)值求積公式精度提高，它必然要增加求積積點(diǎn)，為

8、此將（4）式的右端用求積公式表示為 ， （9）一般來(lái)說(shuō)，點(diǎn)數(shù)越多，精度越高，上式右端相當(dāng)于增量函數(shù)，為得到便于計(jì)算的顯式方法，將公式（9）表示為： （10）其中 （11） 這里均為常數(shù). 為加權(quán)因子，為第段斜率，共有段.我們把（10）和（11）稱為級(jí)顯式Runge-Kutta法，簡(jiǎn)稱為R-K方法.下面給出其中最經(jīng)典最常用的一個(gè)公式： （12）Runge-Kutta方法作為一種重要的單步方法，具有很高的實(shí)用價(jià)值，它關(guān)于初值是穩(wěn)定的，其解連續(xù)地依賴于初值，是一類便于應(yīng)用的單步法，為了計(jì)算，只用到前面一步的值即可，因此每步的步長(zhǎng)可以獨(dú)立取定.常用的Runge-Kutta方法精度較高，為了達(dá)到預(yù)定的精

9、度，與歐拉方法與梯形法相比，步長(zhǎng)可取得大些，求解區(qū)間上的總步數(shù)可以少些.但Runge-Kutta方法也有些缺點(diǎn)，比如四階Runge-Kutta方法每算一步需要四次計(jì)算的值，計(jì)算量較大（對(duì)于復(fù)雜的而言）2 數(shù)值方法的應(yīng)用實(shí)例5-9例1 對(duì)于初值問(wèn)題，分別用歐拉法、改進(jìn)的歐拉法，梯形法求的近似值.解：易得該方程的解析解，為比較，將按不同數(shù)值計(jì)算方法所得結(jié)果列表如下：表 1 三種不同方法的數(shù)值結(jié)果 歐拉法改進(jìn)的歐拉法梯形法0.2-1110.109.7656E-0045.4994E-0050.012.6561-0054.6223-0054.5026-0050.0014.3717-0054.5408-0

10、054.5396-0050.00014.5173-0054.5400-0054.5400-005 圖 1 三種不同方法數(shù)值解與精確解的誤差曲線從表1中可以看出：當(dāng)時(shí)，三種方法均不穩(wěn)定，計(jì)算結(jié)果嚴(yán)重偏離精確值；時(shí)，改進(jìn)后的歐拉和梯形法均穩(wěn)定，但歐拉法效果很差；當(dāng)時(shí)，三種方法均穩(wěn)定，但精確度有區(qū)別可以看出，越小，計(jì)算結(jié)果越好，要想計(jì)算結(jié)果充分接近于解析解還須取較小的值圖1反映的步長(zhǎng)時(shí)，三種數(shù)值方法的所得數(shù)值解與解析解在區(qū)間的誤差曲線，由圖可知，在步長(zhǎng)相同的情況下，梯形法的精確度略高于改進(jìn)的歐拉法；改進(jìn)的歐拉法和梯形法精確度都明顯高于歐拉法例2 用歐拉法、改進(jìn)的歐拉法和Runge-Kutta法求解初

11、值問(wèn)題并比較三種方法的結(jié)果 解：方程為的伯努利方程，可求得解析解為現(xiàn)用MATLAB軟件編程，用題目要求的方法求解，可得如下圖示結(jié)果：圖2 （a）步長(zhǎng)為0.2時(shí)R-K法和解析解比較圖2 （b） 步長(zhǎng)為0.2時(shí)改進(jìn)的Euler法和解析解比較圖2 （c）步長(zhǎng)為0.2時(shí)歐拉法和解析解比較上圖2(a)，(b)，(c)描述的是步長(zhǎng)為0.2時(shí)，用歐拉法、改進(jìn)的歐拉法，Runge-Kutta法求解方程所得的數(shù)值解與解析解之間的對(duì)比圖由圖可知，Runge-Kutta法所得數(shù)值解曲線和解析解曲線吻合的很好，改進(jìn)的歐拉法和歐拉法隨著計(jì)算的進(jìn)行，數(shù)值解和解析解之間誤差逐步增大，但改進(jìn)的歐拉法效果要好于歐拉法圖3 (a

12、) 步長(zhǎng)為0.1時(shí)Euler法和解析解比較 圖3 (b) 步長(zhǎng)為0.1時(shí)改進(jìn)的Euler法和解析解比較 圖3 (c) 步長(zhǎng)為0.1時(shí)Runge-Kutta法和解析解比較上圖3 (a)，(b)，(c)描述的是步長(zhǎng)為0.1時(shí)，用歐拉法、改進(jìn)的歐拉法，Runge-Kutta法求解方程所得的數(shù)值解與解析解之間的對(duì)比圖由圖可知，改進(jìn)的歐拉法和Runge-Kutta法所得數(shù)值解曲線和解析解曲線吻合的很好，而歐拉法隨著計(jì)算的進(jìn)行數(shù)值解和解析解之間誤差逐步增大相應(yīng)的程序如下：主程序x=0:0.2:1;jxj=exp(2*x).*(1./exp(4*x) + (2*x)./exp(4*x).(1/2);y=Eu

13、ler(ff,0,1,0.2,1);gy=geuler(ff,0,1,0.2,1);Ry=RK(ff,0,1,0.2,1);figure(1);plot(x,jxj,x,Ry,'*');figure(2);plot(x,jxj,x,gy,'*');figure(3);plot(x,jxj,x,y,'*')歐拉法程序function y=Euler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:ny(i+1)=y(i)+h*feval(f,x(i),y(i);end改進(jìn)的歐拉

14、法程序function gy=geuler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i);yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endgy=y;Runge-Kutta法程序function Ry=RK(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nk1=feval(f,x(i),y(i);k2=feval(f,x(i)+h/2,y(i)

15、+h*k1/2);k3=feval(f,x(i)+h/2,y(i)+h*k2/2);k4=feval(f,x(i+1),y(i)+h*k3);y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;endRy=y;3 微分方程數(shù)值解法在實(shí)際生活中的應(yīng)用3.1應(yīng)用實(shí)例:耐用消費(fèi)新產(chǎn)品的銷售規(guī)律模型 一種新產(chǎn)品進(jìn)入市場(chǎng)以后，常常會(huì)經(jīng)歷銷售量首先慢慢增加然后逐漸慢慢下降的一個(gè)過(guò)程，由此給出的時(shí)間銷售坐標(biāo)系下的曲線稱為產(chǎn)品的生命曲線,它的形狀是鐘形的.不過(guò)對(duì)于較耐用的消費(fèi)品，情況會(huì)有所不一樣，它的生命曲線會(huì)在開(kāi)始有個(gè)小小的高峰，之后是段比較平坦的曲線，先下降，之后再上升，然后達(dá)到高峰，因此

16、它是雙峰型的曲線. 如何理解這種和傳統(tǒng)產(chǎn)品的生命曲線理論相沖突的現(xiàn)象？澳大利亞學(xué)者斯蒂芬斯與莫賽觀察到的購(gòu)買耐用消費(fèi)產(chǎn)品的人大概可分為兩類：其一是非常善于接受新的事物，稱其為“創(chuàng)新型”消費(fèi)者，他們會(huì)經(jīng)常從產(chǎn)品廣告，制造商給出產(chǎn)品的說(shuō)明書與商店樣品中了解了產(chǎn)品功能與性能之后，再?zèng)Q定其否購(gòu)買；其二是消費(fèi)者相對(duì)保守，稱其為“模仿型”消費(fèi)者，他們往往會(huì)根據(jù)大部分已經(jīng)購(gòu)買產(chǎn)品的消費(fèi)者實(shí)際使用的經(jīng)驗(yàn)而提供的信息來(lái)決定是否購(gòu)買產(chǎn)品，下面的例子經(jīng)過(guò)分析，建立相應(yīng)的數(shù)學(xué)模型，對(duì)這種現(xiàn)象給出了科學(xué)解釋.3.1.1 模型的建立 將顧客獲得信息大致可分成兩類，其一稱之為“搜集型”，源自于產(chǎn)告說(shuō)明、廣告，“創(chuàng)新型”顧客

17、獲得這類消息后就可作出其是否購(gòu)買；第二類稱之為“體驗(yàn)型”，即消費(fèi)者使用之后會(huì)獲得真實(shí)體驗(yàn)，常常以口頭的形式散播，“模仿型”類顧客會(huì)在獲得這種信息之后才能決定是否購(gòu)買. 設(shè)是潛在用戶的總數(shù)，與分別為 “創(chuàng)新型”與“模仿型”的人數(shù).再設(shè)是時(shí)刻已經(jīng)購(gòu)買的商品顧客的人數(shù)，而與分別為“創(chuàng)新型”與“模仿型”的顧客的人數(shù),再設(shè)是時(shí)刻中已獲得“搜集型”信息的人數(shù)，由于該部分的信息可直接由外部獲得，或者可從已獲得該類信息的人群中獲得，因此類似巴斯模型，從而建立如下方程： ， 獲得“搜集型”信息的“創(chuàng)新型”消費(fèi)者決定其是否購(gòu)買的行為，滿足如下方程： 對(duì)于“模仿型”的顧客，可從已經(jīng)購(gòu)買該產(chǎn)品“創(chuàng)新型”或者“模仿型”

18、的顧客中獲取信息，于是有 在這里，忽略消費(fèi)者購(gòu)買產(chǎn)品后需一段短暫的使用后才會(huì)散播體驗(yàn)信息的滯后作用綜上，斯蒂芬斯莫賽模型是一常微分方程組的初值問(wèn)題： 而為時(shí)刻購(gòu)買該商品的總?cè)藬?shù).3.1.2 模型的求解對(duì)于斯蒂芬斯莫賽模型中的解析解則不能求出，于是可以用四階公式求得，且在它的精度要求達(dá)到很高情形下求出.用MATLAB軟件求解，相應(yīng)的程序及結(jié)果如下function RK=RKFC(fc,a,b,h,y0)n=(b-a)/h;x=a:h:b;m=length(y0);y=zeros(n+1,m);y(1,:)=y0;for i=1:nk1=feval(fc,x(i),y(i,:);k2=feval(

19、fc,x(i)+h/2,y(i,:)+h*k1/2);k3=feval(fc,x(i)+h/2,y(i,:)+h*k2/2);k4=feval(fc,x(i+1),y(i,:)+h*k3);y(i+1,:)=y(i,:)+h*(k1+2*k2+2*k3+k4)/6;endRK=y;function f=FC(x,y)k1=50; k2=70;al=0.0013;be=0.0013;ga=0.0015;f=(k1-y(1)*(al+be*y(1),ga*(k2-y(2)*(y(1)+y(2);x=0:0.3:24;RK=RKFC(FC,0,24,0.3,0,0)figure(1)；plot(x,RK(:,1),'+',x,RK(:,2),'*')；legend('N1(t)','N2(t)',2)figure(2)；plot(x,RK(:,1)+RK(:,2),'+')；legend('N(t)',2) 圖 4 與時(shí)間關(guān)系圖 圖 5 與0,25時(shí)間段關(guān)系圖由此例