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

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

3、方程為U246+8=0,可算出5=2，%=4.高階特征根可用MATLAB勺roots函數計算.齊次解為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當n±1時，特解可設為yp(n)=D(3)n,代入差分方程求得D=13.42(3)利用邊界條件求完全解.當n=0時迭代求出y(0)=摟，當n>1時，完全解的形式為y(n)=G(1)n十C2(：)n+號©>，選擇求完全解系數的邊界條件可參考文4選y(0),y(1).根

4、據邊界條件求得C1=-號。2=5.注意完全解的表達式只適于特解成立的n取值范圍，其他點要用6(n)及其延遲表示，如果其值符合表達式則可合并處理.差分方程的完全解為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沒有專用的差分方程求解函數，但可調用maple符號運算工具箱中的rsoke函數實現格式y=maple('rsolve(equs,inis,y(n)'),其中:equs為差分方程表達式，inis為

5、邊界條件，y(n)為差分方程中的輸出函數式.rsolve的其他格式可通過mhelprsolve命令了解.在MATLAB用時域經典法求解例1中的全響應和單位樣值響應的程序如下.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雙零法根據雙零響應的定義，按時

6、域經典法的求解步驟可分別求出零輸入響應和零狀態響應.理解了雙零法的求解原理和步驟，實際計算可調用rsolve函數實現.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設差分方程的一般形式為、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),系統函數H(z)=B(z)/A(z),將系統函數進行逆z變換可得單位樣值響應.由差分方程的初始狀態可算出K(z),由激勵信號的初始狀態可算出X°(z)

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

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

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

11、12,x0=0,%輸入初始狀態，長度分別比a、b短1,長度為0時用口xn=(3/4)An,%輸入激勵信號，自動單邊處理，u(n)可用1表示%下面是變換域法求解差分方程的通用程序，極點為有理數時有解析式輸出%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('系統函數H(z):'),sys=filt(b,a),%計算并顯示系統函數hn=iztrans(Hz);disp('單位樣值響應

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

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