差分方程方法

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ù)列xn ，把數(shù)列中的前 n 1項 xi i 0,1,2

2、, n 關(guān)聯(lián)起來所得到的方程。41 常系數(shù)線性差分方程4.1.1 常系數(shù)線性齊次差分方程常系數(shù)線性齊次差分方程的一般形式為xna1xn 1a2 xn 2akxn k 0(4.1)其中 k 為差分方程的階數(shù)， aii 1,2, ,k為差分方程的系數(shù)，且 ak 0 kn 。對應(yīng)的代數(shù)方程kk 1k2ka1 k 1a2 k 2ak0(4.2 )稱為差分方程的（ 4.1）的 特征方程， 其特征方程的根稱為 特征根。常系數(shù)線性齊次差分方程的解主要是由相應(yīng)的特征根的不同情況有不同的形式。下面分 別就特征根為單根、重根和復(fù)根的情況給出差分方程解的形式。1. 特征根為單根設(shè)差分方程（ 4.1）有 k 個單特征

3、根 1, 2, 3, , k ，則差分方程（ 4.1 ）的通解為xnc1 1c2 2ck kn ，其中 c1,c2,ck 為任意常數(shù)，且當(dāng)給定初始條件xii0i 1,2, ,k （4.3）時，可以唯一確定一個特解。2. 特征根為重根設(shè)差分方程（4.1 ）有|個相異的特征根1, 2, 3, I 1 l k重數(shù)分別為lm1,m2 , ,ml 且mi k 則差分方程（ 4.1 ）的通解為i1miXni 1 nC1i n11m2i 1 nQi n 2i 1mli5n同樣的，由給定的初始條件3. 特征根為復(fù)根4.3 ）可以唯一確定一個特解。設(shè)差分方程（4.1 ）的特征根為一對共軛復(fù)根和相異的k 2個單根

4、k ,則差分方程的通解為XnC1n cosnn C2 sin nC4其中arcta n.同樣由給定的初始條件（4.3 ）可以惟一確定一個特解?；蚬曹棌?fù)根和重根的情況,另外，對于有多個共軛復(fù)根和相異實根， 差分方程解的形式。4. 1. 2 常系數(shù)線性非齊次差分方程 常系數(shù)線性非齊次差分方程的一般形式為都可以類似地給出Xna1Xn 1a2Xn 2akXn k f n(4.4)其中k為差分方程的階數(shù),ai i1,2,k為差分方程的系數(shù),ak0 k n f(n)為已知函數(shù)。在差分方程（4.4 ）中，令f (n)0,所得方程Xn a1Xn 1a?XnakXn k 0(4.5)稱為非齊次差分方程（4.4

5、）對應(yīng)的齊次差分方程，即與差分方程（ 求解非齊次差分方程通解的一般方法為4.1 ）的形式相同。首先求對應(yīng)的齊次差分方程（4.5）的通解xn，然后求非齊次差分方程（4.4 ）的一個特解xn0，則XnXn0Xn為非齊次差分方程（4.4 ）的特解。關(guān)于求Xn的方法同求差分方程（4.1 ）的方法相同。對于求非齊次方程（4.4 ）的特解Xn0的方法，可以用觀察法確定，也可以根據(jù) f（n）的特性用待定系數(shù)法確定，具體方法可參照常系數(shù)線性非齊次微分方程求特解的方法。4.2差分方程的平衡點及其穩(wěn)定性而只需一般來說，差分方程的求解是困難，實際中往往不需要求出差分方程的一般解, 要研究它的平衡點及其穩(wěn)定性即可。4

6、.2.1一階線性常系數(shù)差分方程一階線性常系數(shù)差分方程的一般形式為Xk 1 aXk b, k 0,1,2,，其中a, b為常數(shù)，它的平衡點由代數(shù)方程x ax b求解得到，不妨記為 x*.如果lim xkx*，則稱平衡點x*是穩(wěn)定的，否則是不穩(wěn)定的。k k0的平衡點為了便于研究平衡點 x的穩(wěn)定性問題，一般將其轉(zhuǎn)化為求方程 xk 1 axkx0的穩(wěn)定性問題。事實上，由Xk 1 axk0可以解得kxka X。,*-r.曰 于疋x0是穩(wěn)定的平衡點的充要條件是：Ja 14.2.2一階線性常系數(shù)差分方程組一階線性常系數(shù)齊次差分方程組的一般形式為x k 1 Axk B,k 0,1,2,其中x k為n維向量，A

7、為n n階常數(shù)矩陣。,n).它的平衡點x0是穩(wěn)定的充要條件是A的所有特征根都有i 1 (i=1,2,對于一階線性常系數(shù)非齊次差分方程組x(k 1) Ax(k) B， k 0,1,2,的情況同樣給出4.2.3二階線性常系數(shù)差分方程二階線性常系數(shù)齊次差分方程的一般形式為xk 2 alxk 1a2xk0,k0,1,2,其中ai,a2為常數(shù)，其平衡點 x* 0是穩(wěn)定的充要條件是特征方程2 aia? 0 ,的根1, 2滿足11 ,21。對于一般的Xk 2 a1Xk 1 a?Xk 0的平衡點的穩(wěn)定性問題同樣給出。類似地，也可直接推廣到n階線性差分方程的情況。4.2.4一階非線性差分方程一階非線性方程的一般

8、形式為Xk 1f Xk ， k 0,1,2,其中f為已知函數(shù)，其平衡點定義為方程X f X的解X*。事實上，將f Xk在X*處作一階的泰勒展開有 * * *Xk 1f X Xk X f X ，則X*也是一階線性方程 Xk 1 f X* Xk X* f X*的平衡點，故此，平穩(wěn)衡點X*穩(wěn)定的充要條件是I f(X*)1。4.3連續(xù)模型的差分方法4.3.1微分的差分方程已知 f(x)在點 Xk處的函數(shù)值f(Xk)k 0,1,., n 1)， 且ax0x1. Xn 1 b，試求函數(shù)的導(dǎo)數(shù)值f Xk , (k 1,2,.n)。根據(jù)導(dǎo)數(shù)的定義，用差商代替微商，則有下面的差分公式。向前差：f (Xk)f(X

9、k1) f(Xk)(k 1,2,., n),f(Xk)f(Xk)f(Xk 1)XkXk 1(k1,2,., n),中心差：f (Xk)f (Xk 1)f (Xk 1)(k1,2,., n),Xk1Xk 14.3.2定積分的差分方法向后差:已知函數(shù)f(x)在點Xk處的函數(shù)值f(Xk) , (k1,2,., n),且在a,b上可積，試求函b數(shù)在a,b上的積分值 玄f(x)dx。f (x)dxnAk f(xQ根據(jù)定積分的定義，則有一般的求積公式其中Ak為求積系數(shù)，它與Xk的選取方法有關(guān)。取不同的求積系數(shù)，可以得不同的求積公式。對于等距節(jié)點xk akh (k0,1,., n)，其中步長h為很小的數(shù)，則

10、有如nki(2)復(fù)化梯形矩陣;ba f(x)dxh n 1-f(Xk)2 k 0f(Xk 1)hr z、f (a)21f(Xk) f (b)1(3)復(fù)化辛普森(Simpson)矩陣公式;b f (x)dx h f (xk) 4f(xa6k 0k1)2f(Xk1)其中x 1k 2(4)n 1f(a)4 f(xk 1)k 0 k 2f(Xk) f(b)1尹k Xk1)為子區(qū)間Xk,Xk1(k0,1,., n 1)的中點。復(fù)化柯特斯(Cotes)公式;下的求積公式。(1 )復(fù)化矩陣公式;bf(x)dxan 112 f(Xk) 7f(b)k 1baf(x)dx 907f(a) 32k0f(xq)12k

11、0f(xq)32k0f(y其中x 1 , x 1 ,x 3為子空間Xk,Xki(k0,1,.,n 1)中的四等分點。kkk .2444.3.3常微分方程的差分方法1. 一階常微分方程的差分方法 設(shè)一階常微分方程的定解問題為yf(x,y),y(xo) yo,其中函數(shù)f (x, y)關(guān)于y滿足李普希茲條件，即保證問題解的存在唯一性。現(xiàn)在的問題是求方程在一系列節(jié)點x1x2.xn.處的近似數(shù)值解y1, y2,., yn,.不妨假設(shè)步長為hxn 1 xn為常數(shù)。在此，我們根據(jù)微分的差分方法，即用差商來近似代替微商，再利用“步進(jìn)式”方法，可以給出求解問題(4-6)的差分方法。(1) 單步歐拉(Euler)

12、公式用差商y(xn 1)近似代替y(Xn)f(Xn,y(Xn)中的導(dǎo)數(shù)，則可以得差分公式hyn 1yn hf(Xn,yn), n 0,1,2,.2其精度為O(h )階的。(2) 兩步歐拉公式用差商 也山 山山 近似代替y(Xn)f(Xn,y(Xn)中的導(dǎo)數(shù)，則可得差分公式2hyn 1yn 1 2hf (xn, yn), n 0,1,2,.兩步法需要用到前兩步的方信息，一般不能自行起步，需先用單步方法求出y1，其精度是O(h2)階的。(3) 梯形公式對于方程yf (x, y)的兩邊在Xn,Xn 1上求積分得Xn 1y(Xn 1)y(Xn)y(Xn) x f (X,y(x)dX.Xn利用積分的差分

13、方法中梯形公式求解積分xn 1hx f(x,y(x)dx -f (Xn,y(Xn)f(Xn 1,y(Xn 1)n2hy(Xni)y(Xn) -f(Xn, y(Xn) f(X. 1, y(X. J)離散化即可得到微分方程的梯形差分公式hyn 1yN 2f(Xn，yN) f (Xn 1,yN 1), N 0,1,2,.這是一個隱式格式，計算量大，一般不單獨使用。其鏡的也是0(h3)階的。(4) 改進(jìn)的歐拉公式由于單步歐拉公式色精度低，但計算量??；矩形公式精度高，但是計算量大，為此，我們綜合運(yùn)用這兩種方法舅老爺?shù)玫礁倪M(jìn)的歐拉公式，其精度為Oh3階的。預(yù)報：Yn 1Yn hf Xn, yn , n 0

14、,1,校正：hyn 1 yn 2 f Xn,ynf Xn n 1 ,n 0,1,或?qū)懗善骄问?；Ypyn hf Xn,yn ,Yc yn hf Xn 1, Yp ,1Yn 1-(Yp Yc),n 0,1,2,.(5) 龍格一庫塔法龍格一庫塔方法的基本思想：對于微分方程的定解問題(4.6)，考慮差商Y(Xn Y(Xn)，根據(jù)阿格朗日微分中值h定理可得Y(Xn 1)Y XnhyY Xn hf , Y ， XnXn 1，記丫*f , Y ，稱為 Xn, Xn 1上的平均變化率，則 Y(Xn 1) Y XnhY *?，F(xiàn)在的問題只要找到尋找一種好的計算Y*的近似方法。如果取Y* f Xn, Yn丫1，

15、則就是歐拉公式。* 1如果取Yf人，ynf Xn 1, yn 1Y2，則相應(yīng)的就是改進(jìn)的歐拉公式?，F(xiàn)在，我們?nèi)個點Xnih, ynihxn,xn1i 1,2, , m，用f在這m個點的函數(shù)值的加權(quán)平均作為 Y*的近似值，即if Xnih,ynihi 1其中i為權(quán)系數(shù)。則有mYn 1 Yn hiif Xnih,Yn1ih（4.7）其中i， i，i為待定系數(shù)。實際上，適當(dāng)選擇i， i， i，使得公式有更高的精度,這是龍格-庫塔方法的思想。二階龍格-庫塔公式：在 Xn, Xn 1內(nèi)取中點X 1 Xn1h，則可取10，2 1，1代人（4.7）n .222式得到二階龍格-庫塔公式，其精度為0 h3 階

16、。Yn 1Yn hY2,f Xn, Yn ,Y2hh、/c fXn-，Yn二丫，n0,1,2,22三階龍格-庫塔公式：在Xn, Xn 1內(nèi)任取二點Xn pXn ph，Xn q x. qh，0 p q 1 類似的方法可得到三階的龍格-庫塔公式Y(jié)n 1Yn - Yl 4丫2 丫3 ,6丫 f Xn，Yn ,hh篦f Xn Yn 2，0 O1，丫3 f Xn h，Yn h Yi 2丫2其精度是Oh4階的常用的是三階的情況。四階龍格-庫塔公式：類似的方法可以得到四階龍格-庫塔公式，其精度是0 h5階的yn 1 yn 6 丫1 2丫2 2丫3 JY f Xn,yn ,丫2f Xn2，yn1，n，1,2f

17、 Xnh2，yn2丫2丫4 f(Xn1n hYJ2. 一階常微方程組的差分方法相應(yīng)的差分方法即可用于由多個方程組的將前面的單個方程中的變量和函數(shù)視為向量, 一階方程組的情形。對于二個方程的方程組校正:四階龍格-庫塔公式其中ynZnYn 1ynZnynZnf Xn, yn,Znh2 g Xn,yn,Zn-丫1 2丫2 2丫3-Z1 2Z2 2Z36 1 2 3f xn 1,n1 yn 1,Zn 1g Xn 1,yn 1,Z n 1丫4 ,n 0,1,2,Z4y1Zf X.y,z , g x, y, z ,Y XoZ XoYo,Zo,yn, Zn表示函數(shù)在節(jié)點Xnx0 nh,n1,2,上的近似解,

18、預(yù)報：yn 1Ynhf Xn , Yn , ZnZ n 1Znhf Xn , Yn , Zn(4.8)則有改進(jìn)的歐拉公式:設(shè)以丫1fXn, yn,Zn :i乙gXn,yn,Zn7hhh丫2fxn,ynY1 , Zn乙222hh、，h “Z 2gXn,yr丫1 , ZnZ1222hh、/h “丫3fXn,y丫2, ZnZ2222hh、，h “Z3gXn尹丫2 ,Zn22Z2丫4fXn1, y nhY3,Zn hZ3Z4gXn1 , y nhY3,Zn hZ3其他的公式也都可以類似得到， 即相當(dāng)于同時求解多個一階方程， 從方法上沒有本質(zhì)的差別。3. 高階常微分方程的差分方法對于某些高階方法的定解問

19、題，原則上可以轉(zhuǎn)化為一階方程組來求解。譬如，對于如下的二階微分方程的定解問題y f x, y, yIIy（xo）yo,y xy若令zy，則可以化為一階方程組的定解問題z f x, y, zy乙（4.9）y（x。）y,y xy。實際上,（4.9）式可以視為（4.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個月，平均每條

20、4齡魚的產(chǎn)卵量為52齡和1齡魚不產(chǎn)卵。卵孵化并成活為1.109* 10 （個），3齡魚的產(chǎn)卵量為這個數(shù)的一半,1齡魚，成活率（1齡魚條數(shù)與產(chǎn)卵量n之比）為111.22 1011.22 10 n漁業(yè)部門規(guī)定，每年只允許在產(chǎn)卵孵化期前的8個月內(nèi)進(jìn)行捕撈作業(yè)。 如果每年投入的捕撈能力固定不變，即固定努力量捕撈，這時單位時間捕撈量將與各年齡組魚群條數(shù)成正比， 比例系數(shù)稱為捕撈強(qiáng)度系數(shù)。通常使用13mm網(wǎng)眼的拉網(wǎng)，這種網(wǎng)只能捕撈3、4齡魚，其兩個捕撈系數(shù)之比為 0.42 : 1.要解決的問題是：建立數(shù)學(xué)模型，分析如何實現(xiàn)可持續(xù)性捕撈（即每年開始捕撈時漁場中各年齡組魚群 條數(shù)不變），并且在此前提下得到最高

21、的年收獲量（總質(zhì)量）。432模型的假設(shè)與符號說明1. 模型的假設(shè)（1）只考慮魚的繁殖和捕撈的變化，不考慮魚群遷入與遷出；（2）各齡魚在一年的任何時間都會發(fā)生自然死亡；（3）所有魚都在每年最后四個月內(nèi)（后 1/3年）完成產(chǎn)卵孵化的過程，成活的幼魚在 下一年初成為1齡魚；（4）產(chǎn)卵發(fā)生于后四個月之初，產(chǎn)卵魚的自然死亡發(fā)生與產(chǎn)卵之后；（5） 相鄰兩個年齡組的魚群在相鄰兩年之間變化是連續(xù)的，即第k年底i齡魚的條數(shù) 等于第k 1年初i 1齡魚的條數(shù)；（6）4齡以上的魚全部死亡；（7）采用固定努力量捕撈意味著捕撈的速率正比于捕撈時各齡魚群的條數(shù)，比例系數(shù) 為捕撈強(qiáng)度系數(shù)。2. 符號的說明用Xj（t）表示t

22、時刻（年）i齡魚的條數(shù)；r表示魚的平均自然死亡率，即 r 0.8 ；人A表示齡i魚的產(chǎn)卵數(shù)，即（ff2,f3, f4） （0,0, , A），A 1.109 105 ； Wi表示齡i魚群的捕撈強(qiáng)度系數(shù)，即（w1, w2 ,w3, w4）（5.07,11.55,17.86,22.99） ； qi表示i齡魚群的捕撈2強(qiáng)度系數(shù)，即（0,0,0.42E,E）， E為捕撈努力量；t -表示產(chǎn)卵開始的3Y月份；Yi表示i齡魚的捕撈量；Ci表示i齡魚的捕撈率，即 Ci L。Xi4.4.3 模型的建立與求解dxi tdt1. 無捕撈時魚群的自然增長模型 由假設(shè)（1）和（2 ）得rxi t ， i 1,2,3,

23、4 ； k t k 1， k 0,1,2,又由假設(shè)（3）和（4）得x-! k 11.22* 1011,X0 k1.22* 1011 x0 k tx k tX3Ax4 k t由假設(shè)（5）和（6 ）得Xi 0x1 ,xi 1 k 1Xi k 1 lim xi (t)ii 1it k 1 i(i 1,2,3;k0,1,)2。固定努力量捕撈魚群的增長和捕撈模型由假設(shè)知，捕撈期為k t k t，則有dxi trxi t qi Exi t , k t k t，(4.10)dtdxi trxi t ， kt t k 1，(4.11)dtXi 0X1，Xi 1 k1xi k 1 ， i 1,2,3，(4.12

24、)11x1 k 1(4.13)1.22*10 , Xo k111.22*10x0 k t則有Ax0 k tx3 k tAx4 k t , k 0,1,2,2（1）魚群的增長規(guī)律求解方程（4.10 ）和（4.11），并利用連續(xù)條件（4.12 ）式可得Xi 1 ksli E1,2,3,(4.15)其中sxMk1)b x0 kASt|3EX。X3Astl4 E(4.16)X4 k t ,(4.17)0.4993, |1 EI2 E1,I30.42 tEe , A1.109 105, b 1.22 1011。捕撈量單位時間第i齡魚的捕撈量（條數(shù)）為yi tqi E Xi t , k t k t ,第k

25、年全年（8個月）第i齡魚的捕撈量（條數(shù)）為Yi kt0yit dt0qiE Xi t dt蟲三（1r qi Estli E ）Xi k于是，第k年總捕撈量（質(zhì)量）為W（k） WsY3（k） W4Y4（k）（3）可持續(xù)性捕撈模型可持續(xù)捕撈，即意味著由于自然死亡和捕撈使魚群減少，而通過產(chǎn)卵繁殖補(bǔ)充，使得魚群能夠在每年初開始捕撈時保持平衡不變，這樣的捕撈策略就可以年復(fù)一年地一直持續(xù)下 去。因此，可持續(xù)捕撈的魚群數(shù)應(yīng)是（4.15 ）、（4.16 ）、（4.17 ）式的平衡解，即模型不依賴于時間 t 的解 X*（i 0,123,4）。求解（4.15 ）、（ 4.16 ）、（ 4.17 ）得X* 1sli

26、 Ex*， i 1,2,3，即* *2*3*X2 slx1, X3SX2s X1, X4 sb（E）X3 s l3 E X1* bx0X1*,b bx0* + * + *x00.5As l3（E）x3 As l4（E）x4將（4.18 ）式代入（4.20 ）式得8* *X0 0.5 sl4（E）As3l3（E）X1代入（4.19 ）式有8x*b0.5 sl4（E）As3l3（E）x；X18b 0.5 sl4（E）As?l3（E）x；求解可得* 1X1b 1 -B E代入（4.18 ）式得到X2sb 1*X3sb 1*3X4s 13 E b 1其中B E80.5 sl4（E）As3l3（E）.當(dāng)B（E） 1時,xi0即意味著捕撈過度，致使魚群滅絕。當(dāng)B（E） 1時,Eo 31.4稱之為過度捕撈力量，因此，可以在EEo的范圍內(nèi)找最優(yōu)捕撈策略。在可持續(xù)性捕撈的條件下，第i齡魚的年捕撈量（條數(shù)）為t*qi E 1 s li E XiY, i 3,4 。r qi E整個魚群的年捕撈量（重量）為Y EW3 y3W4 y40.42 1 stli Ew3r 0