差分方程的解法分析及MATLAB 实现（程序）
摘自：张登奇, 彭仕玉. 差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言
线性常系数差分方程是描述线性时不变离散时间系统的数学模型, 求解差分方程是分析离散时间系统的重要内容. 在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域
[1]法.
1 迭代法
y (n -1) +y (n -2) =x (n ) +x (n -1) , 激励信号为例1 已知离散系统的差分方程为y (n ) -x (n ) =() n u (n ) , 初始状态为y (-1) =4，y (-2) =12. 求系统响应. , y (1) =. 根据激励信号和初始状态, 手工依次迭代可算出y (0) =利用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),%迭代计算输出和后段等效初始条件
2 时域经典法
用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项, 根据自由项形
[3]式求特解；然后根据边界条件求完全解. 用时域经典法求解例1的基本步骤如下.
2α+=0, 可算出α1=, α2=. 高阶特征根可用MATLAB 的roots （1）求齐次解. 特征方程为α-) n +C 2() n , n ≥0 . 函数计算. 齐次解为y h (n ) =C 1() n u (n ) 代入差分方程右端得自由项为 （2）求方程的特解. 将x (n ) =( n =0⎧⎪1, n n -1() u (n ) +⋅() u (n -1) =⎨n ⋅() , n ≥1⎪⎩3n 当n ≥1时, 特解可设为y p (n ) =D () , 代入差分方程求得D =. 4
（3）利用边界条件求完全解. 当n =0时迭代求出y (0) =, 当n ≥1时, 完全解的形式为y (n ) =C 1() n +C 2() n +⋅() n , 选择求完全解系数的边界条件可参考文[4]选y (0), y (-1) . 根据边界δ(n ) 及其延, C =条件求得C 1=-2. 注意完全解的表达式只适于特解成立的n 取值范围, 其他点要用
迟表示, 如果其值符合表达式则可合并处理. 差分方程的完全解为
y (n ) =δ(n ) +[-⋅() n +⋅() n +⋅() n ]u (n -1) =[-⋅() n +⋅() n +⋅() n ]u (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))'),
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))'),
4 变换域法
设差分方程的一般形式为∑a k y (n -k ) =∑b r x (n -r ) .
k =0r =0N M
对差分方程两边取单边z 变换, 并利用z 变换的位移公式得
∑k =0N a k z [Y (z ) +∑y (l ) z ]=∑b r z [X (z ) +-k -l -r l =-k r =0-1M m =-r ∑-1x (m ) z -m ]
整理成A (z ) Y (z ) +Y 0(z ) =B (z ) X (z ) +X 0(z ) 形式有
A (z ) =a 0+a 1z -1+ +a N z -N , B (z ) =b 0+b 1z -1+ +b M z -M .
Y 0(s ) =∑∑a k y (l ) z
k =1l =-k N -1-k -l M -1 , X 0(s ) =∑∑b r x (m ) z -r -m . r =1m =-r
可以看出, 由差分方程可直接写出A (z ) 和B (z ) , 系统函数H (z ) =B (z ) /A (z ) , 将系统函数进行逆z 变换可得单位样值响应. 由差分方程的初始状态可算出Y 0(z ) , 由激励信号的初始状态可算出X 0(z ) , 将激励信号进行z 变换可得X (z ) , 求解z 域代数方程可得输出信号的象函数 B (z ) X (z ) +X 0(z ) -Y 0(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), ⋅⋅⋅, x (-M ) 算出Y 0(z ) 和X 0(z ) ；
（4）在z 域求解代数方程A (z ) Y (z ) +Y 0(z ) =B (z ) X (z ) +X 0(z ) 得输出象函数Y (z ) , Y (z ) 的逆变换即为全响应；
（5）分析响应象函数的极点来源及在z 平面中的位置, 确定自由响应与强迫响应, 或瞬态响应与稳态
响应；
（6）根据零输入响应和零状态响应的定义, 在z 域求解双零响应的象函数, 对双零响应的象函数进行逆z 变换, 得零输入响应和零状态响应.
用变换域法求解例1的基本过程如下. z -1+z -2, B (z ) =1+z -1. 系统函数的极点为, . 根据差分方程直接写出A (z ) =1-) 对激励信号进行z 变换得X (z ) =z /(z -. 激励象函数的极点为3/4. +z -1. 根据激励信号的初始状态算出X 0(z ) =0 . 根据差分方程的初始状态算出Y 0(z ) =-z 3-z 2+z ) /(z 3-z 2+z -) . 对z 域代数方程求解, 得全响应的象函数Y (z ) =(⋅() n +⋅() n +⋅() n ]u (n ) 进行逆z 变换得全响应为y (n ) =[-其中, 与系统函数的极点对应的是自由响应; 与激励象函数的极点对应的是强迫响应. 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),% 计算并显示零状态响应 该程序的运行过程与手算过程对应, 显示在命令窗的运行结果与手算结果相同.
差分方程的解法分析及MATLAB 实现（程序）
摘自：张登奇, 彭仕玉. 差分方程的解法分析及其MATLAB 实现[J]. 湖南理工学院学报.2014(03) 引言
线性常系数差分方程是描述线性时不变离散时间系统的数学模型, 求解差分方程是分析离散时间系统的重要内容. 在《信号与系统》课程中介绍的求解方法主要有迭代法、时域经典法、双零法和变换域
[1]法.
1 迭代法
y (n -1) +y (n -2) =x (n ) +x (n -1) , 激励信号为例1 已知离散系统的差分方程为y (n ) -x (n ) =() n u (n ) , 初始状态为y (-1) =4，y (-2) =12. 求系统响应. , y (1) =. 根据激励信号和初始状态, 手工依次迭代可算出y (0) =利用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),%迭代计算输出和后段等效初始条件
2 时域经典法
用时域经典法求解差分方程:先求齐次解；再将激励信号代入方程右端化简得自由项, 根据自由项形
[3]式求特解；然后根据边界条件求完全解. 用时域经典法求解例1的基本步骤如下.
2α+=0, 可算出α1=, α2=. 高阶特征根可用MATLAB 的roots （1）求齐次解. 特征方程为α-) n +C 2() n , n ≥0 . 函数计算. 齐次解为y h (n ) =C 1() n u (n ) 代入差分方程右端得自由项为 （2）求方程的特解. 将x (n ) =( n =0⎧⎪1, n n -1() u (n ) +⋅() u (n -1) =⎨n ⋅() , n ≥1⎪⎩3n 当n ≥1时, 特解可设为y p (n ) =D () , 代入差分方程求得D =. 4
（3）利用边界条件求完全解. 当n =0时迭代求出y (0) =, 当n ≥1时, 完全解的形式为y (n ) =C 1() n +C 2() n +⋅() n , 选择求完全解系数的边界条件可参考文[4]选y (0), y (-1) . 根据边界δ(n ) 及其延, C =条件求得C 1=-2. 注意完全解的表达式只适于特解成立的n 取值范围, 其他点要用
迟表示, 如果其值符合表达式则可合并处理. 差分方程的完全解为
y (n ) =δ(n ) +[-⋅() n +⋅() n +⋅() n ]u (n -1) =[-⋅() n +⋅() n +⋅() n ]u (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))'),
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))'),
4 变换域法
设差分方程的一般形式为∑a k y (n -k ) =∑b r x (n -r ) .
k =0r =0N M
对差分方程两边取单边z 变换, 并利用z 变换的位移公式得
∑k =0N a k z [Y (z ) +∑y (l ) z ]=∑b r z [X (z ) +-k -l -r l =-k r =0-1M m =-r ∑-1x (m ) z -m ]
整理成A (z ) Y (z ) +Y 0(z ) =B (z ) X (z ) +X 0(z ) 形式有
A (z ) =a 0+a 1z -1+ +a N z -N , B (z ) =b 0+b 1z -1+ +b M z -M .
Y 0(s ) =∑∑a k y (l ) z
k =1l =-k N -1-k -l M -1 , X 0(s ) =∑∑b r x (m ) z -r -m . r =1m =-r
可以看出, 由差分方程可直接写出A (z ) 和B (z ) , 系统函数H (z ) =B (z ) /A (z ) , 将系统函数进行逆z 变换可得单位样值响应. 由差分方程的初始状态可算出Y 0(z ) , 由激励信号的初始状态可算出X 0(z ) , 将激励信号进行z 变换可得X (z ) , 求解z 域代数方程可得输出信号的象函数 B (z ) X (z ) +X 0(z ) -Y 0(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), ⋅⋅⋅, x (-M ) 算出Y 0(z ) 和X 0(z ) ；
（4）在z 域求解代数方程A (z ) Y (z ) +Y 0(z ) =B (z ) X (z ) +X 0(z ) 得输出象函数Y (z ) , Y (z ) 的逆变换即为全响应；
（5）分析响应象函数的极点来源及在z 平面中的位置, 确定自由响应与强迫响应, 或瞬态响应与稳态
响应；
（6）根据零输入响应和零状态响应的定义, 在z 域求解双零响应的象函数, 对双零响应的象函数进行逆z 变换, 得零输入响应和零状态响应.
用变换域法求解例1的基本过程如下. z -1+z -2, B (z ) =1+z -1. 系统函数的极点为, . 根据差分方程直接写出A (z ) =1-) 对激励信号进行z 变换得X (z ) =z /(z -. 激励象函数的极点为3/4. +z -1. 根据激励信号的初始状态算出X 0(z ) =0 . 根据差分方程的初始状态算出Y 0(z ) =-z 3-z 2+z ) /(z 3-z 2+z -) . 对z 域代数方程求解, 得全响应的象函数Y (z ) =(⋅() n +⋅() n +⋅() n ]u (n ) 进行逆z 变换得全响应为y (n ) =[-其中, 与系统函数的极点对应的是自由响应; 与激励象函数的极点对应的是强迫响应. 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),% 计算并显示零状态响应 该程序的运行过程与手算过程对应, 显示在命令窗的运行结果与手算结果相同.