脈沖微分方程的一個修正block-by-block數值格式?

曹俊英, 馬群長, 王自強

(1-貴州民族大學理學院，貴陽 550025;2-北京大學數學科學學院，北京 100871)

1 引言

隨著科學技術的不斷發展，許多反映客觀現實的物理模型都具有這樣的特征：在發展的某些階段，由于受到外部的作用或系統內部自身的原因，使得系統瞬間改變原有的狀態，這種瞬間改變原有狀態的現象通常稱為脈沖現象．脈沖現象在現代科技各領域的實際問題中是廣泛存在的，其數學模型往往可歸結為脈沖微分系統．脈沖微分系統能夠更深刻、更精確地反映事物的變化規律．隨著科學技術的突飛猛進，人們越來越認識到脈沖微分系統的重要性以及在實踐中的應用價值．例如，具有脈沖影響系統應用在物理、生物科技、工業、機器人、輻射能、藥物動力學、種群動態、生態、最優控制、微生物的研究、繁殖、經濟學、生產理論等方面．正是由于脈沖微分系統能夠更好地反映瞬間突變事物的變化規律．因此，它已經吸引了國內外眾多的專家和學者對脈沖微分方程的理論進行了研究[1-8]．

由于脈沖微分方程起源較早，研究脈沖微分方程的專家和學者已經從各個方面對脈沖微分方程進行了研究，從而使得脈沖微分方程理論取得了豐富的成果．但是對于脈沖微分方程數值方法的研究還剛剛起步，關于脈沖微分方程數值方法的成果還較少．2000年，Ranbelovi′c等人[9]基于著名的Runge-Kutta數值方法提出了一種計算脈沖微分方程的新算法．2005年，Hosseini[10]利用修正的譜方法求解具有脈沖項的常微分方程．2007年，Liu等人[11]研究了一個線性脈沖微分方程，并指出若將Runge-Kutta方法應用于該系統，那么該方法是數值穩定的．2008年，Ran等人[12]主要研究了線性脈沖微分系統的Runge-Kutta方法的數值解的穩定性．2008年，Liang和Liu[13]主要針對一類非線性的脈沖微分方程進行研究，在f,Ik分別滿足單邊Lipschitz條件和經典Lipschitz條件下得到了其解析解及數值解的穩定性．2010年，Ding和Wu[14]主要將定步長的Euler方法應用于帶延遲的線性脈沖微分方程，并且證明了該方法的收斂階是1階．2014年，Zhang和Liang[15]研究了脈沖微分方程的配置點法．2014年，Liu等人[16]研究了脈沖微分方程的迭代學習控制．

經典的b lock-by-block方法是求解積分方程的一個很有效的數值方法，2013年，Cao和Xu[17]基于經典的block-by-block方法，提出了一種求解了分數階微分方程的修正的block-by-block算法．2014年，王自強和曹俊英[18]利用修正的b lock-by-b lock方法，對非線性二維Volterra積分方程構造了一個高階數值格式．這里，我們將利用修正的block-by-block來求解脈沖微分方程，得到了一個新的高階數值格式，該格式除了每段的前兩層外，其余的未知量不需要耦合求解，并且對該格式進行收斂性分析，最后通過數值算例驗證該格式具有很好的收斂效果．

2 數值格式的構造

我們考慮脈沖微分方程如下

其中f:R+×?→R,Ik:?→R,??R是開區間，τk<τk+1，通常情況下?x(t)=x(t+0)?x(t),x(t+0)代表x在t點的右極限．

當f,Ik分別滿足單邊Lipschitz條件和經典Lipschitz條件下，問題(1)可以等價轉化為如下的積分方程

為了敘述方便，我們不妨假設脈沖點的個數是有限的，且相鄰的兩個脈沖點之間的距離是相等的，即τn+1?τn= λ,n=0,1,2,···,M；λ為常數，為了構造高階數值格式，將區間[τn,τn+1],n=0,1,2,···,m進行2N 等分，設= τn+jh,n=0,1,2,···,M,j=0,1,···,2N，其中h=則方程(2)在點上的數值解記

因此，由(1)我們可以得到關于脈沖項的式子如下

根據文獻[17]的離散思想，數值格式的具體構造思路如下：為了能夠計算最初的兩步解，我們需要先確定x(t)在的值．利用Sim pson求積公式，x(t)在上的逼近式為

其中的近似解．利用二次拉格朗日插值來逼近半點，則有

將(5)代入(4)式，有

相似地，計算x(t)在上的值，我們在上來計算則有

因此

從(6)和(7)我們發現，計算i=1,2需要知道f在,i=1,2上的值．所以，解x在,i=1,2的值需要同時求解(6)和(7)．

現在，我們構造下一步格式．假設,j=0,1,···,2m已經求出，下面給出和的逼近，按照相同的思路，有

因此，我們得到如下

可得如下

3 收斂性分析

設f(t,x)和Ik(x)關于變量x均滿足Lipschitz條件，即存在常數L1和L2，使得

定理 數值格式(10)是收斂的，而且收斂階數為4階．

這里的表示二次拉格朗日插值余項，表示Sim pson求積公式的誤差，且有

同理，有

聯立(13)和(14)，并且根據=0，得

同樣的方法，可得

當10時，則

利用G ronwall不等式[19]，有

同理，有

同理，有

聯立 (20)–(22)，有

同樣的方法，我們得

當1?0時，則

利用Gronwall不等式[19]，有

同理，有

因此，定理得到證明．

4 數值算例

本節，我們要進行一系列的數值試驗，來驗證理論結果的正確性．準確地說，我們的主要目的是驗證數值解關于步長的收斂行為．

例1 考慮如下的脈沖微分方程

它的精確解為

其中t和[t]分別代表t的小數部分和最大整數部分．

例2 考慮如下的脈沖微分方程

此方程的精確解為

其中t和[t]分別代表t的小數部分和最大整數部分．

在圖1和圖2中，取T=3和h=0.1時，利用數值格式(10)所得到的曲線圖形，圖1和圖2分別為例1和例2的數值解與精確解隨著t的變化情況．由圖1和圖2可以看出，我們所構造的數值格式(10)具有很好的收斂性．

圖1:T=3和h=0.1時，例1的數值解與精確解的比較

圖2:T=3和h=0.1時，例2的 數值解與精確解的比較

為了刻畫數值解的精度，考慮例1和例2的精確解與利用數值格式(10)所求得的數值解之間的最大誤差，即為了數值模擬方便，式(28)和式(29)中取T=3．從表1和表2可以看出，最大誤差隨步長h的減小而減小，且收斂階接近4階，這與理論預測結果相吻合．

表1: 例1最大誤差隨步長h的變化情況

表2: 例2最大誤差隨步長h的變化情況

根據文獻[12]可知求解脈沖微分方程(1)的θ-方法的數值格式為

取T=3，利用式(30)分別對例1和例2進行數值求解．表3和表4分別列出了θ和h取不同值時最大誤差的變化情況，并列出了相應的收斂階數．從表3和表4中可知：當θ=0.2時，數值格式(30)的收斂階接近1，當θ=0.5時，收斂階數接近2．

表3: 式(30)中θ=0.2和θ=0.5求解例1的最大誤差隨步長h的變化情況

表4: 式(30)中θ=0.2和θ=0.5求解例2的最大誤差隨步長h的變化情況

從上述表中可以看出，無論是最大誤差或收斂階數，我們所構造出來的數值格式(10)優于利用θ-方法(30)來求解此類脈沖微分方程．

參考文獻：

[1]Ballinger G,Liu X Z.Existence and uniqueness resu lts for Impulsive delay d iff erential equations[J].Dynam ics of Continuous D iscrete and Impulsive System s,1999,5(1-2):579-591

[2]Berezansky L,B raverm an E.Exponential boundedness of solutions for Impulsive delay d iff erential equations[J].App lied mathem atics Letters,1996,9(6):91-95

[3]Liu B,Liu X Z,Teo K,et al.Razum ikhin-type theorem on exponential stability of Impulsive delay system s[J].IM A Jou rnal of App lied mathem atics,2006,71(1):47-61

[4]Shen J H,Yan J R.Razum ikhin-type stability theorem s for Impulsive function d iff erential equations[J].Non linear Analysis,Theory,methods and App lications,1998,33(5):519-531

[5]Wang Q,Liu X Z.Exponentialstability for im pulsive delay diff erentialequations by Razum ikhin method[J].Jou rnal of mathem atical Analysis and App lications,2005,309(2):462-473

[6]Li Y K,X ing W Y,Lu L H.Existence and global exponential stability of period ic solution of a class of neural networks with im pu lses[J].Chaos,Solitons and Fractals,2006,27(2):437-445

[7]Bainov D D,Sim eonov pS.System s with Impulsive Effect:Stability,Theory,and App lications[M].Chichester:E llis Horwood Ltd,1989

[8]Akhm et mU.On the general problem of stability for im pulsive differential equations[J].Journal of mathem atical Analysis and App lications,2003,288(1):182-196

[9]Ranbelovi′c B M,Stefanovi′c L V,Dankovi′c B M.Num erical solu tion of Impulsive d iff erential equations[J].Facta Universitatis,2000,15(15):101-111

[10]Hosseini mM.Num erical solution of ord inary differential equations with im pu lse solu tion[J].App lied mathem atics and Com putation,2005,163(1):373-381

[11]Liu mZ,Liang H,Yang ZW.Stability of Runge-Ku ttam ethods in the num ericalsolution of linear Impulsive differential equations[J].App lied mathem atics and Com pu tation,2007,192(2):346-357

[12]Ran X J,Liu mZ,Zhu Q Y.Num ericalm ethod for Impulsive differential equation[J].mathem atical and Com puter modeling,2008,48(1-2):46-55

[13]Liang H,Liu mZ.Stability of Runge-K utta methods in the num erical solu tion of non linear Impulsive differential equations[J].Journal of Natural Science Heilongjiang University,2008,23(4):469-472

[14]Ding X H,W u K N.The Eu ler schem e and its convergence for Impulsive delay d iff erential equations[J].App lied mathem atics and Com pu tation,2010,216(5):1566-1570

[15]Zhang Z H,Liang H.Collocation methods for Impulsive differential equations[J].App lied mathem atics and Com putation,2014,228:336-348

[16]Liu S D,Wang J R,W eiW.A study on iterative learning control for Impulsive differential equations[J].Comm unications in Non linear Science and Num erical Simulation,2014,24(1-3):4-10

[17]Cao J Y,Xu C J.A high order schem a for the num erical solu tion of the fractional ord inary d iff erential equations[J].Journal ofComputational Physics,2013,238(1):154-168

[18]王自強，曹俊英.非線性二維Volterra積分方程的一個高階數值格式[J].高校應用數學學報，2014,29(4):397-411 Wang Z Q,Cao J Y.A high order schem e for the num erical solution of the nonlinear two-dim ensional Volterra integral equations[J].App lied mathem atics A Journal of Chinese University,2014,29(4):397-411

[19]Dixon J,mcK ee S.W eak ly singu lar discrete G ronwall inequalities[J].Jou rnal of App lied mathem atics and mechanics,1986,66(11):535-544

[20]王自強，曹俊英.分數階微分方程的b lock-by-b lock算法的最優階收斂性分析[J].工程數學學報，2015,32(4):533-545 Wang Z Q,Cao J Y.Op tim al convergence order analysis of a b lock-by-b lock algorithm for fractional differential equations[J].Chinese Journal of Engineering mathem atics,2015,32(4):533-545