差分方程方法

1、ii=1第四章差分方程方法在實際中，許多問題所研究的變量都是離散的形式，所建立的數(shù)學(xué)模型也是離散的，譬如，像政治、經(jīng)濟(jì)和社會等領(lǐng)域中的實際問題。有些時候，即使所建立的數(shù)學(xué)模型是連續(xù)形式，例如像常見的微分方程模型、積分方程模型等等，但是，往往都需要用計算機(jī)求數(shù)值解。這就需要將連續(xù)變量在一定條件下進(jìn)行離散化，從而將連續(xù)型模型轉(zhuǎn)化為離散型模型，因此，最后都?xì)w結(jié)為求解離散形式的差分方程解的問題。關(guān)于差分方程理論和求解方法在數(shù)學(xué)建模和解決實際問題的過程中起著重要作用。下面就不同類型的差分方程進(jìn)行討論。所謂的差分方程是指：對于一個數(shù)列，把數(shù)n列中的前n+1項x(i二0,1,2,n)關(guān)聯(lián)起來所得到的方程。i

2、41常系數(shù)線性差分方程常系數(shù)線性齊次差分方程常系數(shù)線性齊次差分方程的一般形式為x+ax+ax+.+ax=0(4.1)n1n-12n-2kn-k其中k為差分方程的階數(shù)，a(i=1,2,k)為差分方程的系數(shù)，且a豐0(kn)。對應(yīng)的ik代數(shù)方程九k+a九k-1+a九k-2+a=0(4.2)12k稱為差分方程的(4.1)的特征方程，其特征方程的根稱為特征根。常系數(shù)線性齊次差分方程的解主要是由相應(yīng)的特征根的不同情況有不同的形式。下面分別就特征根為單根、重根和復(fù)根的情況給出差分方程解的形式。特征根為單根設(shè)差分方程(4.1)有k個單特征根九，九,九，,九，則差分方程(4.1)的通解為TOC o 1-5 h

3、 z123kx=C九n+C九n+C九n，n1122kk其中c,c，,c為任意常數(shù)，且當(dāng)給定初始條件12k HYPERLINK l bookmark10 o Current Document x=X(0)(i=1,2,k)(4.3)ii時，可以唯一確定一個特解。特征根為重根設(shè)差分方程(4.1)有l(wèi)個相異的特征根九，九，九，,九Glk)重數(shù)分別為123lm1，m2，mi且m=k則差分方程()的通解為x二區(qū)cni-1九n+芳cni-1九nF瓦Cni-1九nn1i12i2lili=1i=1i=1同樣的，3.由給定的初始條件（4.3）可以唯一確定一個特解。特征根為復(fù)根設(shè)差分方程（4.1）的特征根為一對共

4、軛復(fù)根九，九=aip和相異的k-2個單根12九，九,九，則差分方程的通解為34kx=cpncosn0+cpnsinn0+c心+c心HFc九n,n123344kki其中i其中p=尹2+p2,9=arctanP同樣由給定的初始條件（4.3）可以惟一確定一個特解。另外，對于有多個共軛復(fù)根和相異實根，或共軛復(fù)根和重根的情況，都可以類似地給出差分方程解的形式。412常系數(shù)線性非齊次差分方程常系數(shù)線性非齊次差分方程的一般形式為(4.4)x+ax+axFFax(4.4)n1n-12n-2kn-k其中k為差分方程的階數(shù)，川二1,2，,為差分方程的系數(shù)，ak豐-If（n）為已知函數(shù)。在差分方程（4.4）中，令f

5、（）=0,所得方程x+ax+axFFax=0（4.5）n1n-12n-2kn-k稱為非齊次差分方程（4.4）對應(yīng)的齊次差分方程，即與差分方程（4.1）的形式相同。求解非齊次差分方程通解的一般方法為首先求對應(yīng)的齊次差分方程（4.5）的通解x*，然后求非齊次差分方程（4.4）的一個n特解x（o），貝ynx=x*+x（o）nnn為非齊次差分方程（4.4）的特解。關(guān)于求x*的方法同求差分方程（4.1）的方法相同。對于求非齊次方程（4.4）的特解nx（0）的方法，可以用觀察法確定，也可以根據(jù)f（n）的特性用待定系數(shù)法確定，具體方法可n參照常系數(shù)線性非齊次微分方程求特解的方法。4.2差分方程的平衡點及其穩(wěn)

6、定性一般來說，差分方程的求解是困難，實際中往往不需要求出差分方程的一般解，而只需要研究它的平衡點及其穩(wěn)定性即可。一階線性常系數(shù)差分方程一階線性常系數(shù)差分方程的一般形式為x+ax=b,k=0,1,2,k+1k其中a,b為常數(shù)，它的平衡點由代數(shù)方程x+ax二b求解得到，不妨記為x*.如果limx=x*，則稱平衡點x*是穩(wěn)定的，否則是不穩(wěn)定的。kT8*為了便于研究平衡點x*的穩(wěn)定性問題，一般將其轉(zhuǎn)化為求方程xk+1+叮0的平衡點x*=0的穩(wěn)定性問題。事實上，由x+ax-0k+1k可以解得x=Cax,k0于是x*=0是穩(wěn)定的平衡點的充要條件是：|a|Y1.4.2.2一階線性常系數(shù)差分方程組一階線性常系

7、數(shù)齊次差分方程組的一般形式為x(k+1)+Ax(k)=B,k=0,1,2,.其中x(k)為口維向量，A為nxn階常數(shù)矩陣。它的平衡點x*=0是穩(wěn)定的充要條件是A的所有特征根都有卜|1(i=l,2,n).i對于一階線性常系數(shù)非齊次差分方程組x(k+1)+Ax(k)=B，k=0,1,2,.的情況同樣給出4.2.3二階線性常系數(shù)差分方程二階線性常系數(shù)齊次差分方程的一般形式為x+ax+ax=0,k=0,1,2,k+21k+12k其中a,a為常數(shù)，其平衡點x*=0是穩(wěn)定的充要條件是特征方程九2+a九+a二0,1212的根九,九2滿足卜J1，卜J1。對于一般的x+ax+ax二0的平衡點的穩(wěn)定性問題同樣給出

8、。類似地，也可直k+21k+12k接推廣到n階線性差分方程的情況。4.2.4一階非線性差分方程一階非線性方程的一般形式為x=f(x)，k=0,1,2,，k+1k其中f為已知函數(shù)，其平衡點定義為方程x=f(x)的解x*。-f(x*)(-x-f(x*)(-x*)+f(*),kx則x*也是一階線性萬程xk+1已f(x*)(-x*)+f(x*)的平衡點，故此，平穩(wěn)衡點x*穩(wěn)定kk+1的充要條件是|f(x*)1。4.3連續(xù)模型的差分方法微分的差分方程已知f(x)在點x處的函數(shù)值f(x)k=0,1,.,n+1)，且kka=xx.x=b，試求函數(shù)的導(dǎo)數(shù)值f(x),(k=1,2,.n)。01n+1k根據(jù)導(dǎo)數(shù)的

9、定義，用差商代替微商，則有下面的差分公式。向前差：f(、f(x)f(x)(k12)f(x)沁屮k(k=1,2,.,n),kx-xk+1k向后差：f(xk)-f向后差：f(xk)-f(x)-f(x)kx-xkk-1(k=1,2,.,n),中心差：f(xk)-f(x)-f(x)(k=1,2,.,n),kk-1定積分的差分方法已知函數(shù)f(X)在點x處的函數(shù)值f(x),(k=1,2,.,n),且在a,b上可積，試求函kk數(shù)在a,b上的積分值ibf(x)dx。a根據(jù)定積分的定義，則有一般的求積公式ibf(x)dxibf(x)dxakkk=0其中A為求積系數(shù)，它與x的選取方法有關(guān)。取不同的求積系數(shù)，可以得

10、不同的求積公式。kkba對于等距節(jié)點x=a+kh(k=0,1,.,n)，其中步長h=為很小的數(shù)，則有如kn下的求積公式。復(fù)化矩陣公式；2)復(fù)化梯形矩陣；2)復(fù)化梯形矩陣；Jbf(xJbf(x)dxqh藝f(x)+f(x)=a2kk+1k=0復(fù)化辛普森(Simpson)矩陣公式；ibf(x)dxQa+2刃f(x)+f(b)kk=1呢f(xk)+4f(xk)+f(xk+1)k=0k+2=f(a)+邏f(xk+1盪f(f(b)k=0k2k=1其中xk+其中xk+12(4)=2(xk+x)為子區(qū)間x,x(k=0,1,.,n1)的中點。k+1kk+1復(fù)化柯特斯(Cotes)公式;(1)a+k+hI2丿J

11、ibf(x)dxqh藝fak=032藝f(xk+3)+藝f(7f(b)k=0k4k=1丄.32藝f(xk+3)+藝f(7f(b)k=0k4k=1丄.f(x丄)+k+k+4k=02k=0其中x,x,x為子空間x,x(k二0,1,.,n-1)中的四等分點。,Jkk+1k+k+k+2444.3.3常微分方程的差分方法一階常微分方程的差分方法設(shè)一階常微分方程的定解問題為y=f(x,y),y(x)=y,00其中函數(shù)f(x,y)關(guān)于y滿足李普希茲條件，即保證問題解的存在唯一性?，F(xiàn)在的問題是求方程在一系列節(jié)點xx.x.處的近似數(shù)值解12ny,y,y,不妨假設(shè)步長為h二x-x為常數(shù)。在此，我們根據(jù)微分的差分方

12、法，12nn+1n即用差商來近似代替微商，再利用“步進(jìn)式”方法，可以給出求解問題(4-6)的差分方法。單步歐拉(Euler)公式用差商竺*丄近似代替y(x)二f(x,y(x)中的導(dǎo)數(shù)，則可以得差分公式hnnny二y+hf(x,y),n二0,1,2,.n+1nnn其精度為O(h2)階的。兩步歐拉公式用差商蘭化+1二近似代替y(x)二f(x,y(x)中的導(dǎo)數(shù)，則可得差分公式2hnnny二y+2hf(x,y),n二0,1,2,.n+1n-1nn兩步法需要用到前兩步的方信息，一般不能自行起步，需先用單步方法求出y，其精度是O(h2)階的。梯形公式對于方程y二f(x,y)的兩邊在x,x上求積分得nn+1

13、y(x)二y(x)二y(x)+ixn+if(x,y(x)dx.n+1nnxn利用積分的差分方法中梯形公式求解積分,y(x)n+1h,y(x)n+1x+1f(x,y(x)dx沁-f(x,y(x)+f(x2nnn+1hy(x)-y(x)+-f(x,y(x)+f(x,y(x)n+1n2nnn+1n+1離散化即可得到微分方程的梯形差分公式-y=y+-f(x,y)+f(x,y),N=0,1,2,n+1N2nNn+1N+1這是一個隱式格式，計算量大，一般不單獨使用。其鏡的也是O(-3)階的。改進(jìn)的歐拉公式由于單步歐拉公式色精度低，但計算量?。痪匦喂骄雀?，但是計算量大，為此我們綜合運用這兩種方法舅老爺?shù)?/p>

14、到改進(jìn)的歐拉公式，其精度為O3)階的。預(yù)報：yn+yn+1,y),n=0,1,nn校正：yn+1,ynn,yyn+1,ynn,yxin+1-n+1丿,n=0,1,或?qū)懗善骄问?；cnncnn+1py=三(y+y),n=0,1,2,.n+12pc龍格庫塔法龍格庫塔方法的基本思想：y(x)-y(x)對于微分方程的定解問題(4.6),考慮差商八n+1：丿n，根據(jù)阿格朗日微分中值-定理可得gx，n+1y(x)=y(x)+-Y*?，F(xiàn)在n+1ny(x)gx，n+1y(x)=y(x)+-Y*?，F(xiàn)在n+1nn+1nn記Y*二f(g,y(g)的平均變化率,nn+1的問題只要找到尋找一種好的計算Y*的近似方法。

15、如果取Y*沁f(x,y)=Y，則就是歐拉公式。nn1如果取Y*沁f(x,y)+fC,y力=Y，則相應(yīng)的就是改進(jìn)的歐拉公式。2nnn+1n+12現(xiàn)在，我們?nèi)個點(x+ah,y+Ph)gEx,x(=1,2,m)，用f在這m個nininn+1點的函數(shù)值的加權(quán)平均作為Y*的近似值，即33藝藝f(x+ah,y+i=10h)ni其中為權(quán)系數(shù)。則有in+1+h藝f(xii=1+ah,y+0n+1+h藝f(xii=1+ah,y+0h)ini4.7)其中a，i為待定系數(shù)。i實際上，適當(dāng)選擇a，i0，使得公式有更咼的精度，這是龍格庫塔方法的思i想。二階龍格-庫塔公式在lx,x內(nèi)取中點xnn+1=xn+211+-

16、h，則可取巴=0,篤=1，a=卩=-代人(4.7)22式得到二階龍格-庫塔公式，其精度為O(h3階y=y+hY,n+1n2Y=f(x,y),1nnhhjn2n21丿y2=fX七,y+：Y1,n=0,1,2,三階龍格-庫塔公式：在lx,x內(nèi)任取二點x=x+ph，x=x+qh，(0pq1)類似的方法nn+1n+pnn+qn可得到三階的龍格-庫塔公式其精度是O(4)階的常用的是三階的情況。y=y+-(Y+4Y+Y)n+1n6123=f(x,y),nn/hh)=fIx+,y+Y(n2n21丿=f(x+h,y+hCY+2Y)nn12,n=0,1,2,四階龍格庫塔公式：類似的方法可以得到四階龍格-庫塔公式

17、，其精度是O(h5y=y+-(Y+2Y+2Y+Y)TOC o 1-5 h zn+1n61234Y=f(x,y),nn“/hh“)小-宀Y=fx+,y+Y,n=0,1,2,丿In2n2i丿+h+hY22Y=f(x,y+hY)4n+1n3一階常微方程組的差分方法將前面的單個方程中的變量和函數(shù)視為向量，相應(yīng)的差分方法即可用于由多個方程組的一階方程組的情形。對于二個方程的方程組00(4.8)y=f(x.y,zz=g(X,y,z00(4.8)設(shè)以y,z表示函數(shù)在節(jié)點x=x+nh,n=1,2,上的近似解，則有改進(jìn)的歐拉公式:nnn0預(yù)報:校正:y=yn+1nn+1hy,預(yù)報:校正:y=yn+1nn+1hy

18、,=y+2n+1n2,y,znnn,y,znnnf(x,y,zn+1,yn+1,zn+1g(xg(x,y,znnn(gX,yn+1n+1z_n+1丿四階龍格-庫塔公式y(tǒng)=y+-(Y+2Y+2Y+Y)n+1n61234C1c,n=0,1,2,z=z+-(Z+2Z+2Z+Z)n+1n61234其中Y=f(xZ=gY=f(xZ=g(x,y,z;nnn,y,znnn);hhh、 HYPERLINK l bookmark125 o Current Document n2n21n21丿hhh、n2n21n21丿hhh、n2n22n22丿hhh、n2n22n22/Y=fa,y+hY,z+hZ)4n+1n3n

19、3Z=gvc,y+hY,z+hZ丿4n+1n3n3其他的公式也都可以類似得到，即相當(dāng)于同時求解多個一階方程，從方法上沒有本質(zhì)的差別。3.高階常微分方程的差分方法對于某些高階方法的定解問題，原則上可以轉(zhuǎn)化為一階方程組來求解。譬如，對于如下的二階微分方程的定解問題Y=fx+,y+Y,z+Z2Z=gx+,y+Y,z+Z2ccfVY=fx+,y+Y,z+-Z3Z=gx+,y+-Y,z+-Z3丄)y=fW,y,y丿y(x)=y,y(x)=y0000若令z二y，則可以化為一階方程組的定解問題z=fC,y,z)(4.9)Vy=z,y(x)=y,y(x)=y0000(4.9)實際上,(4.9)式可以視為(4.

20、8)式的特例，類似地可以得到相應(yīng)的求解差分公式4.4最優(yōu)捕魚問題4.3.1問題的提出假設(shè)鯷魚可分為4個年齡組：稱1、2、3、4齡魚。各年齡組每條魚的平均重量分別為5.07，11.55，17.86，22.99(g)；各年齡組魚的自然死亡率均為0.8(1/年)；這種魚為季節(jié)性集中產(chǎn)卵繁殖，產(chǎn)卵孵化期為每年的最后4個月，平均每條4齡魚的產(chǎn)卵量為1.109*105(個)，3齡魚的產(chǎn)卵量為這個數(shù)的一半，2齡和1齡魚不產(chǎn)卵。卵孵化并成活為1.22x10111齡魚，成活率(1齡魚條數(shù)與產(chǎn)卵量n之比)為.22x1011+n漁業(yè)部門規(guī)定，每年只允許在產(chǎn)卵孵化期前的8個月內(nèi)進(jìn)行捕撈作業(yè)。如果每年投入的捕撈能力固定

21、不變，即固定努力量捕撈，這時單位時間捕撈量將與各年齡組魚群條數(shù)成正比比例系數(shù)稱為捕撈強(qiáng)度系數(shù)。通常使用13mm網(wǎng)眼的拉網(wǎng)，這種網(wǎng)只能捕撈3、4齡魚，其兩個捕撈系數(shù)之比為0.42：1.要解決的問題是：建立數(shù)學(xué)模型，分析如何實現(xiàn)可持續(xù)性捕撈(即每年開始捕撈時漁場中各年齡組魚群條數(shù)不變)，并且在此前提下得到最高的年收獲量(總質(zhì)量)。模型的假設(shè)與符號說明模型的假設(shè)只考慮魚的繁殖和捕撈的變化，不考慮魚群遷入與遷出；各齡魚在一年的任何時間都會發(fā)生自然死亡；所有魚都在每年最后四個月內(nèi)(后1/3年)完成產(chǎn)卵孵化的過程，成活的幼魚在下一年初成為1齡魚；產(chǎn)卵發(fā)生于后四個月之初，產(chǎn)卵魚的自然死亡發(fā)生與產(chǎn)卵之后；相鄰

22、兩個年齡組的魚群在相鄰兩年之間變化是連續(xù)的，即第k年底i齡魚的條數(shù)等于第k+1年初i+1齡魚的條數(shù)；4齡以上的魚全部死亡；采用固定努力量捕撈意味著捕撈的速率正比于捕撈時各齡魚群的條數(shù)，比例系數(shù)為捕撈強(qiáng)度系數(shù)。符號的說明用x(t)表示t時刻(年)i齡魚的條數(shù)；r表示魚的平均自然死亡率，即r=0.8；fiiA表示齡i魚的產(chǎn)卵數(shù)，即(f,f,f,f)=(0,0,A)，A=1.109x105；w表示齡i魚群TOC o 1-5 h z12342i的捕撈強(qiáng)度系數(shù)，即(w,w,w,w)二(5.07,11.55,17.86,22.99)；q表示i齡魚群的捕撈1234i-2強(qiáng)度系數(shù)，即(q,q,q,q)二(0,

23、0,0.42e,e)，e為捕撈努力量；t=表示產(chǎn)卵開始的12343Y月份；Y表示i齡魚的捕撈量；C表示i齡魚的捕撈率，即C。iiixi模型的建立與求解無捕撈時魚群的自然增長模型由假設(shè)(1)和(2)得dx(/dx(/)idt=-rxC),ii=1,2,3,4；ktk+1，k=0,1,2,又由假設(shè)(3)和(4)得x(kx(k+1)=11.22*10111.22*10ii+xk+-x(k+tA=x(k+t+Ax(k+10(丿23(丿4(丿由假設(shè)(5)和(6)得x(0)=x,i1limx(t)tk+11x(Kx(0)=x,i1limx(t)tk+11I+1i+(I=1,2,3;k=0,1,)2。固定努

24、力量捕撈魚群的增長和捕撈模型由假設(shè)知，捕撈期為ktk+1，則有則有則有dxC)-i=-rxdtI(t)-q(E)xC),dxC)-i=-rxdtI(t)-q(E)xC),kIItk+t-(4.10)dx(t)i-dt-rx(t)k+ttk+1(4.11)x(0)=xI1x(k+1)=x(k+1)，I=1,2,3I+1I+(4.12)x(k+1)=11.22*1011(1.22*1011+xk+1(4.13)(A(A(-)k+t=xk+1+Ax1k+1V丿23V丿4V丿x0k（1）魚群的增長規(guī)律求解方程（1）魚群的增長規(guī)律求解方程（4.10）和（4.11），并利用連續(xù)條件（4.12）式可得x(k

25、+1)=sl(E)X(k),I=1,2,3,I+1II(4.15)x(k+1)=1Bx(k+1)=1Bb+x(k)0(k),(4.16)x(k+1)=0-s-I(E)xAs-tl(E)x44(4.17)=1.22x1011。其中s=e-r=0.4993,I(E)=I(E)=1,/(E)=e-0.421-,A=1.109=1.22x1011。123（2）捕撈量單位時間第i齡魚的捕撈量（條數(shù)）為yC）=q（E）xC）,ktk+1,iii第k年全年（8個月）第i齡魚的捕撈量（條數(shù)）為()(1-sH(Ex(k)iY(k)=Jty(t)dt=f()(1-sH(Ex(k)ii0i0ii于是，第k年總捕撈量

26、（質(zhì)量）為W（k）=wY（k）+wY（k）3344（3）可持續(xù)性捕撈模型可持續(xù)捕撈，即意味著由于自然死亡和捕撈使魚群減少，而通過產(chǎn)卵繁殖補(bǔ)充，使得魚群能夠在每年初開始捕撈時保持平衡不變，這樣的捕撈策略就可以年復(fù)一年地一直持續(xù)下去。因此，可持續(xù)捕撈的魚群數(shù)應(yīng)是（4.15）、（4.16）、（4.17）式的平衡解，即模型不依賴于時間t的解x*（i=0,1,2,3,4）。求解（4.15）、（4.16）、（4.17）得ix*i+1=si(E)x*，i=1,2,3x*i+1iix*2xx*2x*0將（4.18）式代入 HYPERLINK l bookmark179 o Current Document =

27、six*,x*=sx*=s2x*,x*=si(E)x*=s3i(E)x*132143331bx*=0,b+bx*0=0.5As-ti(E)x*+Ast-i(E)x*33444.20)式得代入（4.19）式有求解可得8x代入（4.19）式有求解可得8x*=0.5+si(E)As31(E)x0*31x*18b0.5+si(E)As31(E)x*=43，8b+0.5+si(E)As3i(E)x*431x*x*1(1)=巾-B(E)代入（4.18）式得到x*x*2(1)=sb1-Be)J，x*3x*3x*=s3l（Eb431、B（E）其中B（E）=0.5+si（E）As31（E）.43當(dāng)B（E）1時，x*0即意味著捕撈過度，致使魚群滅絕。當(dāng)B（E）=1時,0條數(shù)）為E=31.4稱之為過度捕撈力量，因此，可以在EE的范圍內(nèi)找最優(yōu)捕撈策略。0條數(shù)）為在可持續(xù)性捕撈的條件下，第i齡魚的年捕撈量qiY=qiY=i（E（1s-l（E）；+q（E）ix*ii=3,4。整個魚群的年捕撈量（重量）為Y（E）=wy+wy33440.42