張登奇, 彭仕玉

(湖南理工學(xué)院 信息與通信工程學(xué)院, 湖南 岳陽 414006)

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

張登奇, 彭仕玉

(湖南理工學(xué)院 信息與通信工程學(xué)院, 湖南 岳陽 414006)

差分方程是描述離散時間系統(tǒng)的數(shù)學(xué)模型, 求解差分方程是分析離散時間系統(tǒng)的重要內(nèi)容, 常用的求解方法有迭代法、時域經(jīng)典法、雙零法和變換域法. 文章根據(jù)各種方法的求解原理, 分別介紹了不同方法的求解步驟, 結(jié)合實(shí)例列出了這些方法的求解過程及MATLAB實(shí)現(xiàn)程序.

離散時間系統(tǒng); 差分方程; 時域分析法; 變換域分析法; MATLAB

引言

線性常系數(shù)差分方程是描述線性時不變離散時間系統(tǒng)的數(shù)學(xué)模型, 求解差分方程是分析離散時間系統(tǒng)的重要內(nèi)容. 在《信號與系統(tǒng)》課程中介紹的求解方法主要有迭代法、時域經(jīng)典法、雙零法和變換域法[1]. 迭代法可手工逐次代入求解, 也可編程用計(jì)算機(jī)求解, 該方法原理簡單, 缺點(diǎn)是只能得到數(shù)值解.時域經(jīng)典法先求齊次解和特解, 再用邊界條件確定待定系數(shù)得完全解, 該方法數(shù)學(xué)過程清晰, 但求解過程麻煩. 雙零法分別求零輸入響應(yīng)和零狀態(tài)響應(yīng), 再通過疊加得到全響應(yīng). 該方法物理意義清晰, 但求解過程依然麻煩. 變換域法利用z變換將差分方程變換成代數(shù)方程求解, 該方法簡便高效, 是求解差分方程的重要方法. 本文根據(jù)不同方法的求解原理, 分別介紹各種方法的求解步驟, 結(jié)合實(shí)例列出這些方法的求解過程和MATLAB實(shí)現(xiàn)程序.

1 迭代法

差分方程本身就是一個遞推方程, 根據(jù)初始狀態(tài)和激勵信號依次迭代就可算出輸出序列. 迭代法是解差分方程的基礎(chǔ)方法, 如果所需輸出序列個數(shù)較少（如計(jì)算邊界條件）可手工直算, 如需計(jì)算大量輸出可利用計(jì)算機(jī)編程實(shí)現(xiàn).現(xiàn)結(jié)合實(shí)例介紹迭代法的計(jì)算過程.

例1已知離散系統(tǒng)的差分方程為激勵信號為初始狀態(tài)為y(?1)=4,y(?2)=12.求系統(tǒng)響應(yīng).

利用MATLAB中的filter函數(shù)實(shí)現(xiàn)迭代過程的m程序如下:

clc;clear;format compact;

a=[1, -3/4, 1/8], b=[1, 1/3, 0], %輸入差分方程系數(shù)向量, 不足補(bǔ)0對齊

n=0:10;xn=(3/4).^n, %輸入激勵信號

zx=[0, 0], zy=[4, 12], %輸入初始狀態(tài)

zi=filtic(b, a, zy, zx), %計(jì)算等效初始條件

[yn, zf]=filter(b, a, xn, zi), %迭代計(jì)算輸出和后段等效初始條件

MATLAB提供的filter函數(shù)是一個內(nèi)建函數(shù), 用type命令看不到程序代碼.為了理解迭代思想, 下面根據(jù)圖1所示的直接Ⅰ型結(jié)構(gòu)[2], 重寫實(shí)現(xiàn)迭代法的m程序.

clc;clear;format compact;

a=[1, -3/4, 1/8], b=[1, 1/3, 0], %輸入差分方程系數(shù)向量, 不足補(bǔ)0對齊

n=0:10;x=(3/4).^n, %輸入激勵信號

zx=[0, 0], zy=[4, 12], %輸入初始狀態(tài)

% 下面是按直接Ⅰ型結(jié)構(gòu)迭代的通用程序 %

N=length(a)-1, %計(jì)算數(shù)據(jù)存儲長度

L=length(x), %計(jì)算激勵信號長度

y=zeros(1, L);%輸出信號初始化

for i=1:L; %逐個計(jì)算輸出信號

for n=1:N;z(n)=b(n+1)*zx(n)-a(n+1)*zy(n);end %分算過去

zz=sum(z);%計(jì)算輸出中的過去分量

y(i)=b(1)*x(i)+zz;% 計(jì)算當(dāng)前輸出y(n)

for n=N:-1:2, zx(n)=zx(n-1);zy(n)=zy(n-1);end%過去數(shù)據(jù)下移

zx(1)=x(i);zy(1)=y(i);%當(dāng)前的激勵和輸出變?yōu)檫^去, 以便算下一個輸出

end

%理解filter函數(shù)中zf參數(shù)的意義

zf=zeros(1, N);%初始化zf

for k=1:N; %逐個計(jì)算zf參數(shù)

for n=1:N;z(n)=b(n+1)*zx(n)-a(n+1)*zy(n);end %算z(n)

zf(k)=sum(z);% 計(jì)算第k個zf

for k=N:-1:2, zx(k)=zx(k-1);zy(k)=zy(k-1);end;%過去數(shù)據(jù)下移

zx(1)=0;zy(1)=0;% 沒有當(dāng)前的激勵和輸出變?yōu)檫^去, 算下一個zf

end

y, zf, %顯示輸出和zf參數(shù)

該程序采用直接Ⅰ型結(jié)構(gòu)實(shí)現(xiàn)迭代計(jì)算, 與filter函數(shù)計(jì)算的結(jié)果相同, zf參數(shù)可用于激勵信號分段處理時，調(diào)用filter函數(shù)計(jì)算下一段輸出所需的等效初始條件.

圖1 直接Ⅰ型網(wǎng)絡(luò)結(jié)構(gòu)圖

2 時域經(jīng)典法

用時域經(jīng)典法求解差分方程與高等數(shù)學(xué)中求解微分方程的過程類似: 先求齊次解; 再將激勵信號代入方程右端化簡得自由項(xiàng), 根據(jù)自由項(xiàng)形式求特解; 然后根據(jù)邊界條件求完全解[3]. 用時域經(jīng)典法求解例1的基本步驟如下.

當(dāng)n≥1時, 特解可設(shè)為代入差分方程求得

(3) 利用邊界條件求完全解.當(dāng)n=0時迭代求出當(dāng)n≥1時, 完全解的形式為選擇求完全解系數(shù)的邊界條件可參考文[4]選y(0),y(?1). 根據(jù)邊界條件求得注意完全解的表達(dá)式只適于特解成立的n取值范圍, 其他點(diǎn)要用δ(n)及其延遲表示, 如果其值符合表達(dá)式則可合并處理.差分方程的完全解為

MATLAB沒有專用的差分方程求解函數(shù), 但可調(diào)用maple符號運(yùn)算工具箱中的rsolve函數(shù)實(shí)現(xiàn)[5], 格式為y=maple('rsolve({equs, inis}, y(n))'), 其中equs為差分方程表達(dá)式, inis為邊界條件,y(n)為差分方程中的輸出函數(shù)式. rsolve的其他格式可通過mhelp rsolve命令了解. 在MATLAB中用時域經(jīng)典法求解例1中的全響應(yīng)和單位樣值響應(yīng)的程序如下:

clc;clear;format compact;

yn=maple('rsolve({y(n)-3/4*y(n-1)+1/8*y(n-2)=(3/4)^n+1/3*(3/4)^(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))'),

如同時域經(jīng)典法一樣, 應(yīng)用rsolve函數(shù)時要特別注意邊界條件的選擇問題，且邊界條件要連續(xù)取點(diǎn).本例求單位樣值響應(yīng)時的邊界條件要選用y(0)和y(1)[6], 在命令窗中顯示的結(jié)果也要分析n的取值范圍.

3 雙零法

雙零法是將完全解分解成物理意義明顯的零輸入響應(yīng)和零狀態(tài)響應(yīng)分別計(jì)算. 零輸入響應(yīng)是激勵為零, 由系統(tǒng)的初始狀態(tài)所產(chǎn)生的響應(yīng); 零輸入響應(yīng)要求差分方程右端為零, 故特解為零; 完全解為齊次解形式, 系數(shù)可直接由初始狀態(tài)確定. 零狀態(tài)響應(yīng)是初始狀態(tài)為零, 由激勵信號所產(chǎn)生的響應(yīng). 計(jì)算零狀態(tài)響應(yīng)可用時域經(jīng)典法, 也可用卷積法. 根據(jù)雙零響應(yīng)的定義, 按時域經(jīng)典法的求解步驟可分別求出零輸入響應(yīng)和零狀態(tài)響應(yīng). 理解了雙零法的求解原理和步驟, 實(shí)際計(jì)算可調(diào)用rsolve函數(shù)實(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)^n+1/3*(3/4)^(n-1), y(0)=1, y(-1)=0}, y(n))'),

邊界條件和n的取值范圍這里不再贅述, 各響應(yīng)的計(jì)算結(jié)果在命令窗都有顯示, 零輸入響應(yīng)和零狀態(tài)響應(yīng)之和與前面計(jì)算的全響應(yīng)相等.

4 變換域法

變換域法是以z變換為數(shù)學(xué)工具, 先將時域差分方程變換成z域代數(shù)方程, 再在z域求解響應(yīng)的象函數(shù), 最后將響應(yīng)象函數(shù)逆變換成時域原函數(shù), 是一種間接求解差分方程的計(jì)算方法.z變換的線性和位移性是變換域法求解差分方程的重要性質(zhì), 這里只列出雙邊序列的單邊z變換位移公式[1]:

對差分方程兩邊取單邊z變換, 并利用z變換的位移公式, 得

整理成A(z)Y(z)+Y0(z)=B(z)X(z)+X0(z)形式, 有

可以看出, 由差分方程可直接寫出A(z)和B(z), 系統(tǒng)函數(shù)將系統(tǒng)函數(shù)進(jìn)行逆z變換可得單位樣值響應(yīng). 由差分方程的初始狀態(tài)可算出Y0(z), 由激勵信號的初始狀態(tài)可算出X0(z), 將激勵信號進(jìn)行z變換可得X(z), 求解z域代數(shù)方程可得輸出信號的象函數(shù)

對輸出象函數(shù)進(jìn)行逆z變換可得輸出信號的原函數(shù)y(n). 為了計(jì)算簡單和便于對結(jié)果進(jìn)行對比分析,例1列舉的實(shí)例中, 系統(tǒng)函數(shù)的極點(diǎn)均為單根有理數(shù)極點(diǎn), 且系統(tǒng)函數(shù)H(z)的零極點(diǎn)與激勵象函數(shù)X(z)的零極點(diǎn)均不相同. 參照s域求解微分方程的方法[7], 利用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),???,x(?M)算出Y0(z)和X0(z);

(4) 在z域求解代數(shù)方程A(z)Y(z)+Y0(z)=B(z)X(z)+X0(z)得輸出象函數(shù)Y(z),Y(z)的逆變換即為全響應(yīng);

(5) 分析響應(yīng)象函數(shù)的極點(diǎn)來源及在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)的定義, 在z域求解雙零響應(yīng)的象函數(shù), 對雙零響應(yīng)的象函數(shù)進(jìn)行逆z變換, 得零輸入響應(yīng)和零狀態(tài)響應(yīng).

用變換域法求解例1的基本過程如下.

對激勵信號進(jìn)行z變換得激勵象函數(shù)的極點(diǎn)為

對z域代數(shù)方程求解, 得全響應(yīng)的象函數(shù)

進(jìn)行逆z變換得全響應(yīng)為

其中, 與系統(tǒng)函數(shù)的極點(diǎn)對應(yīng)的是自由響應(yīng); 與激勵象函數(shù)的極點(diǎn)對應(yīng)的是強(qiáng)迫響應(yīng).Y(z)的極點(diǎn)都在z平面的單位圓內(nèi)故都是瞬態(tài)響應(yīng). 零輸入響應(yīng)和零狀態(tài)響應(yīng)可按定義參照求解.

上述求解過程可借助MATLAB的符號運(yùn)算編程實(shí)現(xiàn). 實(shí)現(xiàn)變換域法求解差分方程的m程序如下:

clc;clear;format compact;

syms z n %定義符號對象

% 輸入差分方程、初始狀態(tài)和激勵信號%

a=[1, -3/4, 1/8], b=[1, 1/3], %輸入差分方程系數(shù)向量

y0=[4, 12], x0=[0], %輸入初始狀態(tài), 長度分別比a、b短1, 長度為0時用[]

xn=(3/4)^n, %輸入激勵信號, 自動單邊處理, u(n)可用1^n表示

% 下面是變換域法求解差分方程的通用程序，極點(diǎn)為有理數(shù)時有解析式輸出 %

N=length(a)-1;M=length(b)-1;%計(jì)算長度

Az=poly2sym(a, 'z')/z^N;Bz=poly2sym(b, 'z')/z^M;%計(jì)算A(z)和B(z)

Hz=Bz/Az;disp('系統(tǒng)函數(shù)H(z):'), sys=filt(b, a), %計(jì)算并顯示系統(tǒng)函數(shù)

hn=iztrans(Hz);disp('單位樣值響應(yīng)h(n)='), pretty(hn), %計(jì)算并顯示單位樣值響應(yīng)

Hzp=roots(a);disp('系統(tǒng)極點(diǎn):');Hzp, %計(jì)算并顯示系統(tǒng)極點(diǎn)

Xz=ztrans(xn);disp('激勵象函數(shù)X(z)='), pretty(Xz), %激勵信號的單邊z變換

Y0z=0;%初始化Y0(z), 求Y0(z)注意系數(shù)標(biāo)號與變量下標(biāo)的關(guān)系

for k=1:N;

for l=-k:-1;Y0z = Y0z+a(k+1)*y0(-l)*z^(-k-l);end

end

disp('初始Y0(z)'), Y0z, %系統(tǒng)初始狀態(tài)的z變換

X0z=0;%初始化X0(z), 求X0(z)注意系數(shù)標(biāo)號與變量下標(biāo)的關(guān)系

for r=1:M;

for m=-r:-1;X0z = X0z+b(r+1)*x0(-m)*z^(-r-m);end

end

disp('初始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), % 計(jì)算并顯示全響應(yīng)

Yziz=-Y0z/Az;disp('零輸入象函數(shù)Yzi(z)='), pretty(Yziz), %零激勵響應(yīng)的z變換

yzin=iztrans(Yziz);disp('零輸入響應(yīng)yzi(n)='), pretty(yzin), % 計(jì)算并顯示零輸入響應(yīng)

Yzsz=(Bz*Xz+X0z)/Az;disp('零狀態(tài)象函數(shù)Yzs(z)='), pretty(Yzsz), %零狀態(tài)響應(yīng)的z變換

yzsn=iztrans(Yzsz);disp('零狀態(tài)響應(yīng)yzs(n)='), pretty(yzsn), % 計(jì)算并顯示零狀態(tài)響應(yīng)

該程序的運(yùn)行過程與手算過程對應(yīng), 顯示在命令窗的運(yùn)行結(jié)果與手算結(jié)果相同.

5 結(jié)束語

求解差分方程有多種方法. 迭代法計(jì)算原理簡單, 調(diào)用filter函數(shù)時要注意初始狀態(tài)的等效變換. 經(jīng)典法數(shù)學(xué)思路清晰, 調(diào)用rsolve函數(shù)時要注意邊界條件的選擇問題. 雙零法物理意義清楚, 可根據(jù)定義用經(jīng)典法計(jì)算. 變換域法簡便高效, 適合計(jì)算各類響應(yīng). 本文列舉的變換域法求解實(shí)例可作為課程教學(xué)的補(bǔ)充材料, 編寫的變換域法求解程序?qū)虒W(xué)科研也有一定的應(yīng)用價值.

[1] 鄭君里, 應(yīng)啟珩, 楊為理. 信號與系統(tǒng)[M]. 第2版. 北京: 高等教育出版社, 2000:16, 64

[2] 高西全, 丁玉美. 數(shù)字信號處理[M]. 第3版. 西安: 西安電子科技大學(xué)出版社, 2008:129

[3] 陳從顏, 翟軍勇. 信號與系統(tǒng)基礎(chǔ)[M]. 北京: 機(jī)械工業(yè)出版社, 2009:19

[4] 朱玲贊. 差分方程的邊界條件和離散系統(tǒng)的初始狀態(tài)[J]. 電氣電子教學(xué)學(xué)報(bào), 2001(03): 35～37

[5] 張志勇. 精通MATLAB6.5[M]. 北京: 北京航空航天大學(xué)出版社, 2003:229

[6] 史林杰. 離散時間系統(tǒng)邊界條件的確定準(zhǔn)則[J]. 電工教學(xué), 1997(02): 22～24

[7] 張登奇, 張 璇. 連續(xù)時間系統(tǒng)的s域分析及MATLAB實(shí)現(xiàn)[J]. 湖南理工學(xué)院學(xué)報(bào)(自然科學(xué)版), 2012(02): 26～29

Analysis of Difference Equation's Solution and Realization Based on MATLAB

ZHANG Deng-qi, PENG Shi-yu
(College of Information and Communication Engineering, Hunan Institute of Science and Technology, Yueyang 414006, China)

The difference equation is mathematical model to describe discrete-time systems, to solve the differential equation is an important part to analyze discrete-time systems, commonly used methods are: iterative method, classical time-domain method, double zero response method and the transform-domain method. This paper introduces the solving steps of these different methods according to the corresponding principles with examples and lists these corresponding MATLAB programs.

discrete-time system; difference equation; time-domain method; transform-domain method; MATLAB

TN911.72; O175.7

A

1672-5298(2014)03-0028-05

2014-06-25

張登奇(1968? ), 男, 湖南臨湘人, 碩士, 湖南理工學(xué)院信息與通信工程學(xué)院副教授. 主要研究方向: 信號與信息處理