用matlab现线性常系数差分方程求解Word文档下载推荐.docx

1、与连续时间信号的微分及积分运算相对应，离散时间信号有差分及序列求和运算。设有序列f（k）,则称,f（k+2）,f（k+1）,f（k1）,f（k2）,为f（k）的移位序列。序列的差分可以分为前向差分和后向差分。一阶前向差分定义为 （3.11）一阶后向差分定义为 （3.12）式中和称为差分算子。由式（3.11）和式（3.12）可见，前向差分与后向差分的关系为 （3.13）二者仅移位不同，没有原则上的差别，因而它们的性质也相同。此处主要采用后向差分，并简称其为差分。由查分的定义，若有序列、和常数，则 （3.14）这表明差分运算具有线性性质。 二阶差分可定义为 （3.15）类似的，可定义三阶、四阶、n

2、阶差分。一般地，n阶差分 （3.16）式中 （3.17）为二项式系数序列f（k）的求和运算为 （3.18）差分方程是包含关于变量k的未知序列y（k）及其各阶差分的方程式，它的一般形式可写为 （3.19a）式中差分的最高阶为n阶，称为n阶差分方程。由式（3.16）可知，各阶差分均可写为y（k）及其各移位序列的线性组合，故上式常写为 （3.19b）通常所说的差分方程是指式（3.19b）形式的方程。若式（3.19b）中，y（k）及其各移位序列均为常数，就称其为常系数差分方程；如果某些系数是变量k的函数，就称其为变系数差分方程。描述LTI离散系统的是常系数线性差分方程。差分方程是具有递推关系的代数方程

3、，若一直初始条件和激励，利用迭代法渴求的差分方程的数值解。2. 差分方程的经典解一般而言，如果但输入单输出的LTI系统的激励f（k），其全响应为y（k），那么，描述该系统激励f（k）与响应y（k）之间关系的数学模型式n阶常系数线性差分方程，它可写为 （3.110a）式中、都是常数。上式可缩写为 （3.110b）与微分方程的经典解类似，上述差分方程的解由齐次解和特解两部分组成。齐次解用表示，特解用表示，即 （3.111）a.齐次解当式（3.110）中的f（k）及其各移位项均为零时，齐次方程 （3.112）的解称为齐次解。首先分析最简单的一阶差分方程。若一阶差分方程的齐次方程为 （3.113）它可

4、改写为y（k）与y（k1）之比等于a表明，序列y（k）是一个公比为a的等比级数，因此y（k）应有如下形式 （3.114）式中C式常数，有初始条件确定。对于n阶齐次差分方程，它的齐次解由形式为的序列组合而成，将代入到式（3.112），得由于C0，消去C；且0，以除上式，得（3.115）上式称为差分方程式（3.110）和式（3.112）的特征方程，它有n个根，称为差分方程的特征根。显然，形式为的序列都满足式（3.112），因而它们是式（3.110）方程的齐次解。依特征根取值的不同，差分方程齐次解的形式见表31,其中、等为待定常数表31 不同特征根所对应的齐次解特征根齐次解单实根重实根一对共轭复根重

5、共轭复跟b.特解特解的函数形式与激励的函数形式有关，表32列出了集中典型的激励f（k）所对应的特解。选定特解后代入原差分方程，求出其待定系数等，就得出方程的特解。表32 不同激励所对应的特解激励特解 所有特征根均不等于1时 当有重等于1时的特征根时 当不等于特征根时 当是特征单根时 当是重特征根时或所有特征根均不等于c.全解 式（3.110）的线性差分方程的全解是齐次解与特解之和。如果方程的特征根均为单根，则差分方程的全解为（3.116）如果特征根为重根，而其余nr个特征根为单根时，差分方程的全解为（3.117）式中各系数由初始条件确定。如果激励信号是在k=0时接入的，差分方程的解适合于k0。

6、对于n阶差分方程，用给定的n个初始条件y（0）,y（1）,y（n1）就可确定全部待定系数。如果差分方程的特解都是单根，则方程的全解为式（3.116），将给定的初始条件y（0）,y（1）,y（n1）分别代入到式（3.116），可得（3.118）由以上方程可求得全部待定系数。2.1零输入响应系统的激励为零，仅由系统的初始状态引起的响应，称为零输入响应，用表示。在零输入条件下，式（3.110）等号右端为零，化为齐次方程，即 （3.125）一般设定激励是在k=0时接入系统的，在k0时，激励尚未接入，故式（3.125）的几个初始状态满足 （3.126）式（3.126）中的y（1）,y（2）,y（n）为系

7、数的初始状态，由式（3.125）和式（3.126）可求得零输入响应。2.2零状态响应当系统的初始状态为零，仅由激励f（k）所产生的响应，称为零状态响应，用 表示。在零状态情况下，式（3.110）仍是非齐次方程，其初始状态为零，即零状态响应满足（3.130）的解。若其特征根均为单根，则其零状态响应为 （3.131）式中为待定常数，为特解。需要指出，零状态响应的初始状态为零，但其初始值不一定等于零。3.线性常系数差分方程3.1一个N 阶线性常系数差分方程可用下式表示： （1.4.1）或者 （1.4.2）式中，x（n）和y（n）分别是系统的输入序列和输出序列，ai和bi均为常系数，式中y（n-i）和

8、x（n-i）项只有一次幂，也没有相互交叉相乘项，故称为线性常系数差分方程。差分方程的阶数是用方程y（n-i）项中i的最大取值与最小取值之差确定的。在（1.4.2）式中，y（n-i）项i最大的取值N，i的最小取值为零，因此称为N阶差分方程。4 线性常系数差分方程的求解已知系统的输入序列，通过求解差分方程可以求出输出序列。求解差分方程的基本方法有以下三种：（1）经典解法。这种方法类似于模拟系统中求解微分方程的方法，它包括齐次解与特解，由边界条件求待定系数，上节已作简单介绍，这里不作介绍。（2）递推解法。这种方法简单，且适合用计算机求解，但只能得到数值解，对于阶次较高的线性常系数差分方程不容易得到封

9、闭式（公式）解答。（3）变换域方法。这种方法是将差分方程变换到z域进行求解，方法简便有效。当然还可以不直接求解差分方程，而是先由差分方程求出系统的单位脉冲响应，再与已知的输入序列进行卷积运算，得到系统输出。但是系统的单位脉冲响应如果不是预先知道，仍然需要求解差分方程，求其零状态响应解。（4）卷积法：由差分方程求出系统的h（n），再与已知的x（n） 进行卷积，得到y（n）。观察（1.4.1）式，求n时刻的输出，要知道n时刻以及n时刻以前的输入序列值，还要知道n时刻以前的N个输出序列值。因此求解差分方程在给定输入序列的条件下，还需要确定N个初始条件。如果求n0时刻以后的输出，n0时刻以前N个输出值

10、y（n0-1）、y（n0-2）、y（n0-N）就构成了初始条件。（1.4.1）式表明，已知输入序列和N个初始条件，则可以求出n时刻的输出；如果将该公式中的n用n+1代替，可以求出n+1时刻的输出，因此（1.4.1）式表示的差分方程本身就是一个适合递推法求解的方程。三设计过程1.用MATLAB求解差分方程 MATLAB信号处理工具箱提供的filter函数实现线性常系数差分方程的递推求解，调用格式如下： yn=filter（B,A.xn） 计算系统对输入信号向量xn的零状态响应输出信号向量yn，yn与xn长度相等，其中，B和A是（1.4.2）式所给差分方程的系数向量，即 B=b0,b1,bM, A

11、=a0,a1,aN其中a0=1，如果a01，则filter用a0对系数向量B和A归一化。 yn=filter（B,A.xn,xi） 计算系统对输入信号向量xn的全响应输出信号yn。所谓全响应，就是由初始状态引起的零输入响应和由输入信号xn引起的零状态响应之和。其中，xi是等效初始条件的输入序列，所以xi是由初始条件确定的。MATLAB信号处理工具箱提供的filtic就是由初始条件计算xi的函数，其调用格式如下：xi=filtic（B,A,ys,xs）其中，ys和xs是初始条件向量：ys=y（-1）,y（-2）,y（-3）,y（-N）,xs=x（-1）,x（-2）,x（-3）,x（-M）。如果x

12、n是因果序列，则xs=0.调用时可缺省xs。 例1.4.1的MATLAB求解程序ep141.m如下：%1.4.1.m:调用MATLAB解差分方程y（n）-0.8y（n-1）=x（n）a=0.8；ys=1; %设差分方程系数a=0.8，初始状态：y（-1）=1xn=1,zeros（1,30）; %x（n）=单位脉冲序列，长度N=31B=1;A=1,-0.8; %差分方程系数xi=filtic（B,A,ys）; %由初始条件计算等效初始条件的输入序列xiyn=filter（B,A,xn,xi）; %调用fiter解差分方程，求系统输出信号y（n）n=0:length（yn）-1;stem（n,yn

13、,.）title（时域波形图）;xlabel（nylabel（y（n）程序中取查分方程系数a=0.8时，得到系统输出y（n）如图1.4.1（a）所示，与例1.4.1的解析递推结果完全相同。如果令初始条件y（-1）=0（仅修改程序中ys=0），则得到系统输出y（n）=h（n），如图1.4.1（b）所示。（a）（b）图（a）为a=0.8，y（-1）=1时，系统输出时域波形图，图（b）为a=0.8，y（-1）=0时，系统输出时域波形图。四设计代码及结果MATLAB源程序源程序如下%1.m: %初始状态： %x（n）=单位脉冲序列，长度N=31 %由初始条件计算等效初始条件的输入序列xi %调用filter解差分方程，求系统输出信号y（n） %n的取值范围） %画出时域波形图） %x轴、y轴分别代表n，x（n）n=-5:5;xn=0.5.n; %xn=0.5.nstem（n,xn,fill）,grid on %画出时域波形图）,ylabel（x（n）, title（时域波形图 %n的取值范围a=0.5;