3.3數(shù)值積分法穩(wěn)定性分析,在進(jìn)行數(shù)字仿真時(shí)，常常會(huì)出現(xiàn)這樣一種情況：一個(gè)系統(tǒng)本來是穩(wěn)定的，可是仿真結(jié)果是發(fā)散的。這種情況通常是由積分步長選得不合適造成的。那么，為什么計(jì)算步長選得不適會(huì)引起數(shù)值解不穩(wěn)定呢?這就需要分析各種數(shù)值解法的穩(wěn)定性。,3.3.1 數(shù)值解法穩(wěn)定性含義,由于動(dòng)力學(xué)系統(tǒng)的數(shù)值仿真模型的初始數(shù)據(jù)和仿真模型的運(yùn)行過程中都引進(jìn)誤差（如初始誤差、舍入誤差、截?cái)嗾`差等）。這些誤差在模型的運(yùn)行過程中都要向后傳播，影響以后的仿真結(jié)果。如果這些誤差對(duì)以后的影響不是無限增長的。則稱這個(gè)算法對(duì)所用的步長是計(jì)算穩(wěn)定的，否則稱為計(jì)算不穩(wěn)定的。,首先看一個(gè)例子，考慮如下一階系統(tǒng),要求采用Euler法求其數(shù)值解。設(shè)計(jì)算步長為h，則前向Euler遞推公式為：,當(dāng),時(shí)，,，遞推結(jié)果顯然是發(fā)散的；,(2)當(dāng),時(shí)，數(shù)值解等幅振蕩；,(3)當(dāng),時(shí)，遞推結(jié)果是收斂的。,現(xiàn)在分析一下原系統(tǒng)的物理含義。首先它是穩(wěn)定的。其時(shí)間常數(shù)是0.1。其解析解為：,用Euler法對(duì)它進(jìn)行仿真，當(dāng)步長,是由于積分步長太大，從而引起截?cái)嗾`差O(h2)過大造成的。,正好是系統(tǒng)時(shí)間常數(shù)的2倍。,時(shí)之所以不穩(wěn)定，,另外，上述的臨界步長,為了使動(dòng)力學(xué)系統(tǒng)的數(shù)值仿真模型的運(yùn)行結(jié)果能夠反映實(shí)際系統(tǒng)的運(yùn)行過程。必須研究數(shù)值仿真模型的計(jì)算穩(wěn)定性問題，以防止計(jì)算誤差淹沒真實(shí)的運(yùn)動(dòng)狀態(tài)。特別要研究計(jì)算穩(wěn)定性與積分步長之間的關(guān)系，給出建立動(dòng)力學(xué)系統(tǒng)數(shù)字仿真模型時(shí)選擇步長的一些準(zhǔn)則。,從前面論述可知，微分方程(組)的數(shù)值解法，實(shí)質(zhì)上就是將微分方程差分化，然后從初值開始進(jìn)行迭代運(yùn)算。不同的數(shù)值解法對(duì)應(yīng)著不同的遞推公式。一個(gè)數(shù)值法是否穩(wěn)定取決于該差分方程的特征根是否滿足穩(wěn)定性要求。,3.3.2穩(wěn)定性分析,以前向Euler法、后向Euler法、梯形公式為例說各種數(shù)值積分方法穩(wěn)定性分析方法。在前幾節(jié)的論述中已知：,(1)前向歐拉法公式,(2)后向歐拉法公式,(1),(2),(3)梯形公式,(3),下面討論上述三式的穩(wěn)定性。以檢驗(yàn)方程,，,為例進(jìn)行討論。,解析解：,(1)前向歐拉公式,要使上述差分方程穩(wěn)定，必須使下式成立,(4),或,圖 前向歐拉法穩(wěn)定區(qū)域,因此，前向歐拉法的穩(wěn)定區(qū)在復(fù)平面上是不等式(4)滿足的點(diǎn),如下圖所示。在仿真運(yùn)行中，選取的步長h應(yīng)使,使計(jì)算穩(wěn)定要求,積分步長h必須小于系統(tǒng)時(shí)間常數(shù)的2倍。,不超,出這個(gè)穩(wěn)定區(qū)域。特別地，當(dāng)系,統(tǒng)有實(shí)根,(,),時(shí)，為了,，即,Xm變成等比級(jí)數(shù)，由級(jí)數(shù)收斂準(zhǔn)則,(2)后向歐拉公式,要是上述差分方程穩(wěn)定，必須使下式成立,(5),由式(5)可見，只要原系統(tǒng)穩(wěn)定，此不等式必然成立。因此，后向歐拉法是恒穩(wěn)定的。,因?yàn)?即,在復(fù)平面的左半平面,(3)對(duì)于梯形公式,按上述推導(dǎo)過程可得出穩(wěn)定條件為,(6),由式(6)可見，只要原系統(tǒng)穩(wěn)定，此不等式必然成立。因此，梯形公式是恒穩(wěn)定的。,(4)其它算法的穩(wěn)定區(qū)域,由給出的曲線可知，除AM1和AM2(隱式Adams法)為恒穩(wěn)法外，其他方法都是條件穩(wěn)定的。就是步長h必須滿足下列不等式,其中,的常數(shù)。,統(tǒng)微分方程或狀態(tài)方程的,為由積分方法確定,相當(dāng)于連續(xù)系,愈大，選取的積分,特征根或閉環(huán)系統(tǒng)的極點(diǎn)。顯然，,步長應(yīng)愈小。上述一些有關(guān)算法穩(wěn)定性的結(jié)論雖然是從,分析檢驗(yàn)方程,用于其他復(fù)雜系統(tǒng)的模型。,的過程的，但它們同樣適,變化的快,3.4 數(shù)值積分法的選擇與計(jì)算步距的確定,為了有效的對(duì)連續(xù)系統(tǒng)進(jìn)行數(shù)字仿真，必須針對(duì)具體問題，合理地選擇算法和計(jì)算步長。這些問題比較復(fù)雜，涉及的因素也比較多，而且直接影響到數(shù)值解的精度、速度和穩(wěn)定性。 能夠做到十分合理地選擇算法和步長并不是一件十分簡單的事情，因?yàn)閷?shí)際系統(tǒng)是干變?nèi)f化的，所以至今尚無一種具體的、確定的、通用的辦法。 一般來說應(yīng)該考慮以下因素：方法本身復(fù)雜程度，計(jì)算量和誤差的大小，步長和易調(diào)整性等。要特別注意穩(wěn)定性的要求，在此給出一些從實(shí)驗(yàn)中獲得的一些經(jīng)驗(yàn)性的方法，可供實(shí)際應(yīng)用參考。,數(shù)值積分方法的選擇應(yīng)從以下幾個(gè)方面加以考慮。 1.精度要求 影響數(shù)值積分精度的因素包括截?cái)嗾`差(同積分方法、方法階次、步長大小等因素有關(guān))，舍入誤差(同計(jì)算機(jī)字長、步長大小、程序編碼質(zhì)量等因素有關(guān))，初始誤差(由初始值準(zhǔn)確度確定)。 當(dāng)步長h取定時(shí)，算法階次越高，截?cái)嗾`差越小；當(dāng)算法階次取定后，多步法（利用前k點(diǎn)信息）精度比單步法高，隱式精度比顯式高。當(dāng)要求高精度仿真時(shí)，可采用高階的隱式多步法，并取較小的步長。但步長h不能太小，因?yàn)椴介L太小會(huì)增加迭代次數(shù)，增加計(jì)算量，同時(shí)也會(huì)加大舍入誤差和積累誤差。 總之，實(shí)際應(yīng)用時(shí)應(yīng)視仿真精度要求合理地選擇方法和階次，并非階次越高、步長越小越好。,3.4.1 積分方法的選擇,2.計(jì)算速度 計(jì)算速度主要取決于每步積分運(yùn)算所花費(fèi)的時(shí)間以及積分的總次數(shù)，每步運(yùn)算量同具體的積分方法有關(guān)。它主要取決于導(dǎo)函數(shù)的復(fù)雜程度以及每步積分應(yīng)計(jì)算導(dǎo)函數(shù)的次數(shù)。在數(shù)值求解中，最費(fèi)時(shí)間的部分往往就是積分變量導(dǎo)函數(shù)的計(jì)算。通常，隱式比顯式慢，高階比低階慢。 為了提高仿真速度，在積分方法選定的前提下，應(yīng)在保證精度的前提下盡可能加大仿真步長，以縮短仿真總時(shí)間。對(duì)于那些對(duì)速度要求特別苛刻的仿真問題，如實(shí)時(shí)仿真，則宜采用實(shí)時(shí)仿真算法。 對(duì)相同的步長h，RK4比四步Adams 預(yù)估校正法慢。而四步Adams預(yù)估校正法又比AB4慢。一般來說，對(duì)于系統(tǒng)階次高、導(dǎo)函數(shù)復(fù)雜、精度要求高的復(fù)雜仿真問題宜采用Adams預(yù)估校正法。,3.數(shù)值解的穩(wěn)定性 保證數(shù)值解的穩(wěn)定性是進(jìn)行數(shù)字仿真的先決條件，否則計(jì)算結(jié)果將失去實(shí)際意義，導(dǎo)致仿真失敗。從前面穩(wěn)定性的分析可知，同階的RK法的穩(wěn)定性比顯式Adams法好，但不如同階次的隱式Adams法好。從數(shù)值解穩(wěn)定性角度考慮，應(yīng)盡量避免使用顯式Adams法。 總之積分方法的選擇具有較大的靈活性，要結(jié)合實(shí)際問題而定。當(dāng)導(dǎo)函數(shù)不是十分復(fù)雜而且要求精度不是很高時(shí)，RK法是合適的選擇；如果導(dǎo)函數(shù)復(fù)雜、計(jì)算量大，則最好采用Adams預(yù)估校正法；對(duì)于那些實(shí)時(shí)仿真問題，則必須采用實(shí)時(shí)仿真算法。,3.4.2 積分步長的確定,由前面的分析可知，數(shù)字仿真中，步長的選取是一個(gè)十分關(guān)鍵的因素。步長太大則會(huì)導(dǎo)致較大的截?cái)嗾`差，甚至?xí)霈F(xiàn)數(shù)值解的不穩(wěn)定現(xiàn)象；步長太小，又勢必增加計(jì)算次數(shù)，無形中造成舍入誤差的積累，使總誤差加大?？傊?，步長既不能太大也不能很小。一般來說數(shù)字仿真的總誤差不是步長的單調(diào)函數(shù)，而是一個(gè)具有極值的函數(shù)，如下圖所示。,用E表示總誤差，它與步長h之間的函數(shù)關(guān)系為e，則存在一個(gè)最佳步距h0，使得總誤差最小，即,H為步長h的一個(gè)合理的取值區(qū)間。,以二階線性振蕩系統(tǒng)的仿真為例：,取,求解時(shí)間為10s，表3-5列出幾種方法的精度 比較，,其中N為每周計(jì)算點(diǎn)數(shù)，周期為6.28s。,表 四種計(jì)算方法的精度比較,在實(shí)際仿真中，對(duì)于那些變化比較平穩(wěn)的慢變量，步長的改變對(duì)它們的積分總誤差的影響并不很明顯。然而對(duì)于那些變化劇烈的快變量來說，當(dāng)對(duì)它們進(jìn)行積分運(yùn)算時(shí)所產(chǎn)生的總誤差對(duì)步長的改變卻很敏感。因而，在確定積分方法以后，選擇積分步長時(shí)，需要考慮的一個(gè)重要的因素就是系統(tǒng)的動(dòng)態(tài)響應(yīng)特性，對(duì)變化劇烈的快變量，不僅要選擇高階的計(jì)算方法，而且要取較小的積分步長。為了保證計(jì)算穩(wěn)定性，步長只需限制在系統(tǒng)最小時(shí)間常數(shù)的數(shù)量級(jí)，但是為了保證足夠的仿真精度，把積分總誤差控制在一個(gè)較小的范圍內(nèi)，實(shí)際選用的積分步長要比系統(tǒng)最小時(shí)間常數(shù)小得多。 根據(jù)經(jīng)驗(yàn)，對(duì)于一般工程系統(tǒng)的仿真，若采用RK4法，為保證計(jì)算精度在0.5左右，可采用如下經(jīng)驗(yàn)公式確定步長。,其中,為系統(tǒng)在階躍函數(shù)作用下的過渡過程時(shí)間；,為系統(tǒng)開環(huán)頻率特性的剪切頻率；,為系統(tǒng)的最小時(shí)間常數(shù)。,確定相應(yīng)的,當(dāng)系統(tǒng)中有多個(gè)閉合回路時(shí)，按反應(yīng)最快的那個(gè)閉環(huán)系統(tǒng),，,，,。,但高階仿真系統(tǒng)的,，,，,有時(shí)是很難估計(jì)的，或者,是由于系統(tǒng)的非線性，或者有時(shí)根本就無法估計(jì)上述性能指標(biāo)。另外，系統(tǒng)中最小時(shí)間常數(shù)對(duì)應(yīng)的極點(diǎn)只影響到過渡過程起始段形態(tài)，而系統(tǒng)過渡過程主要由那些靠近虛軸的主導(dǎo)極點(diǎn)所決定。固定步長積分方法的計(jì)算步長是按起始段來選取的，這就不可避免地造成后面階段由于使用過小步長積分而引起計(jì)算量增加和時(shí)間浪費(fèi)。為此可以采取一些變步長策略：,分段變步長，將過渡過程分成幾段，每段使用不同的 步長； 根據(jù)每步積分的誤差，自動(dòng)調(diào)整下一步的積分步長； 最優(yōu)步長法，即使每一步積分步長在保證精度的前提下取最大的步長。 積分步長的自動(dòng)改變往往是通過對(duì)誤差的估計(jì)來進(jìn)行的，為此先討論一下誤差估計(jì)方法，然后再介紹相應(yīng)的變步長策略。,3.4.3 誤差估計(jì)與步長控制,因?yàn)檎`差在計(jì)算的每一步都產(chǎn)生，所以對(duì)誤差的估計(jì)與控制也必須在每一步都進(jìn)行。原則上是估算出每一步的局部截?cái)嗾`差，并設(shè)法使它保持在某一容許值內(nèi)?？刂普`差的辦法通常有三種：,(1),(2),(3),其中,是m步的局部截?cái)嗾`差，,是誤差容許值，,是第m,步的步長，xm是第m步的計(jì)算值。一旦獲得局部誤差的估值，就可以利用它來調(diào)整步長。,為了得到每一步局部誤差的估計(jì)值，在實(shí)踐中常使用兩種不同階次的公式，同時(shí)計(jì)算xm+1，并對(duì)得到的兩個(gè)xm+1進(jìn)行比較，根據(jù)比較結(jié)果自動(dòng)改變步長。Merson在1957年提出的RKM(RungeKuttaMerson)方法是最典型的利用誤差估值來控制步長的變步長方法。,Merson首先給出一個(gè)四階RK公式,(1),用上面的RK系數(shù)還可以導(dǎo)出一個(gè)3階的RK公式,由式(1)和式(2)兩式確定的誤差為,(2),(3),該算法四階精度，三階估計(jì)誤差，因此被稱為RKM34法，是目前廣泛采用的一種數(shù)值積分方法。其絕對(duì)穩(wěn)定域和RK4法相近似。該法的缺點(diǎn)是計(jì)算量大，每次需計(jì)算5次導(dǎo)函數(shù)，比普通的RK4法多出1/4的計(jì)算量。,除了采用以上絕對(duì)誤差,誤差做為每步計(jì)算的局部誤差,以外，通常還使用以下的相對(duì),雖然，當(dāng),很大時(shí)，,是相對(duì)誤差；而當(dāng)