非連續(xù)系統(tǒng)的Simulink仿真方法研究

1、第18卷第7期 系統(tǒng) 仿 真 學(xué) 報(bào)© V ol. 18 No. 72006年7月 Journal of System Simulation July, 2006非連續(xù)系統(tǒng)的Simulink 仿真方法研究向 博1,高丙團(tuán)1,張曉華1,胡廣大2(1.哈爾濱工業(yè)大學(xué)電氣工程學(xué)院,哈爾濱 150001;2.哈爾濱工業(yè)大學(xué)航天學(xué)院,哈爾濱 150001 摘 要:針對(duì)滑模變結(jié)構(gòu)控制、最優(yōu)控制等控制系統(tǒng)中出現(xiàn)的由于不連續(xù)性問(wèn)題而直接影響Simulink 仿真速度和仿真效果的情況,分析了該現(xiàn)象產(chǎn)生的原因,并提出了相應(yīng)的解決策略。結(jié)合吊車(chē)防擺變結(jié)構(gòu)控制系統(tǒng),采用提出的解決方案進(jìn)行了Simulink

2、仿真實(shí)驗(yàn)研究,仿真實(shí)驗(yàn)的結(jié)果證明所提出的策略對(duì)于解決這類(lèi)不連續(xù)性仿真問(wèn)題具有較好的效果。最后,給出了解決和避免這種問(wèn)題的幾點(diǎn)建議。 關(guān)鍵詞:非連續(xù)系統(tǒng);Simulink 仿真;滑模變結(jié)構(gòu)控制;吊車(chē)中圖分類(lèi)號(hào):TP391.9 文獻(xiàn)標(biāo)識(shí)碼:A 文章編號(hào):1004-731X (2006 07-1750-05Simulink Simulation Research of Uncontinuous SystemXIANG Bo, GAO Bing-tuan, ZHANG Xiao-hua, HU Guang-da(1.School of Electrical engineering, Harbin In

3、stitute of Technology, Harbin 150001, China;2.School of Astronautics, Harbin Institute of Technology, Harbin 150001, ChinaAbstract: The Simulink simulation speed and effect are decreased because of the uncontinuous functions in the uncontinuous systems, such as sliding mode variable structure contro

4、l systems, optimal control systems and so on. The reasons caused the problem were analyzed in detail and corresponding resolving methods were also developed. Based on the presented methods, the Simulink simulation on variable structure control of an overhead crane system was studied, and the simulat

5、ion results demonstrate the effectiveness of the prompted Simulink simulation schemes. Finally, some suggestions on the problem were summed up.Key words: uncontinuous system; Simulink simulation; sliding mode variable structure control; overhead crane引 言 Simulink 是MathWorks 公司于1990年推出的產(chǎn)品,是用于MATLAB 下

6、建立系統(tǒng)框圖和仿真的環(huán)境。該環(huán)境剛推出時(shí)的名字叫Simulab ,由于其名字很類(lèi)似于當(dāng)時(shí)的一個(gè)很著名的語(yǔ)言Simula 語(yǔ)言,所以次年更名為Simulink 。從名字上看,立即就能看出該程序有兩層含義1,首先,“Simu ”一詞表示它可以用于計(jì)算機(jī)仿真,而“Link ”一詞表明它能進(jìn)行系統(tǒng)連接,即把一系列模塊連接起來(lái),構(gòu)成復(fù)雜的系統(tǒng)模型。正是由于它的這兩大功能和特色,使得它成為仿真領(lǐng)域首選的計(jì)算機(jī)環(huán)境。Simulink 是一個(gè)用來(lái)對(duì)動(dòng)態(tài)系統(tǒng)進(jìn)行建模、仿真分析的軟件包,它支持連續(xù)、離散及兩者混合的線性和非線性系統(tǒng),也支持具有多種采樣速率的多速率系統(tǒng)。其利用數(shù)學(xué)模型來(lái)獲取系統(tǒng)的一些重要特性參數(shù),這

7、些數(shù)學(xué)模型通常是以時(shí)間為變量的常微分方程來(lái)描述,并用數(shù)值方法進(jìn)行計(jì)算機(jī)仿真求解。所以,Simulink 軟件包常被用于運(yùn)動(dòng)學(xué)仿真來(lái)反復(fù)求解 機(jī)構(gòu)運(yùn)動(dòng)約束方程2,通過(guò)積分獲得最終的速度(或加速度,收稿日期:2005-05-25 修回日期:2006-04-03作者簡(jiǎn)介:向博(1980-,男,四川廣元人,碩士生,研究方向?yàn)榛诨Ｗ兘Y(jié)構(gòu)控制的吊車(chē)防擺系統(tǒng);張曉華(1961-,男,黑龍江哈爾濱人,教授,博導(dǎo),哈爾濱工業(yè)大學(xué)電氣工程學(xué)院常務(wù)副院長(zhǎng),研究方向?yàn)楣艿罊C(jī)器人,MATLAB 仿真技術(shù)。從而確定機(jī)構(gòu)運(yùn)動(dòng)的位置(和速度。圖1為運(yùn)動(dòng)學(xué)仿真的流程圖,它描述了諸如MATLAB 仿真工具箱Simulink

8、仿真軟件包的仿真步驟和流程2。圖1 運(yùn)動(dòng)學(xué)仿真流程框圖利用Simulink 仿真軟件包來(lái)求解機(jī)構(gòu)運(yùn)動(dòng),比通常求解機(jī)構(gòu)運(yùn)動(dòng)的方法具有許多優(yōu)點(diǎn)。最顯著的優(yōu)點(diǎn)是位置求解問(wèn)題可以用軟件包隱式求解,使用者只需為仿真提供適當(dāng)?shù)某跏紬l件即可以求解機(jī)構(gòu)在任意時(shí)刻的位置問(wèn)題。但是,如果在仿真模型中一旦出現(xiàn)了非連續(xù)模塊,則可能出現(xiàn)系統(tǒng)仿真停滯的情形。本文針對(duì)這一現(xiàn)象,以“基于滑模變結(jié)構(gòu)控制的橋式吊車(chē)防擺系統(tǒng)”為例,詳細(xì)分析了該現(xiàn)象產(chǎn)生的原因,并提出了相應(yīng)的解決策略。1 問(wèn)題的提出本文以“基于滑模變結(jié)構(gòu)控制的吊車(chē)防擺系統(tǒng)”為例,來(lái)研究非連續(xù)系統(tǒng)的Simulink 仿真方法。2006年7月 向 博,等:非連續(xù)系統(tǒng)的S

9、imulink 仿真方法研究 July, 2006根據(jù)文獻(xiàn)3,吊車(chē)的物理模型如圖2所示:x 圖2 吊車(chē)的物理模型重物通過(guò)繩索與小車(chē)相連,小車(chē)在行走電機(jī)的水平驅(qū)動(dòng)力f 的作用下在水平軌道上運(yùn)動(dòng)。小車(chē)的質(zhì)量為M ,重物的質(zhì)量為m ,繩索的長(zhǎng)度為l ,小車(chē)與水平軌道的摩擦阻尼系數(shù)為D ,g 為重力加速度時(shí)間常數(shù)。利用拉格朗日方程,我們可以得到吊車(chē)的數(shù)學(xué)模型為: f Mx Dx mg x l g =+=+(1 1.1 滑模變結(jié)構(gòu)控制器的設(shè)計(jì)設(shè)1x =,2x = ,1l m =,29.8/g m s =,并取控制:a x =則根據(jù)式(1,可寫(xiě)出如下二階系統(tǒng):12219.8xx x x a =+ (2 取

10、滑模面函數(shù)為:112s c x x =+ 10c > (3采用等速趨近率控制策略:sgn(dss dt= 0> (4 設(shè)李亞甫諾夫函數(shù):212V s =(5 則可得:sgn(0Vs s s s s =< (6 由于在切換面鄰域內(nèi)函數(shù)式(5是正定義的,而式(6中2s 的導(dǎo)數(shù)是負(fù)半定義的。所以,在s =0附近V 是一個(gè)非增函數(shù),因此,如果滿(mǎn)足條件式(6,則定義函數(shù)式(5是系統(tǒng)的一個(gè)條件李亞普諾夫函數(shù)。系統(tǒng)本身也就穩(wěn)定于條件0s =。系統(tǒng)將沿滑模面在t 趨于無(wú)窮大時(shí)歸于平衡點(diǎn)。由(2、(3、(4式變換可得:1121129.8sgn(a x c x c x x =+(7 式(2、(7

11、即構(gòu)成了具有自動(dòng)防擺功能的二維滑模變結(jié)構(gòu)控制器。1.2 系統(tǒng)仿真針對(duì)前面所設(shè)計(jì)的滑模變結(jié)構(gòu)控制系統(tǒng)進(jìn)行仿真,搭建的仿真模型如圖3所示。在此模型中,我們?nèi)坎捎肧imulink 工具箱中的模塊搭建模型。圖3 滑模變結(jié)構(gòu)控制器仿真模型取11c =,1=,1(01x =,2(00x =,仿真時(shí)間T=10s , 采用Variable-step 的ODE45算法。仿真結(jié)果如圖4所示:圖4 系統(tǒng)時(shí)域仿真結(jié)果由上圖可以看出,當(dāng)系統(tǒng)仿到1s 左右時(shí),速度非常慢,停滯不前,出現(xiàn)了圖形未完全仿真完成的情形,這顯然是我們不希望看到的。2 問(wèn)題的分析為了更好的分析該問(wèn)題產(chǎn)生的原因以及提出相應(yīng)的解決策略,首先,我們認(rèn)為

12、有必要從Simulink 仿真運(yùn)行原理開(kāi)始進(jìn)行分析。2.1 Simulink 仿真運(yùn)行原理對(duì)Simulink 系統(tǒng)模型的仿真主要包含兩個(gè)階段4:初始化階段和運(yùn)行階段1.初始化階段在初始化階段主要要完成以下工作:(1 每個(gè)模塊的所有參數(shù)都傳遞給MATLAB 進(jìn)行求值,得到的數(shù)值作為實(shí)際的參數(shù)使用;(2 展開(kāi)模型的層次結(jié)構(gòu),每個(gè)子系統(tǒng)被他們所包含的模塊替代,帶有觸發(fā)和使能模塊的子系統(tǒng)被視為原子單元進(jìn)行處理;2006年7月系統(tǒng)仿真學(xué)報(bào) July, 2006(3 檢查信號(hào)的寬度和模塊的連接情況,提取狀態(tài)和輸入/輸出依賴(lài)關(guān)系方面的信息,確定模塊的更新順序;(4 確定狀態(tài)的初值和采樣時(shí)間;完成這些工作后就

13、可以進(jìn)行仿真了。2.運(yùn)行階段初始化之后,仿真進(jìn)入運(yùn)行階段。仿真是由求解器控制的,它計(jì)算模塊的輸出,更新離散狀態(tài),計(jì)算連續(xù)狀態(tài)。在采用變步長(zhǎng)求解器時(shí),求解器還會(huì)確定時(shí)間步長(zhǎng)。計(jì)算連續(xù)狀態(tài)包含下面兩個(gè)步驟:(1 首先,求解器為待更新的系統(tǒng)提供當(dāng)前狀態(tài)、時(shí)間和輸入值。反過(guò)來(lái),求解器需要狀態(tài)導(dǎo)數(shù)的值;(2 然后,求解器對(duì)狀態(tài)的導(dǎo)數(shù)進(jìn)行積分,計(jì)算新的狀態(tài)的值。狀態(tài)計(jì)算完成后,再進(jìn)行一次模塊的輸出更新,這時(shí),一些模塊可能會(huì)發(fā)出過(guò)零的警告,促使求解器探測(cè)出發(fā)生過(guò)零的準(zhǔn)確時(shí)間。Simulink的仿真過(guò)程是在Simulink求解器和系統(tǒng)相互作用之下完成的,系統(tǒng)和求解器在仿真過(guò)程中的對(duì)話作用如圖5所示。 圖5 系

14、統(tǒng)與求解器之間的對(duì)話在圖5中,求解器的作用是傳遞模塊的輸出,對(duì)狀態(tài)導(dǎo)數(shù)進(jìn)行積分,并確定采樣時(shí)間。系統(tǒng)的作用是計(jì)算模塊的輸出,狀態(tài)進(jìn)行更新,計(jì)算狀態(tài)的導(dǎo)數(shù),生成過(guò)零事件、從求解器傳遞給系統(tǒng)的信息包括時(shí)間、輸入和當(dāng)前狀態(tài);反過(guò)來(lái),系統(tǒng)為求解器提供模塊的輸出、狀態(tài)的更新和狀態(tài)的導(dǎo)數(shù)。2.2 過(guò)零檢測(cè)過(guò)零檢測(cè)5是通過(guò)在系統(tǒng)和求解器之間建立對(duì)話的方式工作,對(duì)話包含的一個(gè)內(nèi)容是事件通知,即系統(tǒng)告知求解器在前一時(shí)間步發(fā)生了個(gè)事件。事件由過(guò)零表示,過(guò)零在下列兩個(gè)條件下產(chǎn)生:(1 信號(hào)在上一個(gè)時(shí)間步改變了符號(hào)(含變?yōu)镺和離開(kāi)O;(2 模塊在上一個(gè)時(shí)間步改變了模式(例如積分器進(jìn)入了飽和區(qū)段。過(guò)零是一個(gè)重要的事件,

15、表征系統(tǒng)中的不連續(xù)性。例如,響應(yīng)中的跳變。如果仿真中不對(duì)過(guò)零進(jìn)行檢測(cè),可能會(huì)導(dǎo)致不準(zhǔn)確的仿真結(jié)果。當(dāng)采用變步長(zhǎng)求解器時(shí),Simulink 能夠檢測(cè)到過(guò)零(使用固定步長(zhǎng)的求解器,Simulink不檢測(cè)過(guò)零。當(dāng)一個(gè)模塊通知系統(tǒng)前一時(shí)間步發(fā)生了過(guò)零,變步長(zhǎng)求解器就會(huì)縮小步長(zhǎng),即便絕對(duì)誤差和相對(duì)誤差是可接受的?？s小步長(zhǎng)的目的是判定事件發(fā)生的準(zhǔn)確時(shí)間。當(dāng)然,這樣會(huì)降低仿真的速度,但這樣做對(duì)有些模塊來(lái)講是至關(guān)重要和必要的,因?yàn)檫@些模塊的輸出可能表示了一個(gè)物理值,它的零值有著重要的意義。事實(shí)上,只有少量的模塊能夠發(fā)出過(guò)零事件通知。每個(gè)模塊發(fā)出專(zhuān)屬于自己的事件通知,而且可能與不止一個(gè)類(lèi)型的事件發(fā)生關(guān)聯(lián)。能夠產(chǎn)

16、生過(guò)零通知的Simulink模塊如圖6所示:圖6 能夠產(chǎn)生過(guò)零通知的SIMULINK模塊各個(gè)模塊的過(guò)零的類(lèi)型是有差異的。例如,abs模塊在輸入改變符號(hào)時(shí)產(chǎn)生個(gè)事件,而saturation模塊能夠生成兩個(gè)不同的過(guò)零,一個(gè)用于下飽和,一個(gè)用于上飽和。signals&system庫(kù)中的HitCrossing模塊輸入穿過(guò)零點(diǎn)時(shí)產(chǎn)生一個(gè)過(guò)零,可以用來(lái)為不帶過(guò)零能力的模塊提供過(guò)零檢測(cè)的能力。一些過(guò)零只是用來(lái)通知求解器模式已經(jīng)發(fā)生了改變,另外一些則與信號(hào)相關(guān),用于觸發(fā)其他模塊。觸發(fā)包含三種類(lèi)型:(1 上升沿信號(hào)上升到或穿過(guò)零或者信號(hào)離開(kāi)零變正;(2 下降沿信號(hào)下降到或穿過(guò)零或者信號(hào)離開(kāi)零變負(fù);(3

17、雙邊沿上升或下降兩者之一發(fā)生。注意,Fcn模塊不支持過(guò)零,結(jié)果一些拐角點(diǎn)被漏掉了。2.3 算法的選擇Simulink提供了許多微分方程的解法6,如Runge-Kutta 法、Adams法、Gear法、Euler法、Linsim法等等。需要說(shuō)明的是,對(duì)于不同的仿真問(wèn)題,需要按照系統(tǒng)的特性,選擇不同的仿真方法。一般來(lái)說(shuō):Runge-Kutta方法(包含rk45、rk23變步長(zhǎng)方法適合于高度非線性或不連續(xù)的系統(tǒng),不適合于剛性系統(tǒng)(即有快變特性又有慢變特性的系統(tǒng);Adams方法適合于非線性小、時(shí)間常數(shù)變化小的系統(tǒng);Gear方法是專(zhuān)門(mén)用于剛性系統(tǒng)仿真的方法,但對(duì)非剛性系統(tǒng)較差;Euler方法比較差,盡量

18、避免使用;Linsim方法適合于近似線性的系統(tǒng),對(duì)線性剛性系統(tǒng)有很大的優(yōu)越性。一般來(lái)說(shuō),使用變步長(zhǎng)的自適應(yīng)算法是較好的選擇。這類(lèi)算法會(huì)依照給定的精確度在各積分區(qū)內(nèi)自適應(yīng)地尋找各自的最大步長(zhǎng)進(jìn)行積分,從而使得速率最高。Simulink的變步2006年7月向博,等:非連續(xù)系統(tǒng)的Simulink仿真方法研究July, 2006長(zhǎng)解法能夠把積分段分到足夠細(xì)以至到滿(mǎn)足精度要求的解。2.4 系統(tǒng)Simulink仿真停滯原因分析綜上所述,我們可以發(fā)現(xiàn):由于系統(tǒng)存在不連續(xù)模塊signum,所以,當(dāng)系統(tǒng)于1s左右到達(dá)滑模面(s=0時(shí),signum 模塊向系統(tǒng)發(fā)出過(guò)零通知。而當(dāng)采用變步長(zhǎng)求解器時(shí), Simulin

19、k能夠檢測(cè)到過(guò)零。所以,當(dāng)Simulink檢測(cè)到過(guò)零時(shí),便自動(dòng)縮小步長(zhǎng),可是,在下一仿真步里,系統(tǒng)繼續(xù)過(guò)零。如圖7和如圖8所示,因?yàn)榛Ｃ嬖?s處不能正常歸零,所以signum模塊就反復(fù)過(guò)零,同時(shí)一直向求解器發(fā)出過(guò)零通知。求解器便相應(yīng)的一直不停的縮小步長(zhǎng)。 圖7 滑膜面和signum模塊的時(shí)域響應(yīng) 圖8 1s左右時(shí)的滑模面和signum響應(yīng)的局部放大圖如圖9所示,系統(tǒng)大約經(jīng)過(guò)12個(gè)仿真步到達(dá)1s處時(shí),步長(zhǎng)急劇縮小至接近于零。這樣,由于仿真步長(zhǎng)太小,系統(tǒng)便在不連續(xù)處形成了過(guò)多的點(diǎn),超出了系統(tǒng)可用的內(nèi)存和資源,使得系統(tǒng)進(jìn)展緩慢,仿真停滯不前。 圖9 系統(tǒng)仿真步長(zhǎng)3 問(wèn)題的對(duì)策通過(guò)以上分析,我們總結(jié)

20、、概括得出了以下系統(tǒng)仿真速度慢的原因:(1 系統(tǒng)方程中存在不連續(xù)函數(shù)sgn(s;(2 Simulink仿真模型中存在能夠產(chǎn)生過(guò)零通知的signum模塊;(3 采用的變步長(zhǎng)求解器具有過(guò)零檢測(cè)并自動(dòng)調(diào)整步長(zhǎng)的功能。3.1 解決策略因此,基于以上原因,我們提出了以下4種解決策略:(1 取消Zero crossing detection功能;(2 采用fixed-step求解器;(3 采用不能夠產(chǎn)生過(guò)零通知的Fcn函數(shù)模塊;(4 柔化sgn(s函數(shù),使其連續(xù)化。以上所提出的4種解決策略中,采用單獨(dú)任何一種或幾種都可行。以下,我們嘗試采用第3種策略來(lái)觀察系統(tǒng)的仿真效果。即采用Fcn函數(shù)模塊來(lái)代替圖3中的

21、其他模塊所表示的數(shù)學(xué)關(guān)系,搭建的仿真模型如圖10所示: 圖10 采用Fcn函數(shù)模塊后的系統(tǒng)仿真模型對(duì)于函數(shù)fun,我們?nèi)?x,1x,a分別作為輸入u1, u2, u3。根據(jù)式(2可得函數(shù)fun的表達(dá)式為:-9.8*u2+u3對(duì)于函數(shù)fun1,我們?nèi)?x,1x,1c,分別作為輸入u1, u2, u3, u4。根據(jù)式(7可得函數(shù)fun1的表達(dá)式為: 9.8*u2-u3*u1-u4*sgn(u3*u2+u1取11c=,1=,系統(tǒng)仿真結(jié)果如圖11所示: 圖11 系統(tǒng)時(shí)域圖2006年7月 系 統(tǒng) 仿 真 學(xué) 報(bào) July, 2006由于Fcn 模塊不支持過(guò)零,所以,系統(tǒng)在不連續(xù)的情況下仍然能迅速完成仿真

22、。除此之外,我們還采用了fixed-step 求解器和變步長(zhǎng)下置Zero crossing detection 為off 的仿真方法,其仿真結(jié)果亦如圖11所示?？梢?jiàn),這前三種方法的仿真結(jié)果一樣。由上面的仿真結(jié)果可以看出:(1 取消系統(tǒng)的過(guò)零檢測(cè)功能之后,仿真速度快; (2 但是由于仍然存在不連續(xù)模塊signum ,所以加速度存在較大的抖振。3.2 sgn(s函數(shù)連續(xù)化 由于加速度存在劇烈的抖振,使得控制系統(tǒng)難于工程實(shí)現(xiàn),為了能消除系統(tǒng)因?yàn)椴贿B續(xù)性而存在的抖振現(xiàn)象,同時(shí)加快仿真速度,且不影響系統(tǒng)的仿真效果,我們?cè)O(shè)計(jì)了第四種控制策略-sgn(s函數(shù)連續(xù)化。連續(xù)化前后的函數(shù)sgn(s函數(shù)和sat(s

23、函數(shù)如圖12所示: 圖12 sgn(s和sat(s函數(shù)其中,sat(s函數(shù)的表達(dá)式為:1(1s sat s s = <(8 為大于零的正數(shù),且。當(dāng)取無(wú)窮小時(shí),sat(s 函數(shù)便非常逼近sgn(s函數(shù)。采用sat(s函數(shù)來(lái)代替sgn(s函數(shù)后的系統(tǒng)仿真模型如圖13所示: 圖13 柔化sgn(s函數(shù)后的系統(tǒng)仿真模型根據(jù)式(8可得函數(shù)fun 的表達(dá)式為:sin(u/0.1414。其表示,s 時(shí)sat(s的值。仿真結(jié)果為:(I 當(dāng)0=時(shí),連續(xù)函數(shù)sat(s即轉(zhuǎn)化為符號(hào)函數(shù)sgn(s,仿真時(shí)將出現(xiàn)如前面所述的問(wèn)題(仿真速度慢、抖振等。(II 當(dāng)0.1=時(shí),其它參數(shù)與圖3中相同,可得系統(tǒng)仿真結(jié)果如圖

24、14所示:圖14 采用連續(xù)函數(shù)(0.1=后的系統(tǒng)時(shí)域圖仿真結(jié)果說(shuō)明:采用連續(xù)函數(shù)sat(s代替不連續(xù)函數(shù)sgn(s后,系統(tǒng)仿真速度加快,不影響系統(tǒng)控制效果,且加速度不存在抖振現(xiàn)象,易于工程實(shí)現(xiàn),滿(mǎn)足我們所期望的結(jié)果。(III 當(dāng)0.0001=時(shí),系統(tǒng)的仿真結(jié)果如圖15所示:圖15 0.0001=時(shí)的系統(tǒng)仿真結(jié)果仿真結(jié)果說(shuō)明:當(dāng)0.0001=時(shí),加速度出現(xiàn)劇烈的抖振,性能接近理想的滑模變結(jié)構(gòu)控制。通過(guò)以上的仿真結(jié)果,我們可以總結(jié)如下:(1 采用sgn(s函數(shù)連續(xù)化的方法,系統(tǒng)仿真速度快; (2 越小,系統(tǒng)仿真結(jié)果越逼近原理想的滑模變結(jié)構(gòu)控制系統(tǒng),但抖振也越大,相應(yīng)的仿真速度亦越慢;反之,越大,仿

25、真結(jié)果失真也越大,但系統(tǒng)抖振越小,相應(yīng)的仿真速度亦越快。因此,我們建議:(1 如果希望既不影響系統(tǒng)的仿真速度,又不損失精度,0.0001為宜;第 18 卷第 7 期 2006 年 7 月 Vol. 18 No. 7 系 統(tǒng) 仿 真 學(xué) 報(bào) 表4 步長(zhǎng) h=0.1 求解例 2 在 x July, 2006 4 數(shù)值試驗(yàn) 例 1（參見(jiàn)5） 常數(shù)延遲微分方程初值問(wèn)題 = 10 處的整體誤差. x = 10 處整體誤差 1.6235864E-4 6.70327E-4 y '(x = 1000y(x + 997e y(x 1 +1000 997e 3x y(x =1+ e 3 3 0 x4 x0

26、 h=0.2 上面兩個(gè)例子驗(yàn)證了我們所構(gòu)造的方法的收斂階， 并且 驗(yàn)證了方法在求解剛性問(wèn)題時(shí)是有效的。 問(wèn)題的精確解是 y ( x = 1 + e 的整體誤差列在表 3 中。 表3 步長(zhǎng) h=0.1 h=0.01 3 x 。 用表 2 列出三階方法求解例 1，在積分區(qū)間0,4內(nèi)得到 參考文獻(xiàn): 1 H Podhaisky, B A Schmitt, R Weiner. Design, analysis and testing of some parallel two-step W-methods for stiff system J. Applied Numerical Mathematics

27、(S0168-9274, 2002, 42: 381-395. 2 3 H. Podhaisky, R. Weiner, B.A. Schmitt, Two-step W-methods for stiff ODE systems J. Vietnam J. Math. (S0866-7179, 2002, 30: 591-603. Bernhard A. Schmitt , Rüdiger Weiner, Parallel Two-Step W-Methods with 4 Peer Variables J. SIAM Journal on Numerical Analysis(S

28、0036-1429, 2004, 42: 265-282. K J int Hout, M N Spijker. Stability analysis of numerical methods for delay differential equations J. Numer. Math. (S0006-3835, 1991, 59: 807-814. 5 6 曹學(xué)年，劉德貴，李壽佛. 求解延遲微分方程的 ROSENBROCK 方法的漸近穩(wěn)定性J. 2002，14(3：290-292. E Hairer, S P NØrsett, G.Wanner. Solving Ordinary

29、 Differential Equations I. Nonstiff problems M. Berlin: Springer-Verlag , 1993. 求解例 1 在區(qū)間0,4上的整體誤差 區(qū)間0,4上的整體誤差 1.1762236E-05 1.5437829E-08 例 2（參見(jiàn)6） y1' ( x = y1 ( x y2 ( x 1 + y2 ( x 10 ' y2 ( x = y1 ( x y2 ( x 1 y2 ( x ' y3 ( x = y2 ( x y2 ( x 10 100 y3 ( x 當(dāng) x 0 時(shí)， y1(x = 5 ， y2 ( x = 0.1 ， y3 ( x = 1 ，積分區(qū)間 是 0,10 ，我們用表 1 列出 GP-穩(wěn)定的二階方法求解例 2。 例 2 是文獻(xiàn)6的 348 頁(yè)中傳染病模型的例子， 我們?cè)谶@ 里增加了-100y3(x這一項(xiàng)，使問(wèn)題成為剛性的。表 4 列出了 例 2 的計(jì)算結(jié)果。 （上接第 1754 頁(yè)） (