张登奇, 彭仕玉

(湖南理工学院 信息与通信工程学院, 湖南 岳阳 414006)

差分方程的解法分析及其MATLAB实现

张登奇, 彭仕玉

(湖南理工学院 信息与通信工程学院, 湖南 岳阳 414006)

差分方程是描述离散时间系统的数学模型, 求解差分方程是分析离散时间系统的重要内容, 常用的求解方法有迭代法、时域经典法、双零法和变换域法. 文章根据各种方法的求解原理, 分别介绍了不同方法的求解步骤, 结合实例列出了这些方法的求解过程及MATLAB实现程序.

离散时间系统; 差分方程; 时域分析法; 变换域分析法; MATLAB

引言

线性常系数差分方程是描述线性时不变离散时间系统的数学模型, 求解差分方程是分析离散时间系统的重要内容. 在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法[1]. 迭代法可手工逐次代入求解, 也可编程用计算机求解, 该方法原理简单, 缺点是只能得到数值解.时域经典法先求齐次解和特解, 再用边界条件确定待定系数得完全解, 该方法数学过程清晰, 但求解过程麻烦. 双零法分别求零输入响应和零状态响应, 再通过叠加得到全响应. 该方法物理意义清晰, 但求解过程依然麻烦. 变换域法利用z变换将差分方程变换成代数方程求解, 该方法简便高效, 是求解差分方程的重要方法. 本文根据不同方法的求解原理, 分别介绍各种方法的求解步骤, 结合实例列出这些方法的求解过程和MATLAB实现程序.

1 迭代法

差分方程本身就是一个递推方程, 根据初始状态和激励信号依次迭代就可算出输出序列. 迭代法是解差分方程的基础方法, 如果所需输出序列个数较少（如计算边界条件）可手工直算, 如需计算大量输出可利用计算机编程实现.现结合实例介绍迭代法的计算过程.

例1已知离散系统的差分方程为激励信号为初始状态为y(−1)=4,y(−2)=12.求系统响应.

利用MATLAB中的filter函数实现迭代过程的m程序如下:

clc;clear;format compact;

a=[1, -3/4, 1/8], b=[1, 1/3, 0], %输入差分方程系数向量, 不足补0对齐

n=0:10;xn=(3/4).^n, %输入激励信号

zx=[0, 0], zy=[4, 12], %输入初始状态

zi=filtic(b, a, zy, zx), %计算等效初始条件

[yn, zf]=filter(b, a, xn, zi), %迭代计算输出和后段等效初始条件

MATLAB提供的filter函数是一个内建函数, 用type命令看不到程序代码.为了理解迭代思想, 下面根据图1所示的直接Ⅰ型结构[2], 重写实现迭代法的m程序.

clc;clear;format compact;

a=[1, -3/4, 1/8], b=[1, 1/3, 0], %输入差分方程系数向量, 不足补0对齐

n=0:10;x=(3/4).^n, %输入激励信号

zx=[0, 0], zy=[4, 12], %输入初始状态

% 下面是按直接Ⅰ型结构迭代的通用程序 %

N=length(a)-1, %计算数据存储长度

L=length(x), %计算激励信号长度

y=zeros(1, L);%输出信号初始化

for i=1:L; %逐个计算输出信号

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

zz=sum(z);%计算输出中的过去分量

y(i)=b(1)*x(i)+zz;% 计算当前输出y(n)

for n=N:-1:2, zx(n)=zx(n-1);zy(n)=zy(n-1);end%过去数据下移

zx(1)=x(i);zy(1)=y(i);%当前的激励和输出变为过去, 以便算下一个输出

end

%理解filter函数中zf参数的意义

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

for k=1:N; %逐个计算zf参数

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

zf(k)=sum(z);% 计算第k个zf

for k=N:-1:2, zx(k)=zx(k-1);zy(k)=zy(k-1);end;%过去数据下移

zx(1)=0;zy(1)=0;% 没有当前的激励和输出变为过去, 算下一个zf

end

y, zf, %显示输出和zf参数

该程序采用直接Ⅰ型结构实现迭代计算, 与filter函数计算的结果相同, zf参数可用于激励信号分段处理时，调用filter函数计算下一段输出所需的等效初始条件.

图1 直接Ⅰ型网络结构图

2 时域经典法

用时域经典法求解差分方程与高等数学中求解微分方程的过程类似: 先求齐次解; 再将激励信号代入方程右端化简得自由项, 根据自由项形式求特解; 然后根据边界条件求完全解[3]. 用时域经典法求解例1的基本步骤如下.

当n≥1时, 特解可设为代入差分方程求得

(3) 利用边界条件求完全解.当n=0时迭代求出当n≥1时, 完全解的形式为选择求完全解系数的边界条件可参考文[4]选y(0),y(−1). 根据边界条件求得注意完全解的表达式只适于特解成立的n取值范围, 其他点要用δ(n)及其延迟表示, 如果其值符合表达式则可合并处理.差分方程的完全解为

MATLAB没有专用的差分方程求解函数, 但可调用maple符号运算工具箱中的rsolve函数实现[5], 格式为y=maple('rsolve({equs, inis}, y(n))'), 其中equs为差分方程表达式, inis为边界条件,y(n)为差分方程中的输出函数式. rsolve的其他格式可通过mhelp rsolve命令了解. 在MATLAB中用时域经典法求解例1中的全响应和单位样值响应的程序如下:

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))'),

如同时域经典法一样, 应用rsolve函数时要特别注意边界条件的选择问题，且边界条件要连续取点.本例求单位样值响应时的边界条件要选用y(0)和y(1)[6], 在命令窗中显示的结果也要分析n的取值范围.

3 双零法

双零法是将完全解分解成物理意义明显的零输入响应和零状态响应分别计算. 零输入响应是激励为零, 由系统的初始状态所产生的响应; 零输入响应要求差分方程右端为零, 故特解为零; 完全解为齐次解形式, 系数可直接由初始状态确定. 零状态响应是初始状态为零, 由激励信号所产生的响应. 计算零状态响应可用时域经典法, 也可用卷积法. 根据双零响应的定义, 按时域经典法的求解步骤可分别求出零输入响应和零状态响应. 理解了双零法的求解原理和步骤, 实际计算可调用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)^n+1/3*(3/4)^(n-1), y(0)=1, y(-1)=0}, y(n))'),

边界条件和n的取值范围这里不再赘述, 各响应的计算结果在命令窗都有显示, 零输入响应和零状态响应之和与前面计算的全响应相等.

4 变换域法

变换域法是以z变换为数学工具, 先将时域差分方程变换成z域代数方程, 再在z域求解响应的象函数, 最后将响应象函数逆变换成时域原函数, 是一种间接求解差分方程的计算方法.z变换的线性和位移性是变换域法求解差分方程的重要性质, 这里只列出双边序列的单边z变换位移公式[1]:

对差分方程两边取单边z变换, 并利用z变换的位移公式, 得

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

可以看出, 由差分方程可直接写出A(z)和B(z), 系统函数将系统函数进行逆z变换可得单位样值响应. 由差分方程的初始状态可算出Y0(z), 由激励信号的初始状态可算出X0(z), 将激励信号进行z变换可得X(z), 求解z域代数方程可得输出信号的象函数

对输出象函数进行逆z变换可得输出信号的原函数y(n). 为了计算简单和便于对结果进行对比分析,例1列举的实例中, 系统函数的极点均为单根有理数极点, 且系统函数H(z)的零极点与激励象函数X(z)的零极点均不相同. 参照s域求解微分方程的方法[7], 利用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),⋅⋅⋅,x(−M)算出Y0(z)和X0(z);

(4) 在z域求解代数方程A(z)Y(z)+Y0(z)=B(z)X(z)+X0(z)得输出象函数Y(z),Y(z)的逆变换即为全响应;

(5) 分析响应象函数的极点来源及在z平面中的位置, 确定自由响应与强迫响应, 或瞬态响应与稳态响应;

(6) 根据零输入响应和零状态响应的定义, 在z域求解双零响应的象函数, 对双零响应的象函数进行逆z变换, 得零输入响应和零状态响应.

用变换域法求解例1的基本过程如下.

对激励信号进行z变换得激励象函数的极点为

对z域代数方程求解, 得全响应的象函数

进行逆z变换得全响应为

其中, 与系统函数的极点对应的是自由响应; 与激励象函数的极点对应的是强迫响应.Y(z)的极点都在z平面的单位圆内故都是瞬态响应. 零输入响应和零状态响应可按定义参照求解.

上述求解过程可借助MATLAB的符号运算编程实现. 实现变换域法求解差分方程的m程序如下:

clc;clear;format compact;

syms z n %定义符号对象

% 输入差分方程、初始状态和激励信号%

a=[1, -3/4, 1/8], b=[1, 1/3], %输入差分方程系数向量

y0=[4, 12], x0=[0], %输入初始状态, 长度分别比a、b短1, 长度为0时用[]

xn=(3/4)^n, %输入激励信号, 自动单边处理, u(n)可用1^n表示

% 下面是变换域法求解差分方程的通用程序，极点为有理数时有解析式输出 %

N=length(a)-1;M=length(b)-1;%计算长度

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

Hz=Bz/Az;disp('系统函数H(z):'), sys=filt(b, a), %计算并显示系统函数

hn=iztrans(Hz);disp('单位样值响应h(n)='), pretty(hn), %计算并显示单位样值响应

Hzp=roots(a);disp('系统极点:');Hzp, %计算并显示系统极点

Xz=ztrans(xn);disp('激励象函数X(z)='), pretty(Xz), %激励信号的单边z变换

Y0z=0;%初始化Y0(z), 求Y0(z)注意系数标号与变量下标的关系

for k=1:N;

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

end

disp('初始Y0(z)'), Y0z, %系统初始状态的z变换

X0z=0;%初始化X0(z), 求X0(z)注意系数标号与变量下标的关系

for r=1:M;

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

end

disp('初始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)='), pretty(yzin), % 计算并显示零输入响应

Yzsz=(Bz*Xz+X0z)/Az;disp('零状态象函数Yzs(z)='), pretty(Yzsz), %零状态响应的z变换

yzsn=iztrans(Yzsz);disp('零状态响应yzs(n)='), pretty(yzsn), % 计算并显示零状态响应

该程序的运行过程与手算过程对应, 显示在命令窗的运行结果与手算结果相同.

5 结束语

求解差分方程有多种方法. 迭代法计算原理简单, 调用filter函数时要注意初始状态的等效变换. 经典法数学思路清晰, 调用rsolve函数时要注意边界条件的选择问题. 双零法物理意义清楚, 可根据定义用经典法计算. 变换域法简便高效, 适合计算各类响应. 本文列举的变换域法求解实例可作为课程教学的补充材料, 编写的变换域法求解程序对教学科研也有一定的应用价值.

[1] 郑君里, 应启珩, 杨为理. 信号与系统[M]. 第2版. 北京: 高等教育出版社, 2000:16, 64

[2] 高西全, 丁玉美. 数字信号处理[M]. 第3版. 西安: 西安电子科技大学出版社, 2008:129

[3] 陈从颜, 翟军勇. 信号与系统基础[M]. 北京: 机械工业出版社, 2009:19

[4] 朱玲赞. 差分方程的边界条件和离散系统的初始状态[J]. 电气电子教学学报, 2001(03): 35～37

[5] 张志勇. 精通MATLAB6.5[M]. 北京: 北京航空航天大学出版社, 2003:229

[6] 史林杰. 离散时间系统边界条件的确定准则[J]. 电工教学, 1997(02): 22～24

[7] 张登奇, 张 璇. 连续时间系统的s域分析及MATLAB实现[J]. 湖南理工学院学报(自然科学版), 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− ), 男, 湖南临湘人, 硕士, 湖南理工学院信息与通信工程学院副教授. 主要研究方向: 信号与信息处理