第四章連續(xù)系統(tǒng)的離散化方法

第四章連續(xù)系統(tǒng)的離散化方法第一頁，共三十八頁，2022年，8月28日其一般公式為稱為截斷誤差例4－1用歐拉法求下述微分方程的數(shù)值解。第二頁，共三十八頁，2022年，8月28日解：因歐拉法的遞推公式為所以則遞推公式為：若取步長由t=0開始計算可得上式的精確解是數(shù)值計算與精確解的比較見表t00.10.20.31.0精確解x(t)10.90909090.83333330.76923070.510.90.8190.57190.4627810第三頁，共三十八頁，2022年，8月28日4.1.3龍格－庫塔法將在點處作臺勞級數(shù)展開，并取線性部分可得

第四頁，共三十八頁，2022年，8月28日將

代入式比較各項系數(shù)得待定系數(shù)個數(shù)超過方程個數(shù)，必須先設(shè)定一個系數(shù)，然后即可求得其參數(shù)。一般有以下幾種取法：1、則一般形式

第五頁，共三十八頁，2022年，8月28日2、則

一般形式3、則

一般形式以上幾種遞推公式均稱為二階龍格－庫塔公式。是較典型的幾個常用算法。其中的方法3又稱為預估－校正法，或梯形法。第六頁，共三十八頁，2022年，8月28日其意義如下：用歐拉法以斜率先求取一點，再由此點求得另一斜率然后，從點開始，既不按該點斜率變化，也不按預估點斜率變化，而是取兩者平均值求得校正點，即：。第七頁，共三十八頁，2022年，8月28日四階龍格－庫塔法的計算公式為：對于用狀態(tài)方程表示的高階線性系統(tǒng)

其中狀態(tài)變量為用四階龍格－庫塔法時，有計算公式為：第八頁，共三十八頁，2022年，8月28日其中，相應的輸出為按上式，取不斷遞推，

即可求得所需時刻的狀態(tài)變量和輸出值德國學者Felhberg對傳統(tǒng)的龍格－庫塔法提出了改進，在每一個計算步長內(nèi)對f()函數(shù)進行了六次求值，以保證更高的精度和數(shù)值穩(wěn)定性，假設(shè)當前的步長為，定義下面6個變量下一步的狀態(tài)變量可由下式求出：第九頁，共三十八頁，2022年，8月28日

四階/五階龍格－庫塔法系數(shù)表016/13525/2161/41/4003/83/329/326656/128251408/256512/131932/2197-7200/21977296/219728561/564302197/41041439/216-83680/513-845/4104-9/50-1/51/2-8/272-3544/25651859/4104-11/402/550這一方法又稱為四階/五階龍格－庫塔法。第十頁，共三十八頁，2022年，8月28日4.1.4微分方程數(shù)值解的MATLAB實現(xiàn)1、該指令適用于一階常微分方程組如遇到高階常微分方程，必須先將他們轉(zhuǎn)換成一階微分方程組，即狀態(tài)方程方可使用。2、輸入?yún)?shù)為定義微分方程組

M－函數(shù)文件名，可以在文件名加寫@，或用英文格式單引號界定文件名。3、在編輯調(diào)試窗口中編寫一階常微分方程組的M－函數(shù)文件時，每個微分方程的格式必須與一致，即等號嚴格以“先自變量t，后函數(shù)”的固定順序輸入，表示微分方程的序數(shù)。左邊為帶求函數(shù)的一階導數(shù)，右邊函數(shù)的變量4、輸入?yún)?shù)“Tspan”規(guī)定了常微分方程的自變量取值范圍，它以矩陣[t0,tf]的形式輸入，表示自變量5、輸入?yún)?shù)x0表示初始條件向量，

第十一頁，共三十八頁，2022年，8月28日微分方程組中的方程個數(shù)必須等于初始條件數(shù)，這是求微分方程特解所必須的條件。6、輸入?yún)?shù)options表示選項參數(shù)（包括tol，trace），可缺省，即取默認值，tol是控制結(jié)果精度的選項對ode23()函數(shù)取，對ode45()函數(shù)取。trace為輸出形式控制變量，如果trace不為0，則

會將仿真中間結(jié)果逐步地由頻幕顯示出來，否則將不

顯示中間結(jié)果7、輸出參數(shù)[t,x]為微分方程組解函數(shù)的列表（t和x都是列矩陣），它包含向量t各節(jié)點和與對應向量x的第j個分量值（即第j個方程解），i表示節(jié)點序列數(shù)。8、輸出參數(shù)[t,x]缺省時，輸出解函數(shù)的曲線，即函數(shù)及其各的曲線。階導數(shù)求解微分方程的指令還有ode113（多步解法器），ode15s（基于數(shù)字微分公式的解法器），ode23s（單步解法器），ode23T（梯形規(guī)則的一種自由插值實現(xiàn)），ode23TB（二階隱式龍格－庫塔公式）等。第十二頁，共三十八頁，2022年，8月28日例4－1求解常微分方程

，初始條件為：解：方法1把二階微分方程化成兩個一階微分方程組：令

則：

首先編制M文件，并且函數(shù)名和M文件名相同。

functionxdot=wffc_1(t,x)%定義輸入、輸出變量和函數(shù)文件名xdot=zeros(2,1);%明確xdot的維數(shù)xdot(1)=x(1);%第一個微分方程表示形式xdot(2)=-x(1)+2+t^2/pi;%第二個微分方程表示形式方法2寫出系統(tǒng)的狀態(tài)方程第十三頁，共三十八頁，2022年，8月28日首先編制M文件，并且函數(shù)名和M文件名相同。functionxdot=wffc_1(t,x)%定義t,x,xdot和文件名xdot=[01;-10]*x+[0;1]*(2+t^2/pi);%狀態(tài)方程的表示形式在命令窗口鍵入[t,x]=ode45(@wffc_1,[0,10],[-1;1])，可得微分方程的數(shù)值解，其前10組數(shù)據(jù)如下：t=00.01670.03350.05020.06700.15070.23440.31820.40190.5964------x=

-1.00001.0000-0.98281.0501-0.96481.0999-0.94601.1494-0.92631.1986-0.81581.4395-0.68561.6709-0.53631.8917-0.36912.10070.08302.5349------第十四頁，共三十八頁，2022年，8月28日例4－2求VarderPol微分方程，在初始條件下的數(shù)值。解解：將二階微分方程變換成一階微分方程組則

第十五頁，共三十八頁，2022年，8月28日編制M文件，并且函數(shù)名和M文件名相同。functionxdot=vdpl(t,x)xdot=zeros(2,1);％賦初值，并規(guī)定向量的維數(shù)。xdot(1)=x(2);％對第一個微分方程進行描述xdot(2)=2*(1-x(1)^2)*x(2)-x(1);％對第二個微分方程進行描述在命令窗口鍵入[t,x]=ode45('vdpl',[0,20],[1;1])，可得微分方程的數(shù)值解，其前10組數(shù)據(jù)如下：t=00.05020.10050.15070.20100.31360.42620.53890.65150.7484------x=1.00001.00001.04890.94371.09460.87621.13680.79951.17480.71581.24410.51571.29080.31561.31590.13111.3215-0.02781.3131-0.1431------第十六頁，共三十八頁，2022年，8月28日若在命令窗口輸入ode45(@wffc_1,[0,10],[-1;1])，則可得到系統(tǒng)狀態(tài)曲線，而不輸出數(shù)據(jù)。如圖4－4所示。該系統(tǒng)也可直接用SIMULINK進行仿真。首先建立SIMULINK仿真模型如圖4－5所示。設(shè)置系統(tǒng)參數(shù)如下，將積分器初值設(shè)置為系統(tǒng)要求的初值，如圖4－6所示。第十七頁，共三十八頁，2022年，8月28日XY示波器參數(shù)設(shè)置如圖4－7所示。第十八頁，共三十八頁，2022年，8月28日仿真時間參數(shù)設(shè)置如圖4－8所示。系統(tǒng)仿真曲線如圖4－9和4－10所示。第十九頁，共三十八頁，2022年，8月28日4.2數(shù)值解法的穩(wěn)定性及選擇原則4.2.1數(shù)值解法的穩(wěn)定性歐拉法對按歐拉公式計算得這是一個一階差分方程，Z變換后得，其特征值為則該算法的穩(wěn)定條件為。算法的穩(wěn)定域如圖4－11所示。

第二十頁，共三十八頁，2022年，8月28日2.梯形法可見只要原微分方程是穩(wěn)定的（），應用梯形公式進行數(shù)值計算時）。因此梯形公式是恒穩(wěn)公式。算法的穩(wěn)定域如圖4－12所示必然穩(wěn)定（。梯形公式也是一種隱式公式，所以隱式算法是一種恒穩(wěn)算法。第二十一頁，共三十八頁，2022年，8月28日4.2.2數(shù)值解法的選擇原則一般來說，選擇數(shù)值計算方法從以下原則考慮：1、精度問題數(shù)值計算的精度主要受下面三類誤差影響。截斷誤差、舍入誤差和積累誤差。其中截斷誤差同數(shù)值算法的階次有關(guān)，階次越高，截斷誤差越小，一般減小步長可減小每一步的截斷誤差。舍入誤差則與計算機字長有關(guān)，字長越長，舍入誤差越小。積累誤差是上述兩類誤差積累的結(jié)果，它同積分時間長短有關(guān)。一般積分步長越小，則積累誤差越大（在一定的積分時間下）。所以，在一定的數(shù)值計算方法下，若從總誤差考慮，必定有一個最佳步長值。2、計算速度問題計算速度決定于計算的步數(shù)以及每一步積分所需的時間，而每一步的計算時間同數(shù)值計算方法有關(guān)。例如四階龍格－庫塔法每一步需要求四個斜率值，花費較多的時間（但精度高）。在積分方法確定時，應在保證一定的計算精度下盡量能選用較大的步長（必須滿足數(shù)值穩(wěn)定的條件），以縮短積分時間（有時往往選用變步長的算法）。3、數(shù)值解的穩(wěn)定性數(shù)值算法實際上就是將微分方程化成差分方程進行求解，對一個穩(wěn)定的微分方程組，經(jīng)過變換得到的差分方程不一定穩(wěn)定。不同的數(shù)值計算方法有不同的穩(wěn)定域。從穩(wěn)定性角度看，隱式比顯式好。第二十二頁，共三十八頁，2022年，8月28日4.3數(shù)值算法中的“病態(tài)”問題4.3.1“病態(tài)”微分方程例4－3，已知系統(tǒng)的狀態(tài)方程為其中

，試用MATLAB仿真

解：系統(tǒng)狀態(tài)方程可寫為：可建立系統(tǒng)的SIMULINK仿真模型如圖4－13所示。其中第一個積分器輸出為，第二個積分器輸出為，第三個積分器輸出為。第二十三頁，共三十八頁，2022年，8月28日第二十四頁，共三十八頁，2022年，8月28日采用四階龍格－庫塔法，選取固定步長的仿真方法，取h=0.01，參數(shù)設(shè)置由圖可見，t=0.1s之后，曲線變化趨于平緩，若仍以h=0.01計算，會使計算時間拖得很長。第二十五頁，共三十八頁，2022年，8月28日由結(jié)果可看出，誤差很大，仿真有較大的失真。第二十六頁，共三十八頁，2022年，8月28日因此，受數(shù)值穩(wěn)定性限制，只能取小步長進行仿真，若系統(tǒng)是一個慢變過程，則計算速度大受影響。究其原因是狀態(tài)方程的系數(shù)矩陣A的特征值差異較大?？汕蟪鱿到y(tǒng)的特征值如下：A=[-2119-20;19-2120;40-40-40]第二十七頁，共三十八頁，2022年，8月28日A=-2119-2019-212040-40-40>>eig(A)ans=-2.0000-40.0000+40.0000i-40.0000-40.0000i三個特征值為：

而系統(tǒng)的特征值在實際系統(tǒng)中反映了動態(tài)過渡過程的作用。

各瞬態(tài)分量時間常數(shù)第二十八頁，共三十八頁，2022年，8月28日“病態(tài)（stiff）”方程。表述如下：一般線性常微分方程組：的系數(shù)矩陣A的特征值具有如下特征

（4－13）則稱式（4－13）為病態(tài)方程，相應的系統(tǒng)稱為病態(tài)系統(tǒng)。第二十九頁，共三十八頁，2022年，8月28日4.3.2“病態(tài)”系統(tǒng)的仿真方法采用自動變步長數(shù)值計算方法

對于例4－3，

第三十頁，共三十八頁，2022年，8月28日4.4連續(xù)系統(tǒng)狀態(tài)方程的離散化連續(xù)定常系統(tǒng)狀態(tài)方程為：連續(xù)系統(tǒng)狀態(tài)方程解的一般形式為：若上式中令采樣時刻的狀態(tài)

令

采用零階保持器，到期間保持器的輸出為恒定值。且等于時刻的采樣值,并記:

令第三十一頁，共三十八頁，2022年，8月28日最后得:于是離散化的狀態(tài)空間方程為：c2d()函數(shù)可以立即得出連續(xù)系統(tǒng)離散化的模型來，該函數(shù)的調(diào)用命令為：例4－4已知連續(xù)系統(tǒng)的狀態(tài)方程為：，用MATLAB求T=0.1時，相應的離散狀態(tài)方程。第三十二頁，共三十八頁，2022年，8月28日在MATLAB命令窗口輸入：>>A=[2-1-1;0-10;021];>>B=[7;2;3];>>C=[124];>>D=0;>>[G,H]=c2d(A,B,0.1)G=1.2214-0.1162-0.116200.9048000.20031.1052H=0.74730.19030.3355

則離散狀態(tài)方程為：第三十三頁，共三十八頁，2022年，8月28日MATLAB還提供了c2dm()函數(shù)來作類似的變換，其調(diào)用格式為：[G,H,Cd,Dd]=c2dm(A,B,C,D,T,‘選項’)它與c2d函數(shù)的區(qū)別在于，它可以允許用戶自己選擇變換方法。下面采用兩種不同的方法進行離散化：第一種：>>[G,H,Cd,Dd]=c2dm(A,B,C,D,0.1,'zoh')G=1.2214-0.1162-0.116200.9048000.20031.1052H=0.74730.19030.3355Cd=124Dd=0第二種[G,H,Cd,Dd]=c2dm(A,B,C,D,0.1,'tustin')G=1.2222-0.1170-0.117000.9048000.20051.1053H=7.48541.90483.3584Cd=0.11110.22470.4152Dd=1.2364可看出用零階保持器的離散化得出的結(jié)果和直接調(diào)用c2d()函數(shù)所得的結(jié)果相同。第三十四頁，共三十八頁，2022年，8月28日MATLAB的控制工具箱還給出了離散狀態(tài)方程轉(zhuǎn)化為連續(xù)狀態(tài)方程的函數(shù)d2c()，其調(diào)用格式為：

[A,B]=d2c(G,H,T)如對上面采用c2d()得到的結(jié)果，采用如下的變換可得到原來的連續(xù)狀態(tài)方程[A,B]=d2c(G,H,0.1)A=2.0000-1.0000-1.00000-1.0000002.00001.0000B=7.00002.00003.0000第三十五頁，共三十八頁，2022年，8月28日除了d2c()函數(shù)外，MATLAB還提供了d2cm()函數(shù)，其調(diào)用格式為：

[A,B,C,D]=d2cm(G,H,Cd,Dd,T,‘選項’)其中轉(zhuǎn)換方法的選項意義如表4－2所述。對上面采用Tustin變換得到的結(jié)果，采用如下的變換可得到原來的連續(xù)狀態(tài)方程[A,B,C,D]=d2cm(G,H,Cd,Dd,0.1,'tustin')A=2.0000-1.0000-1.00000-1.0000002.00001.0000B=7.00002.00003.0000C=1.00002.00004.0000D=-2.2204e-016由轉(zhuǎn)換結(jié)果可看出，除了在計算中產(chǎn)生了微小誤差外，所得的結(jié)果和原連續(xù)系統(tǒng)相同。第三十六頁，共三十八頁，2022年，8月28日下面介紹一種間接的變換方法，首先將要轉(zhuǎn)換的連續(xù)傳遞函數(shù)轉(zhuǎn)換成連續(xù)的狀態(tài)方程模型，再對該連續(xù)狀態(tài)方程離散化，然后將相應的離散狀態(tài)方程再轉(zhuǎn)換成傳遞函數(shù)，即得相應的脈沖傳遞函數(shù)?？删幹迫缦碌某绦?，并命名為：tfc2d.m可實現(xiàn)由連續(xù)傳遞函數(shù)表示的系統(tǒng)轉(zhuǎn)換成脈沖傳遞函數(shù)表示的離散系統(tǒng)。function[nd,dd