《差分方程的解法分析及MATLAB实现(程序).docx》由会员分享，可在线阅读，更多相关《差分方程的解法分析及MATLAB实现(程序).docx（3页珍藏版）》请在课桌文档上搜索。
1、差分方程的解法分析及MATLAB实现程序摘自：张登奇,彭仕玉.差分方程的解法分析及其MATLAB实现J.湖南理工学院学报.2023(03)引言线性常系数差分方程是描述线性时不变离散时间系统的数学模型,求解差分方程是分析离散时间系统的重要内容.在信号与系统课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域法Ul.1迭代法例1离散系统的差分方程为必)-分1)+打(-2)=1()+卜(-1),鼓励信号为X(n)=)(),初始状态为y(-l)=4,y(-2)=12.求系统响应.根据鼓励信号和初始状态,手工依次迭代可算出y(0)=9,ND=符.利用MATLAB中的filter函数实现迭代过程的
2、m程序如下：Clcjclearjformatcompact;a=1,-3/4,1/8,b=l,1/3,0,%输入差分方程系数向量,缺乏补0对齐n=0:IOixn=(3/4).n,%输入鼓励信号Zx=0,0,Zy=4,12,%输入初始状态zi=filtic(b,a,zy,zx),%计算等效初始条件yn,ZfrfilteNb,a,xn,zi),与迭代计算输出和后段等效初始条件2时域经典法用时域经典法求解差分方程:先求齐次解；再将鼓励信号代入方程右端化简得自由项,根据自由项形式求特解；然后根据边界条件求完全解.用时域经典法求解例1的根本步骤如下.求齐次解.特征方程为储Ta+1=0,可算出?二+，%高
3、阶特征根可用MATLAB的roots函数计算.齐次解为5)=C1()w+C2w,0.(2)求方程的特解.将文()=)()代入差分方程右端得自由项为当1时,特解可设为yp(n)=Oe),代入差分方程求得。%(3)利用边界条件求完全解.当/2=0时迭代求出义0)=*,当时，完全解的形式为M)=G(+)+g6)+3G)，选择求完全解系数的边界条件可参考文4选y(O),y(-l).根据边界条件求得G=-号，G=*.注意完全解的表达式只适于特解成立的取值范围,其他点要用/)及其延迟表示,如果其值符合表达式那么可合并处理.差分方程的完全解为MATLAB没有专用的差分方程求解函数,但可调用maple符号运算
4、工具箱中的FSOlve函数实现,格式为y=maple(rsolvc(cqus,inis,y(n),),其中：CqUS为差分方程表达式，inis为边界条件,y(n)为差分方程中的输出函数式.rsolve的其他格式可通过mhelprsolve命令了解.在MATLAB中用时域经典法求解例1中的全响应和单位样值响应的程序如下.clc;clear;formatcompact;yn-maple(,rsolve(y(n)-34*y(n-l)+l8*y(n-2)=(34)n+l3*(34)(n-l),y(0)=5/2,y(-1)=4,y(n),),hn=maple(,rsolve(y(n)-34*y(n-l)
5、+l8*y(n-2)=0,y(0)=1,y(l)=1312,y(n),),3双零法根据双零响应的定义,按时域经典法的求解步骤可分别求出零输入响应和零状态响应.理解了双零法的求解原理和步骤,实际计算可调用rsolve函数实现.yzi=maple(,rsolve(y(n)-34*y(nl)+l8*y(n-2)=0,y(-1)=4,y(2)=12,y(n),),yzs=mapleCrsolve(y(n)-34*y(n-l)+l8*y(n-2)=(3/4)-n+l3*(34)(n-l),y(0)=1,y(T)=0),y(n),4变换域法设差分方程的一般形式为S%M-幻=Whrxn-r).*-0r-4)
6、对差分方程两边取单边Z变换,并利用Z变换的位移公式得整理成A(Z)Y(Z)+Y0(Z)=B(Z)X(Z)+Xo(Z)形式有可以看出，由差分方程可直接写出A(Z)和B(Z),系统函数”(Z)=B(Z)A(z),将系统函数进行逆Z变换可得单位样值响应.由差分方程的初始状态可算出(z),由鼓励信号的初始状态可算出XO(Z),将鼓励信号进行Z变换可得X(Z),求解Z域代数方程可得输出信号的象函数对输出象函数进行逆z变换可得输出信号的原函数M).利用Z变换求解差分方程各响应的步骤可归纳如下：(1)根据差分方程直接写出A(z)、B(Z)和H(Z),H(Z)的逆变换即为单位样值响应；(2)根据鼓励信号算出X
7、(Z),如鼓励不是因果序列那么还要算出前M个初始状态值；(3)根据差分方程的初始状态N-I),y(-2),y(-N)和鼓励信号的初始状态X(T),M-2),x(-M)算出%(z)和X0(Z)；(4)在Z域求解代数方程4z)y(z)+%(z)=3(z)X(z)+X0(z)得输出象函数丫，Y(z)的逆变换即为全响应;(5)分析响应象函数的极点来源及在Z平面中的位置,确定自由响应与强迫响应,或瞬态响应与稳态响应；(6)根据零输入响应和零状态响应的定义,在Z域求解双零响应的象函数,对双零响应的象函数进行逆Z变换,得零输入响应和零状态响应.用变换域法求解例1的根本过程如下.根据差分方程直接写出A(Z)=
8、I-?7+卜-2回)=+9-系统函数的极点为*%对鼓励信号进行Z变换得X(Z)=Z/(z-6.鼓励象函数的极点为3/4.根据差分方程的初始状态算出E(Z)=+根据鼓励信号的初始状态算出XO(Z)=O.对Z域代数方程求解,得全响应的象函数y(Z)=(3-z2+Q)/(z3_*z2+Z-).进行逆Z变换得全响应为y(n)=-+jrt+3()其中，与系统函数的极点对应的是自由响应;与鼓励象函数的极点对应的是强迫响应.y(z)的极点都在Z平面的单位圆内故都是瞬态响应.零输入响应和零状态响应可按定义参照求解.上述求解过程可借助MATLAB的符号运算编程实现.实现变换域法求解差分方程的m程序如下：clc;
9、clear;formatcompact;symszn%定义符号对象%输入差分方程、初始状态和鼓励信号a=l,-34,l8,b=l,l3,%输入差分方程系数向量y=4,12,x0=0,%输入初始状态,长度分别比a、b短1,长度为0时用口xn=(34)-n,%输入鼓励信号，自动单边处理,u(n)可用ICn表示%下面是变换域法求解差分方程的通用程序，极点为有理数时有解析式输出%N=Iength(a)-1;M=Iength(b)-l;%计算长度Az=poly2sym(a,z)zN;BZ=PoIy2sym(b,z,)/zM;%计算A(Z)和B(z)HZ=BZ/Az;disp(系统函数H(z):*),sy
10、s=filt(b,a),%计算并显示系统函数hn=iztrans(Hz);disp(单位样值响应h(n)三,),pretty(hn),%计算并显示单位样值响应llzp=roots(a)jdisp(,系统极点：);Hzp,%计算并显示系统极点Xz=Ztrans(xn)disp(,鼓励象函数X(z)=),Pretty(XZ),%鼓励信号的单边z变换YoZ=0;%初始化Yo(Z),求Yo(Z)注意系数标号与变量下标的关系fork=l:N;forl=-k：-l;YOz=YOz+a(k+l)*yO(-l)*z(-k-l);endenddisp(初始YO(z),),Y0z,%系统初始状态的z变换XoZ=0
11、;%初始化Xo(Z),求XO(Z)注意系数标号与变量下标的关系forr=l:M;form=-r：-l;XOz=X0z+b(r+l)*x(-m)*z(-r-m);endenddisp(,初始XO(Z),XOz,%鼓励信号起始状态的z变换Yz=(Bz*Xz+X0z-Y0z)/Az;disp(,全响应的Z变换Y(z),pretty(simple(Yz),yn=iztrans(Yz);dispC全响应y(n)=,),Pretty(yn),%计算并显示全响应Yziz=-YOzzdisp(,零输入象函数YZi(Z)=),PretIy(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),%计算并显示零状态响应该程序的运行过程与手算过程对应,显示在命令窗的运行结果与手算结果相同.
