差分方程的解法分析及MATLAB實現(xiàn)程序

1、差分方程的解法分析及MATLA聯(lián)現(xiàn)(程序)摘自：張登奇，彭仕玉.差分方程的解法分析及其MATLAB實現(xiàn)J.湖南理工學(xué)院學(xué)報.2014(03)引言線性常系數(shù)差分方程是描述線性時不變離散時間系統(tǒng)的數(shù)學(xué)模型，求解差分方程是分析離散時間系統(tǒng)的重要內(nèi)容.在信號與系統(tǒng)課程中介紹的求解方法主要有迭代法、時域經(jīng)典法、雙零法和變換域法山1迭代法例1已知離散系統(tǒng)的差分方程為y(n)3y(n1)+4y(n2)=x(n)+x(n1),激勵信號為4o3x(n)=(3)nu(n),初始狀態(tài)為y(-1)=4,y(0=12.求系統(tǒng)響應(yīng).根據(jù)激勵信號和初始狀態(tài)，手工依次迭代可算出y(0)=2,y(1)=蜀.利用MATLABHf

2、ilter函數(shù)實現(xiàn)迭代過程的m程序如下：clc;clear;formatcompact;a=1,-3/4,1/8,b=1,1/3,0,%輸入差分方程系數(shù)向量，不足補(bǔ)0對齊n=0:10;xn=(3/4).An,%zx=0,0,zy=4,12,%zi=filtic(b,a,zy,zx),%輸入激勵信號輸入初始狀態(tài)計算等效初始條件yn,zf=filter(b,a,xn,zi),%迭代計算輸出和后段等效初始條件2時域經(jīng)典法用時域經(jīng)典法求解差分方程：先求齊次解；再將激勵信號代入方程右端化簡得自由項，根據(jù)自由項形式求特解；然后根據(jù)邊界條件求完全解3.用時域經(jīng)典法求解例1的基本步驟如下.(1)求齊次解.特征

3、方程為U246+8=0,可算出5=2，%=4.高階特征根可用MATLAB勺roots函數(shù)計算.齊次解為yh(n)=C1©)n+C2(4)n,n>0.(2)求方程的特解.將x(n)=(,)nu(n)代入差分方程右端得自由項為(2+1<4)nAu(n-1)=J1,3.(3)n；S.9當(dāng)n±1時，特解可設(shè)為yp(n)=D(3)n,代入差分方程求得D=13.42(3)利用邊界條件求完全解.當(dāng)n=0時迭代求出y(0)=摟，當(dāng)n>1時，完全解的形式為y(n)=G(1)n十C2(：)n+號©>，選擇求完全解系數(shù)的邊界條件可參考文4選y(0),y(1).根

4、據(jù)邊界條件求得C1=-號。2=5.注意完全解的表達(dá)式只適于特解成立的n取值范圍，其他點要用6(n)及其延遲表示，如果其值符合表達(dá)式則可合并處理.差分方程的完全解為v(n)一5(n)17/1)n一5(1)n一13(3)niN(n1)_r17(1)n一5(1)n-13/3niN/n)y(n)2v(n)-3(2)3(4)2(4)u(n')_3(2)3(4)2(4)iu(n)MATLAB沒有專用的差分方程求解函數(shù)，但可調(diào)用maple符號運(yùn)算工具箱中的rsoke函數(shù)實現(xiàn)格式y(tǒng)=maple('rsolve(equs,inis,y(n)'),其中:equs為差分方程表達(dá)式，inis為

5、邊界條件，y(n)為差分方程中的輸出函數(shù)式.rsolve的其他格式可通過mhelprsolve命令了解.在MATLAB用時域經(jīng)典法求解例1中的全響應(yīng)和單位樣值響應(yīng)的程序如下.clc;clear;formatcompact;yn=maple('rsolve(y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)An+1/3*(3/4)A(n-1),y(0)=5/2,y(-1)=4,y(n),hn=maple('rsolve(y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(0)=1,y(1)=13/12,y(n)'),3雙零法根據(jù)雙零響應(yīng)的定義，按時

6、域經(jīng)典法的求解步驟可分別求出零輸入響應(yīng)和零狀態(tài)響應(yīng).理解了雙零法的求解原理和步驟，實際計算可調(diào)用rsolve函數(shù)實現(xiàn).yzi=maple('rsolve(y(n)-3/4*y(n-1)+1/8*y(n-2)=0,y(-1)=4,y(-2)=12,y(n)'),yzs=maple('rsolve(y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)An+1/3*(3/4)A(n-1),y(0)=1,y(-1)=0,y(n)1),4變換域法NM設(shè)差分方程的一般形式為、aky(n-k)brx(n-r).k&r旦對差分方程兩邊取單邊z變換，并利用z變換的位移

7、公式得N1M1、akz*Y(z)二：y(l)z力=.1，brzX(z)x(m)z""k3l-kr蟲m-l整理成A(z)Y(z)+Y°(z)=B(z)X(z)+X0(z)形式有A(z)=a0a1z1，aNz小，B(z)=b0b1z，bMz.N1M-1Y0(s)=、.工akyz”;X0(s)=""brX(m)z，”.k；11-krz1m-_l可以看出，由差分方程可直接寫出A(z)和B(z),系統(tǒng)函數(shù)H(z)=B(z)/A(z),將系統(tǒng)函數(shù)進(jìn)行逆z變換可得單位樣值響應(yīng).由差分方程的初始狀態(tài)可算出K(z),由激勵信號的初始狀態(tài)可算出X°(z)

8、,將激勵信號進(jìn)行z變換可得X(z)，求解夜代數(shù)方程可得輸出信號的象函數(shù)B(z)X(z)X0(z)-Y0(z)Y(z),A(z)對輸出象函數(shù)進(jìn)行逆z變換可得輸出信號的原函數(shù)y(n).利用z變換求解差分方程各響應(yīng)的步驟可歸納如下：(1)根據(jù)差分方程直接寫出A(z)、B(z)和H(z),H(z)的逆變換即為單位樣值響應(yīng)；(2)根據(jù)激勵信號算出X(z),如激勵不是因果序列則還要算出前M個初始狀態(tài)值；(3)根據(jù)差分方程的初始狀態(tài)y(-1),y(-2),-；y(-N)和激勵信號的初始狀態(tài)x(-1),x(-2),7(f)算出Y0(z)和X0(z)；(4)在zM求解代數(shù)方程A(z)Y(z)+Y)(z)=B(z

9、)X(z)+X0(z)得輸出象函數(shù)Y(z),Y(z)的逆變換即為全響應(yīng)；(5)分析響應(yīng)象函數(shù)的極點來源及在z平面中白位置，確定自由響應(yīng)與強(qiáng)迫響應(yīng)，或瞬態(tài)響應(yīng)與穩(wěn)態(tài)響應(yīng);(6)根據(jù)零輸入響應(yīng)和零狀態(tài)響應(yīng)的定義，在城求解雙零響應(yīng)的象函數(shù)，對雙零響應(yīng)的象函數(shù)進(jìn)行逆z變換，得零輸入響應(yīng)和零狀態(tài)響應(yīng).用變換域法求解例1的基本過程如下.根據(jù)差分方程直接寫出A(z)=1-3z-+8z-,B(z)=1+亞，系統(tǒng)函數(shù)的極點為2,1.對激勵信號進(jìn)行z變換得X(z)=z/(z一節(jié).激勵象函數(shù)的極點為3/4.根據(jù)差分方程的初始狀態(tài)算出(z)=3+#z根據(jù)激勵信號的初始狀態(tài)算出X0(z)=0.對z域代數(shù)方程求解，得全響

10、應(yīng)的象函數(shù)Y(z)=(5z3-31z2+3z)/(z3-3z+11z爰).224o21632進(jìn)行逆z變換得全響應(yīng)為y(n)=-號(*)n+5)n+,(4)nu(n)其中，與系統(tǒng)函數(shù)的極點對應(yīng)的是自由響應(yīng)；與激勵象函數(shù)的極點對應(yīng)的是強(qiáng)迫響應(yīng).Y(z)的極點都在z平面的單位圓內(nèi)故都是瞬態(tài)響應(yīng).零輸入響應(yīng)和零狀態(tài)響應(yīng)可按定義參照求解上述求解過程可借助MATLAB勺符號運(yùn)算編程實現(xiàn).實現(xiàn)變換域法求解差分方程的m程序如下：clc;clear;formatcompact;symszn%定義符號對象%輸入差分方程、初始狀態(tài)和激勵信號%a=1,-3/4,1/8,b=1,1/3,%輸入差分方程系數(shù)向量y0=4,

11、12,x0=0,%輸入初始狀態(tài)，長度分別比a、b短1,長度為0時用口xn=(3/4)An,%輸入激勵信號，自動單邊處理，u(n)可用1表示%下面是變換域法求解差分方程的通用程序，極點為有理數(shù)時有解析式輸出%N=length(a)-1;M=length(b)-1;%計算長度Az=poly2sym(a,'z')/zAN;Bz=poly2sym(b,'z')/zAM;%計算A(z)和B(z)Hz=Bz/Az;disp('系統(tǒng)函數(shù)H(z):'),sys=filt(b,a),%計算并顯示系統(tǒng)函數(shù)hn=iztrans(Hz);disp('單位樣值響應(yīng)

12、h(n)='),pretty(hn),%計算并顯示單位樣值響應(yīng)Hzp=roots(a);disp('系統(tǒng)極點:');Hzp,%計算并顯示系統(tǒng)極點Xz=ztrans(xn);disp('激勵象函數(shù)X(z)='),pretty(Xz),%激勵信號的單邊z變換Y0z=0;%初始化Y0(z),求Y0(z)注意系數(shù)標(biāo)號與變量下標(biāo)的關(guān)系fork=1:N;forl=-k:-1;Y0z=Y0z+a(k+1)*y0(-l)*zA(-k-l);endenddisp('初始Y0(z)'),Y0z,%系統(tǒng)初始狀態(tài)的z變換X0z=0;%初始化X0(z),求X0(z

13、)注意系數(shù)標(biāo)號與變量下標(biāo)的關(guān)系forr=1:M;form=-r:-1;X0z=X0z+b(r+1)*x0(-m)*zA(-r-m);endenddisp('初始X0(z)'),X0z,%激勵信號起始狀態(tài)的z變換Yz=(Bz*Xz+X0z-Y0z)/Az;disp(,全響應(yīng)的z變換Y(z)'),pretty(simple(Yz),yn=iztrans(Yz);disp('全響應(yīng)y(n)='),pretty(yn),%計算并顯示全響應(yīng)Yziz=-Y0z/Az;disp('零輸入象函數(shù)Yzi(z)='),pretty(Yziz),%零激勵響應(yīng)的z變換yzin=iztrans(Yziz);disp('零輸入響應(yīng)yzi(n)='