系統(tǒng)時(shí)間響應(yīng)及其仿真PPT學(xué)習(xí)教案

1、會(huì)計(jì)學(xué)1 系統(tǒng)時(shí)間響應(yīng)及其仿真系統(tǒng)時(shí)間響應(yīng)及其仿真 00 ( , ) ( ) dy f t y dt y ty 01N atttb 012N y y yy ( )0,1,2 kk yy tkN 第1頁/共30頁 )1 ( )( ),(,( 0 yty btatytf dt dy o 將區(qū)間a,b分成 N個(gè)小區(qū)間,時(shí)間間隔h（ )也稱為積分步長,在第k個(gè)間隔 t=tk,tk+1 內(nèi)積分： N ab h )2(),( 1 1 k k t t kk dtytfyy 則可用yk（k=0,1,N）作為解y(t)的近似 值，如圖所示。 a b t k y0 yk y t 數(shù)值積分圖解 tk+1 第2頁/共

2、30頁 式(3)是一個(gè)遞推公式。積分值與實(shí)際微分方程解的誤差取決于 步長h和計(jì)算所用的階數(shù)，它是數(shù)值積分的基礎(chǔ)。 ), 2 , 1 , 0( ) 3( 3)3( ! 3 1 2)2( ! 2 1 1 Nk hyhyhyyy kkkkk )2(),( 1 1 k k t t kk dtytfyy 第3頁/共30頁 n顯示法和隱式法 顯示法在計(jì)算yk+1時(shí)所需數(shù)據(jù)均已算出；隱式法在計(jì)算yk+1 時(shí)需用到tk+1時(shí)刻的數(shù)據(jù)，該算法必須借助予估公式。 n定步長和變步長 定步長為積分步長在仿真運(yùn)行過程中始終不變；變步長指在仿真 運(yùn)行過程中自動(dòng)修改步長。 第4頁/共30頁 t0t1t2t3 h y( t)

3、 y 0 y 1 y t 歐拉近似 解 歐拉法 n【說明】歐拉法是用一條過各點(diǎn) 的切線取代曲線來逼近精確解。 該算法簡單，計(jì)算量小，但精度 較低。 )4( 1mmm hyyy )5( ),( ),( ),( 1 1112 00001 0 0 hytfyy hytfyy hytfyhyyy mmmm t y ) 3( 3)3( ! 3 1 2)2( ! 2 1 1 kkkkk hyhyhyyy 第5頁/共30頁 ) 3( 3)3( ! 3 1 2)2( ! 2 1 1 kkkkk hyhyhyyy ),(),( 11 2 1 1 mmmmmm ytfytfhyy )6( ),(),( )(),(

4、 112 1 1 1 （校正） 予估 mmmmmm mmmm ytfytfhyy hytfyy 第6頁/共30頁 ) 3( 3)3( ! 3 1 2)2( ! 2 1 1 kkkkk hyhyhyyy 第7頁/共30頁 ) 3( 3)3( ! 3 1 2)2( ! 2 1 1 kkkkk hyhyhyyy )7(),(),( 2 2 1 2)2( 2 1 1 hytfythfy hyhyyy my f mmmt f mmm mmmm 設(shè)原微分方程(1)式解具有以下形式: )8( )( ),( ),( 22111 1212 1 KaKahyy hKbyhbtfK ytfK mm mm mm 式中

5、，a1，a2，b1，b2為待定系數(shù)。 ),( ),( )2( ytf dt d y ytf dt dy 第8頁/共30頁 )7(),(),( 2 2 1 1 hytfythfyy my f mmmt f mmmm )8( )( ),( ),( 22111 1212 1 KaKahyy hKbyhbtfK ytfK mm mm mm 將(8)式中K2按二元函數(shù) 展開成泰勒級(jí)數(shù)，并取前三 項(xiàng) )9(),( ),( 121 1212 my f mt f mm my f mt f mm Kbbhytf hKbhbytfK 將K1，K2代入(8)式： )10(),(),()( ),(),(),( 21

6、2 221 21211 my f mmmt f mmm my f mmmt f mmmmmm ytfbbhaythfaay ythfbhbytfhaythfayy 比較(6-10)、(6-7)式： )11( 2/1 2/1 1 22 12 21 ba ba aa 第9頁/共30頁 顯然由(11)式并不能唯一確定a1,a2,b1,b2 ，因?yàn)橹挥腥齻€(gè)方程。因此對于同一種算 法可以有不同的表現(xiàn)形式。 n【說明】由于該算法只取到泰勒展開式的二階導(dǎo)數(shù)項(xiàng)，所以稱為 二階龍格-庫塔法。但由(8)(12)式可知，算法并沒有用y的二階 導(dǎo)數(shù)。 若設(shè)a1=a2，則 1 2/1 21 21 bb aa 即二階RK

7、法公式為 )12( )( 2 ),( ),( 211 12 1 KK h yy hKyhtfK ytfK mm mm mm )11( 2/1 2/1 1 22 12 21 ba ba aa 第10頁/共30頁 )13( , 2 , 1 ),( 1 1 1 1 ri KbyhatfK Khyy i j jijmimi i r i imm 式中，i為待定權(quán)系數(shù)，ai,bij為待定系數(shù)，r為使用Ki的個(gè)數(shù)(即級(jí)數(shù)) ，Ki為所取各點(diǎn)導(dǎo)數(shù)f的值。 Ki的個(gè)數(shù)與yk+1泰勒展開式所取的項(xiàng)數(shù)有 關(guān)(即RK算法的階數(shù))，同時(shí)還與計(jì)算區(qū)間內(nèi)所取導(dǎo)數(shù)值的點(diǎn)數(shù)有關(guān)。 第11頁/共30頁 )14( ),( ),(

8、),( ),( 22 34 22 1 23 12 1 22 1 432161 hKyhtfK hKytfK hKytfK ytfK KKKKyy mm m h m m h m mm h mm 由于(14)式在同級(jí)的RK算法中，計(jì)算精度較高，計(jì)算量較少，而 在系統(tǒng)仿真的數(shù)值積分中應(yīng)用十分廣泛。稱之為四階四級(jí)RK公式。 第12頁/共30頁 00 ( )( )( ),( )X tAX tBU tX tX 則稱為“病態(tài)”方程。 )17(50 Remin Remax ), 2 , 1(0Re i i ii ni 第13頁/共30頁 取計(jì)算步長與計(jì)算精度、計(jì)算時(shí)間之 間的矛盾。 第14頁/共30頁 )18

9、( 2 1 1 2 1 , 2 1 1 2 1 2 1 , 2 2 , 2 , 324 213 1 2 1 KKyhthfK KKy h thfK K y h thfK ythfK mm mm mm mm )19( 6 1 2 1 1 3 1 2 1 1 3 1 6 1 43211 KKKKyy mm 第15頁/共30頁 n局部截?cái)嗾`差局部截?cái)嗾`差：由積分方法和階次 的限制而引起的誤差。這種誤差與 h成正比。 截?cái)嗾`差 舍入誤差 總誤差 e h 誤差與積分步長 顯然選擇一個(gè)合適的積分步長可使 總誤差達(dá)到最小。 第16頁/共30頁 n此外，h還應(yīng)與模型的信號(hào)響 應(yīng)情況有關(guān)，例如在穩(wěn)態(tài)時(shí)， 可取較

10、大的步長，見上圖。 6.283212.566418.8496 0.95 1 1.05 1.2843 y = 1 - e - tcos t =0.3 =0.7 鎮(zhèn)定時(shí)間 ts = 9.6133 t y 變步長：變步長：在仿真計(jì)算過程 中根據(jù)計(jì)算誤差的大小來 改變步長。其目的是在保 證一定計(jì)算精度的前提下 ，盡可能選擇較大步長。 第17頁/共30頁 T,Y=ode45(f,tspan,y0) n【說明1】 f 為常微分方程（組）或系統(tǒng)模型的文件名； tspan=t0,tfinal 即積分時(shí)間初值和終值； y0是積分初值； T為計(jì)算時(shí)間點(diǎn)的時(shí)間向量； Y為相應(yīng)的微分方程解數(shù)據(jù)向量或矩陣。 第18頁/

11、共30頁 n【說明2】 對于剛性微分方程（特征值數(shù)值相差較大），可用ode15s，其調(diào)用格式與ode45相同。 ode函數(shù)只能用于求解一階微分方程或一階微分方程組。若系統(tǒng)的數(shù)學(xué)模型為高階微分方程，則應(yīng)將高階微分方程轉(zhuǎn)化成一階微分方程組。因此在用MATLAB的ode函數(shù)求解微分方程時(shí)，應(yīng)首先建立描述系統(tǒng)模型的一階微分方程（組）函數(shù)f 。 第19頁/共30頁 【例2】已知二階 微分方程 求時(shí)間區(qū)間t=0，20微分方程的解。 1)0(, 0)0( 0)1 ( 2 yy yyyy 解：解：（1）將微分方程表示為一階微分方程組 12 2 12 21 1 )1 (yyyy yy yy n【說明】這種描述系

12、統(tǒng)微分方程的函數(shù)與ODE函數(shù)配套使用， 其格式是固定的。dy為2*1數(shù)組，其維數(shù)等于微分方程的階數(shù)。 （2）建立描述系統(tǒng)微分方程的m-函數(shù)文件vdp.m function dy=vdp(t,y) dy=zeros(2,1); % 生成2行1列的零陣 dy(1)= y(2); % dy(2)= (1-y(1)2)*y(2)-y(1); % 21 yy 12 2 12 )1 (yyyy 第20頁/共30頁 （3）編寫MATLAB主程序 T,Y=ode45(VDPd,0 20,0,1); %調(diào)用ode45產(chǎn)生離散點(diǎn)時(shí)間向量和解向量 plot(T,Y(:,1),r-,T,Y(:,2),b:) titl

13、e(Solution) xlabel(time s),ylabel(Position Y) legend(y1,y2) 運(yùn)行結(jié)果如右圖 所示。其中y1(紅線)為 微分方程的解。 第21頁/共30頁 step(sys) step(sys,Tfinal) 其中，sys為系統(tǒng)模型（傳遞函數(shù)）；Tfinal為仿真終止 時(shí)間，若省略則由系統(tǒng)默認(rèn)。 第22頁/共30頁 【例4】已知系統(tǒng)模型 ，求其單位階躍響應(yīng)。 5 1 )( 2 ss s sG sys=tf(1,-1,1,1,5) subplot(1,2,1),step(sys,20) subplot(1,2,2),step(sys) 建立系統(tǒng)模型 指定

14、階躍響應(yīng)時(shí)間 不指定階躍響應(yīng)時(shí) 間 第23頁/共30頁 step(sys1,sys2,) step(sys1,r,sys2,y-,sys3,gx) Y,T=step(sys) 第24頁/共30頁 impulse(sys) impulse(sys,Tfinal) impulse(sys1,sys2,) Y,T=impulse(sys) 第25頁/共30頁 U,T=gensig(Type,Tau) U,T=gensig(Type,Tau,Tf,Ts) 第26頁/共30頁 lsim(sys,U,T) %基本調(diào)用格式 lsim(sys1,sys2,U,T) %繪制多個(gè)系統(tǒng)對同一個(gè)任意輸入響應(yīng)曲線 Ys,Ts=lsim(sys,U,T) %該格式不繪制響應(yīng)曲線 第27頁/共30頁 【例7】已知系統(tǒng)模型 ，計(jì)算系統(tǒng)在周期為5s的方波信號(hào)作用下的響應(yīng)。 sys=tf(3,100,1,10,40,100); u,t=