Chapter 10 基于有限体积法的控制方程离散
10.1 控制方程
稳态气相通用控制方程如下：
()
()()+div ρU div Γgrad S t ρφφφ??=+
()rv u w U U z r r r θU=ui+vj+wk
div ???=??=++
???
()()()() z r r u r v w t z r r r r S
z r r r ?ρφ???
ρφρφρφ????θ
??φ??φ??φ?????θ?θ+++=Γ??????
Γ+Γ++ ? ? ???????
上式中φ表示气相场通用变量，在三维柱坐标下可以是轴向速度u 、径向速度v 、
周向速度w 以及焓h ，φ为1时代表连续方程，Γ表示输运系数，S 为气相源项，U
为速度矢量。密度由完全气体状态方程ρRT p =确定，温度由焓温关系式(h)T f =确定。每个方程中的φ、Γ、S 如表10.1所示。
表10.1 气相方程组
§10.2 方程的离散化
§10.2.1 有限体积法
有限差分法：原理比较简单，数学推演和程序编制的工作量较小，但其灵活性较差；
有限单元法：数学推演和程序编制工作要比有限差分法复杂，因而比较费时，但其灵活性较好；
有限体积法：作为有限差分法和有限单元法的中间产物，继承了两者的优点，克服了各自的缺点，因而在流场计算中占有重要的地位。
有限体积法的基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解微分方程对每一控制体积积分，得出离散方程，其中的未知数是网格点上的因变量φ。有限体积法在对控制体积积分时，必须假定φ值在网格点之间的分布。
由有限体积法导出的离散方程具有明确的物理意义，即因变量φ在有限大小的控制体积内的守恒原理，如同推导微分方程时因变量在无限小的控制体积中的守恒原理。有限体积法即使在粗网格条件下，也精确地满足
积分守恒。
图10-1 三维控制体积
()()()
[]
z r r z z r z r ,,,,,z z
r r
,,,,,i+j+k
r z r r r θ i+j+k S i+S j+S k S +S +S nb w e s n b t
nb
nb w e s n b t
nb
div div dV ndS
θθθθθθ==Φ=ΦΦΦ?Φ?Φ?ΦΦ=++
???Φ=Φ???=ΦΦΦ???=
ΦΦΦ??∑∑
()()
[]
[]()[][]z
z z r ,,,,,z
z z
,,,,,z
z z ,,,,,r r r
r r ,,,,,,,,,, i+0j+0k S i+S j+S k z S z
S z r S r r r r S r θnb w e s n b t nb
nb w e s n b t
nb
nb w e s n b t
nb
nb w e s n b t nb nb w e s n b t
nb
dV V V V V θθ
θθ=====?Φ??=Φ?????Φ=
Φ=
???Φ=Φ??Φ?ΦΦ=Φ?=+???Φ=Φ??∑?∑∑∑∑
有限体积法的四条基本原则是： 原则1：控制体积交界面的一致性。
当一个表面为相邻的两个控制体积所共有时，在这两个控制体积的离
散方程中，通过交界面的通量表达式必须相同。 原则2：正系数。
中心节点系数a P 和相邻节点系数a nb 必须恒为正值。 原则3：源项的负坡线性化。
源项S 常依赖于因变量φ，为便于计算，通常将源项线性化，即将其表达为S S S C P =+φ，此时要求S P 必须小于或等于零。从数值计算的角度看来，可避免出现计算不稳定或结果不合理。
如径向动量方程源项中含有v 2μ
2μv, 0, r
r
C P C P S S S S --=+==
周向动量方程源项中含有ρvw w, 0, ρv C P C P S S S S -=+==- 原则4：相邻节点系数和。
当S P 为零时，中心节点系数a P 应等于各相邻节点系数之和，即满足
a a P n
b =∑。
§10.2.2方程的离散化 将通用控制方程改写为：
()?ρφ???ρφ?φ???ρφ?φ???θρφ?φ?θ t z u z r r r v r r w r S +-?? ???+-?? ?????
??
??
+
-?? ??
?=ΓΓΓ11
（1） 首先按常用方法将源项作线性化处理，即将源项表达为： S S S c p p =+φ （2） 其次对非定常项
()
?ρφ? t
进行离散化，令 ()?ρφ?ρφρφ t t
p p p p
=-00
? （3） 此处采用全隐式，假定“新”数值ρφP P ,控制着整个控制体积。上标0
表示时间步长开始时的“老”数值，不带上标的量均表示时间步长结束时的“新”数值。
接下来将方程（1）对图10-1所示控制体积积分，得：
()
ρφρφφp p p p
e w n s t b c p p t V J J J J J J S S V -??
??
?+-+-+-=+00
? （4） 上式中V 是控制体积的容积，可表达为
()V r r z r n s =+05. ???θ （5）
r r n s ,代表控制体积的上、下半径，?z 、?r 、?θ分别为控制体积在三个坐
标方向的长度。
我们以e ,w ,n ,s ,t ,b 来代表控制体积的6个交界面，其位置如图10-1所示。则：
J u z A J u z A J v A J v A J w r A J w r A w w w e e
e n n n s s
s t t t b b
b =-?? ???=-?? ???=-??
???=-?
? ???=-?? ???=-?
? ???ρφ?φ?ρφ?φ?ρφ?φ?ρφ?φ?ρφ?φ?θρφ?φ?θΓΓΓΓΓΓ
r r
（6）
式中A w 、A e 、A n 、A s 、A t 、A b 分别为控制体积的左、右、上、下、前、后表面积，且有
()A A r r A A r z A r z A r z w e s t b n n s s =+=== =0.5r = ;
, ;
n ????????θθθ;
（7）
将连续方程
()()()?ρ??ρ??ρ??ρ?θ t z r +++=u v w r 0 对控制体积积分，得
()
ρρp
p e w n s t b V t
F F F F F F -+-+-+-=00? （8）
式中F w , F e , F n ,F s ,F t ,F b 分别为控制体积6个表面上的质量流率，即有
()()()()()()F u A F u A F v A F v A F w A F w A e e e w w w n n n s s s t t t b b b ======ρρρρρρ , ;
, ; , ;
（9）
用（4）式减去（8）式与φP 的乘积，得到
()()()()()()()()φφρφφφφφφφp
p
p e
e
p
w
w
p
n
n
p
s
s
p
t
t
p
b
b
p
c
p
p
V
t J F J F J F J F J F J F S S V
-+---+---+---=+00? （10）
图10-2 两节点间的总通量J e
接下来有必要讨论一下通用离散格式的问题。考虑如图10-2所示在Z 向相距为δz 的两个节点P 和E ，令交界面上由对流通量ρφe e e e u A 和扩散通量
-Γe
e
e d dz
A φ组成的总通量为J e ，即 J u d dz A e e
e =-?? ?
??ρφφΓ
（11） 若令J J u e e z e e e z e
*=
=δρδ
ΓΓ , P ,e 则有 ()J P d d z A e
e e e z e *=-
?? ?
?
?φφ （12） 在上式中交界面上的φ值φe 应为φP 和φE 的加权平均，梯度
()
d d z e
z φ应为()φφE P -与某因子的乘积，因此不妨假设
()[]()J A P e e e P E E P *=+---αφαφβφφ1 （13） 式中α和β为依赖于P e 的无因次乘子，该式又可进一步改写为
J A B A e e P E *=-φφ （14）
式中无因次系数A 和B 是P e 的函数；A 与交界面前方的节点E 相联系，B 与交界面后方的节点P e 相联系。
当φφP E =时，扩散通量为零，总通量只是对流通量ρφe e e e u A ，即： *e e P E e P e E J A B A P P φφφφ=-== （15） 故有：
B A P e =+ （16）
当改换坐标轴方向时，P e 换为-P e ，A 与B 互相替换，故函数()e A P 与
()e B P 之间必有对称性关系：
()()A P B P e e -= （17） 或
()()B P A P e e -= （18） 由（16）至（18）式可以看出，对于任意P e 值，有
()()()()[]()()[]
A P
B P P A P P A P P B P A P P e e e e e e e e e e =-=--=+-=+ , 0,0
（19）
表10-2 各类离散格式的函数()A P
将（16）式代入（14）式，可得出
()()J A A P A A P e e e P E P E e P *=+-=-+φφφφφ （20） 即
()J A P A e e e P P E *-=-φφφ （21） 将J J u e e z e e e z e
*=
=δρδ
ΓΓ , P ,e 代入（21）式可得 ()()J F A A P e e P e e
z
e P E -=-φδφφΓ （22）
令D A e e e
z
=
Γδ，派克里特数P F D i i i = , i =e,w,n,s,t,b ，即：
()()[]a D A P D A P F E e e e e e ==+- , 0，故有：
()J F a e e P E P E -=-φφφ （23） 用类似的方法可推导出如下表达式：
()()
()()()
J F a J F a J F a J F a J F a w w P W W P n n P N P N s s P S S P t t P T P T b b P B B P -=--=--=--=--=-φφφφφφφφφφφφφφφ （24）
将（23）与（24）式代入（10）并加以整理，即得到三维通用控制方程（1）式的离散方程如下：
a a a a a a a
b P P E E W W N N S S T T B B φφφφφφφ=++++++ （25） 式中
()[]a D A P F E e e e =+- , 0 （26a ） ()[]a D A P F W w w w =+ , 0 （26b ） ()[]a D A P F N n n n =+- , 0 （26c ）
()[]a D A P F S s s s =+ , 0 （26d ） ()[]a D A P F T t t t =+- , 0 （26e ） ()[]a D A P F B b b b =+ , 0 （26f ） a V
P
P 00=
ρ? t
（26g ）
b S V a C P P =+00
φ （26h ） a a a a a a a a S V P E W N S T B P P =++++++-0 （26i ）
容积V 的计算式见(5 )，对流强度F i (i =e,w,n,s,t,b)的表达式参见(9 )，扩散率则定义如下： D A D A D A r e w z e w n s r n s t b
P t b
,,,,,,,,=??
???=?? ?
??=?? ???ΓΓΓδδδθ （27） 此外，中心网格点P 处的径向坐标值()r r r P n s =+05.。各类离散格式的函数
()A P 陈列在表10-2中。推荐采用如下的幂函数格式：
()()[]
A P P =-01015
,. （28）
§10.2.3同位网格下的压力校正方程
以z 向二维动量方程为例(参照上图)，其压力梯度项计算如下：
()()()
()() sn we z e w z n s sn we r e w r n s p
dV S p p S p p z rp dV S p p S p p r r ?=
-+-??=
-+-???
(29)
sn ew z z S S 、和sn ew r r S S 、分别为控制体积中心面sn 、we 在z 、r 向的投影面积。
将压力梯度项从源项中分离出来后z 向动量方程可改写成以下形式：
()()()()sn we nb nb u z z P e w n s P P P sn we
nb nb v r r P e w
n s P P P A u b S S u p p p p A A A A v b S S v p p p p A A A ??
??∑+=---- ? ?????
??
??∑+=---- ? ?????
(30)
上式中下标nb 代表E, W, S, N 等邻点。在交错网格安排下待解的速度变量存储在控制体积界面上，可直接用于界面流量计算，而在同位网格下所有变量都存储在控制体积中心，控制体积界面上的速度须由其相邻中心节点的速度值插值而来。通常，简单的插值方法容易导致振荡，而高阶的插值方法则将使压力校正方程变得非常复杂，且在大多数情况下数值不稳定。为避免出现振荡，同时仍能使压力校正方程保持简洁的形式，可采用一种将界面上的速度与跨过界面的压力梯度相联系的方法，以界面e 为例，该处速度分量e e u v 、可由E P u u ,表达式内插计算如下：
()()()()we sn nb nb u z n s z e E P P
P e e we sn nb nb v r
n s r e E P P
P e
e A u b S p p S u p p A A A v b S p p S v p p A A ??∑+--??
=-- ? ???????
∑+--??=--
? ????
? (31)
带上划线的部分由其邻点表达式中对应项内插而得。 同理可得界面n 处速度分量n n u v 、为：
()()()()sn we nb nb u z e w z n N P P P n n sn we nb nb v r e w r n N P P P n
n A u b S p p S u p p A A A v b S p p S v p p A A ??∑+--??
=-- ? ???????∑+--??=-- ? ????? (32)
令*'*'
,e e e e e e
u u u v v v =+=+，带上标*者代表近似值，带上标’者代表修正值，由公式(31)可得速度-压力校正关系：
()()()()()()()()
''
''''
'
'
'
'
11e P
E P E e P
E
P
E
sn BC z P P z e e sn
BC r
P P r
e
e
u S A p p A S p p v S
A p
p A S p
p =-≈-=-≈- (33)
上式中P '
代表正确的压力场与猜测压力场P *
之差，则有
P P P =+*' (34)
同理可得界面n 处速度-压力校正关系：
()()()()()()()()
''
'''''
''
'11n P
N P N n P
N
P
N
we DC z P P z n n we DC r P P r
n
n
u S A p
p A S p p v S A p
p A S p
p =-≈-=-≈- (35)
界面e 处对流通量为：
()()()()()()()
*'*'**
''**''**''2
()() 1 1P E P E BC BC BC BC
e e e z e r e e e z e e r BC BC BC
BC BC BC e e e z e r e e P z
z r r e e e P e e F u S v S u u S v v S F u S v S F A p p S
S S S F A p p S ρρρρρ??=+=+++??=++=+-?+?=+- (36)
同理可得界面n 处对流通量为：
()()
**
''21P N n n n P n n F F A p p S ρ=+-
(37)
最终可由连续方程
()0
0p p e w n s V
F F F F t
ρρ-+-+-=?导出如下压力校正方程：
'''
''P P E E W W N N S P a P a P a P a P a P b =++++
(38)
式中
()()()()2
22
2
1111E e P e e W w P w
w N n P n
n S s P s
s P E W N S
a A S a A S
a A S a A S
a a a a a ρρρρ=====+++ (38a)
()0
P
P w e s n V b F F F F t
ρρ-=
+-+-?
(38b)
上式中上标*被省略，交界面上的密度值e n ρρ、等应采用插值方法求出。
§10.2.4 压力校正方程边界条件
将(36) 、(37)改写如下：
()()
*''*'
'P E P
e e E n n
N
N
F F a p p F F a p p =+-=+- (39)
如果E 、N 边界为壁面、对称轴，由于对流通量恒为0，故有对应的0E N a a =、；
对于超音速出口边界，可假定''''
P E P N p p or p p ==，也有0E N a a =、；
对于亚音速出口边界，出口压强等于外界压强，''=0 0E N p or p =，按(38a)计算
E N P a a a 、、，上划线部分按中心点P 计算，然后取0E N a a =、；
对于左端流量入口边界有0W a =。
对于左端给定压力入口边界按(38a)计算W P a a 、，上划线部分按中心点P 计算，然后取0W a =；
对于左端给定总温总压入口边界，左边界压力等于中心点P 压力，应取0W a =； §10.3 同位网格下的SIMPLE 计算程式
SIMPLE 的含义为“求解与压力耦连的方程组的半隐式方法”(Semi －Implicit Method for Pressure －Linked Equations )，由Patankar 和Spalding 首先提出。其执行顺序如下：
① 根据经验给出压力场的初始猜测值P *，由(31)式计算界面速度，由(36)、(37) 式计算对流通量；
② 求解动量方程(30)式，得出**,u v ； ③ 求解压力校正方程(38)式，得出P '； ④ 由(34)式算出压力P ；
⑤ 由(39)式校正界面对流通量，仿照(35)式用如下的速度校正公式(40)求出 , u v ；
()
()
*
''
*''
sn z e w P we
r n s P S u u p p A S v v p p A ??=-- ???
??=-- ???
(40)
界面上的压力校正值可由相邻中心点上的对应值线性插值得到。
⑥ 如果有某些φ的数值通过流体的性质或源项对流场施加影响，如k -ε方程通过μt 与动量方程耦连，则应求解此类φ的离散方程。如果φ的数值对流场不施加影响，则应在得出流场收敛解后再求解这些φ的离散方程；
⑦ 把经过校正的压力P 作为新的猜测压力P *，返回第②步重复整个过程直到得出收敛解。
§10.4 松弛迭代法
由于压力校正方程本身并不精确，如果略去的项过多，有可能导致迭代过程发散，特别是在猜测压力场远离正确压力场时更是如此。经验表明，为使压力校正方程稳定地趋于收敛，有必要采取欠松弛迭代，即在求解动量方程(30)时将P A 改写为P A α，并
在方程右端加上()01P A αφα-，φ0代表前一次迭代的φ值，α就是松弛系数。表述如下：
()()()()()()0
11sn we P
P nb nb u z e w z n s P P
sn we P P nb nb v r e w r n s P P A u A u b S p p S p p A u A v A v b S p p S p p A v αα
αααα
-=∑+----+
-=∑+----+ (30’) 当01<<α时为欠松弛，其作用是减缓未知量的变化，α>1时则加速未知量的变化，为超松弛。
此外，还应对压力校正方程进行欠松弛，即以下式取代(34)式： P P P P =+*'α (34’)
原则上松弛系数α和αP 可任意取值，然而α过小时收敛极慢，甚至于出现假收敛，
α过大时方程仍可能发散。大部分人推荐αα==05
08.,. P 。 §10.5 可压缩形式的压力校正方程
不可压流密度与压力解耦，压力校正只与速度校正相关联；对于可压缩流动则不然，压力不仅出现在动量方程中，而且通过状态方程直接与密度以及能量方程相关联。密度的变化也会影响压力变化，所以必须建立压力-速度-密度的耦合关系式。
令'*'*'*,,e
e e e e e e e e v v v u u u ρρρ+=+=+=，速度-压力校正关系如(33)(35)式所示： 密度-压力校正关系为()''p p ?=ρρ，对于等熵流动有:
()()()211p kRT kp ρρ??=== (41)
则界面e 处对流通量为：
()()()()()()
*'*'*'**
'''**'''**''2*'
*
()()() 1P E BC BC
e e e e e z
e e r BC BC BC BC BC BC e e e z e r e e z e r e e z e r e e P e e e e
e F u u S v v S F u S v S u S v S u S v S F A p p S F p kp ρρρρρρ??=++++??=++++++=+-+ (42) 对界面上的压力修正应用上风格式，即：
()))*'*'**
'*
*
,0,0e e
e
P
e P
E
e E
F p kp p
F kp p
F kp ????=--??
??
(43)
所以
()(){
}()(){}**2*
*
'
*2
**'
1,0 1,0P
E
e e e P e e P
e e
P e
e
E
e
F F A S F kp p
A S
F kp p
ρρ??=++??
??-+-?? (44)
同样可得：
()(){
}
()(){}**2**'
*2
*
*'
1,0 1,0P
N
n n n P n n P n n
P n
n
N
n
F F A S F kp p A S
F kp p
ρρ??=++????-+-??
(45)
最终可由连续方程
()00p
p e w n s V
F F F F t
ρ
ρ-+-+-=?导出如下压力校正方程：
''
'''P P E E W W N N S P a P a P a P a P a P b =++++
(46)
式中
())())())()()()()()(){
}()2
2
2
2
2222
1,01,01,01,01111 ,0,0,0,0
E e P e e E
e W w
P w
w
W w
N
n
P n
n
N
n
S
s
P s
s
S
s
P e P e w P w n P
n s
P
s
e w n s
e
w n s P
a A S F k p
a A S F k p a A S F k p a A S F k p a A S A S A S A S
F F F F k p ρρρρρρρρ??=+-??
??=+????=+-????=+??
=+++????????++-++-????????
(46a)
()0
P
P w e s n V b F F F F t
ρρ-=
+-+-?
(46b)
上式中上标*被省略，交界面上的密度值e n ρρ、等应采用插值方法求出。从(46a)可看出，该方程的系数矩阵不再具有对角占优特性，且为病态方程组，其求解出现了新的困难。
实验3离散系统的差分方程、冲激响应和卷积分析 一 一、实验目的 1 加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。 二、实验原理 离散系统 ] [n x ] [n y Discrete-time systme 其输入、输出关系可用以下差分方程描述：∑=∑=-=-M k k N k k k n x p k n y d 0 ][][ 输入信号分解为冲激信号 ∑-=∞ -∞=m m n m x n x ][][][δ 记系统单位冲激响应 ] [][n h n →δ 则系统响应为如下的卷积计算式 ∑∞ -∞ =-= *=m m n h m x n h n x n y ][][][][][ 当 N k d k ,...2,1,0==时，h[n]是有限长度的（n ：[0，M]），称系统为 FIR 系统；反之，称系统为IIR 系统。 在MATLAB 中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv(x,h)计算卷积。 二、实验内容 编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。 ] 1[][]2[125.0]1[75.0][--=-+-+n x n x n y n y n y ]}4[]3[]2[]1[{25.0][-+-+-+-=n x n x n x n x n y 程序1： A=[1,0.75,0.125];B=[1,-1]; x2n=ones(1,65); x1n=[1,zeros(1,30)]; y1n=filter(B,A,x1n); subplot(2,1,1);y='y1(n)'; stem(y1n,'g','.'); title('单位冲击响应') 课程名称 实验成绩 指导教师 实 验 报 告 院系 班级 学号 姓名 日期
第2章 控制系统的状态方程求解 要点： ① 线性定常状态方程的解 ② 状态转移矩阵的求法 ③ 离散系统状态方程的解 难点： ① 状态转移矩阵的求法 ② 非齐次状态方程的解 一 线性定常系统状态方程的解 1 齐次状态方程的解 考虑n 阶线性定常齐次方程 ? ? ?==0)0()()(x x t Ax t x & （2-1） 的解。 先复习标量微分方程的解。设标量微分方程为 ? ??==0)0(x x ax x & （2-2） 对式（2-2）取拉氏变换得 )()(0s aX X s sX =- 移项 0)()(x s X a s =- 则 a s x s X -= )(
取拉氏反变换，得 00 0!)()(x k at x e t x k k at ∑∞ === 标量微分方程可以认为是矩阵微分方程当n=1时的特征，因此矩阵微分方程的解与标量微分方程应具有形式的不变性，由此得如下定理： 定理2-1 n 阶线性定常齐次状态方程（2-1）的解为 00 0!)()(x k At x e t x k k At ∑∞ === （2-3） 式中，∑∞ ==0 !)(k k At k At e 推论2-1 n 阶线性定常齐次状态方程 ???==00 )()()(x t x t Ax t x & （2-4） 的解为 0)(0 )(x e t x t t A -= （2-5） 齐次状态方程解的物理意义是)(0 t t A e -将系统从初始时刻0t 的初始 状态0x 转移到t 时刻的状态)(t x 。故)(0 t t A e -又称为定常系统的状态转移 矩阵。 (状态转移矩阵有四种求法：即定义（矩阵指数定义）法、拉氏反变换法、特征向量法和凯来-哈密顿（Cayly-Hamilton ）法) 从上面得到两个等式 ∑∞ ==0 !)(k k At k At e ])[(11---=A sI L e At 其中，第一式为矩阵指数定义式，第二式可为At e 的频域求法或拉氏反变换法
信息工程学院实验报告 课程名称：数字信号处理 实验项目名称：离散序列的卷积和系统差分方程的MATLAB 实现 实验时间： 班级：电信131 姓名： 学号：201311404113 一、 实 验 目 的: 熟悉序列的卷积运算及其MATLAB 实现；熟悉离散序列的傅里叶变换理论及其MATLAB 实现；加深对离散系统的差分方程和系统频率响应的理解。 二、实 验 原 理 1、MA TLAB 提供了一个内部函数conv(x,h)来计算两个有限长序列之间的卷积。 2、对于时域离散系统，可用差分方程描述或研究输入、输出之间的关系。对于线性时不变系统，经常用的是线性常系数差分方程。一个N 阶线性常系数差分方程用下式表示： ()() N M i i i i b y n i a x n i ==-=-∑∑ 当 0,1,2,,i b i N == 时，[]h n 是有限长度的，称系统为FIR 系统；反之，称系统为IIR 系统。 在MA TLAB 中，可以用函数filter(a,b,x)求解差分方程，其中参数a,b 分别系统函数的分子和分母多项式的系数。 三、实 验 内 容 与 步 骤 实验内容： 1、已知 1(){1,1,1,1,1}x n =，2(){1,1,1,1,1,1,1}x n =，计算12()()*()y n x n x n =。 2、在0到π区间画出矩形序列 10()R n （其定义见例1-3）的离散时间傅里叶变换（含幅度和相位）。 3、求系统：()0.5((1)(2)(3)(4))y n x n x n x n x n =-+-+-+-的单位冲激响应和阶跃响应。 实验步骤: 1、
第三章 系统的分析——状态方程的解 §3-1线性连续定常齐次方程求解 一、齐次方程和状态转移矩阵的定义 １、齐次方程 状态方程的齐次方程部分反映系统自由运动的状况（即没有输入作用的状况），设系统的状态方程的齐次部分为： )()(t Ax t x =& 线性定常连续系统： Ax x =& 初始条件：00x x t == ２、状态转移矩阵的定义 齐次状态方程Ax x =&有两种常见解法：（1）幂级数法；（2）拉氏变换法。其解为 )0()(x e t x At ?=。其中At e 称为状态转移矩阵（或矩阵指数函数、矩阵指数），记为： At e t =)(φ。 若初始条件为)(0t x ，则状态转移矩阵记为：) (0 0)(t t A e t t -=-Φ 对于线性时变系统，状态转移矩阵写为),(0t t φ，它是时刻t ，t 0的函数。但它一般不能写成指数形式。 （１）幂级数法——直接求解 设Ax x =&的解是t 的向量幂级数 Λ ΛΛΛ+++++=k k t b t b t b b t x 2210)( 式中ΛΛ，，， ，，k b b b b 210都是n 维向量，是待定系数。则当0=t 时， 000b x x t === 为了求其余各系数，将)(t x 求导，并代入)()(t Ax t x =&，得： Λ ΛΛΛ&+++++=-1232132)(k k t kb t b t b b t x )(2210ΛΛΛΛ+++++=k k t b t b t b b A
上式对于所有的t 都成立，故而有： ????? ??????======00 3 230 21201!1!31312121b A k b b A Ab b b A Ab b Ab b K K M 且有：00x b = 故以上系数完全确定，所以有： Λ ΛΛΛ+++++=k k t b t b t b b t x 2210)( ΛΛ++++ +=k k t b A k t b A t Ab b 020200! 1 !21 )0()! 1!21(22x t A k t A At I k k ΛΛ+++++= 定义（矩阵指数或矩阵函数）： ∑∞==+++++=022! 1!1!21K k k k k At t A k t A k t A At I e ΛΛ 则 )0()(x e t x At ?=。 （２）拉氏变换解法 将Ax x =&两端取拉氏变换，有 )()0()(s AX X s sX =- )0()()(X s X A sI =- )0()()(1X A sI s X ?-=- 拉氏反变换，有 )0(])[()(1 1x A sI L t x ?-=--
1.设离散控制系统差分方程为x采样周期T。试求：（1） 系统的脉冲传递函数。（2）系统的频率特性表达式。 解：差分方程两边取Z变换，得 脉冲传递函数 频率特性 2.假设离散系统差分方程为。其中; ，，，。试求:(1)分析系统的稳定性。（2），，。 解：（1）对差分方程两边取Z变换，得 特征方程： 解得：； 由于，即系统稳定。 （2）n=0时， n=1时， n=2时， 3.某离散控制系统的差分方程为，其中： ，，，，，，。试求：（1），。（2）分析稳定性。 解：（1）对差分方程两边Z变换，得 特征方程： 解得：； 由于，所以系统稳定。
（2）n=0时， n=1时。 4.离散控制系统的差分方程为：，其中 ，，时，时。试求：(1)，，。（2）脉冲传递函数。 解：（1）差分方程两边取Z变换，得 特征方程： 解得：； 由于，所以系统稳定。 （2）n=0时， n=1时， n=2时， 5.已知：离散控制系统的差分方程为。试求：脉冲传 递函数。系统频率特性 解：对差分方程Z变换，得 频率特性 6.某离散系统的差分方程为=，其中 ，。试求（1）脉冲传递函数，并分析稳定。（2） ，，。 解：对差分方程两边Z变换，得 （）
特征方程： 解得：； 由于，所以系统稳定。 (2)n=0时， n=1时， n=2时，y 7.已知离散系统的差分方程为，试求：（1）脉冲传递 函数。（2）分析系统稳定性 解：（1）对差分方程两边Z变换，得 (2)特征方程：=0 解得：； 由于，所以系统临界稳定。 8.离散系统差分方程为，其中 ，；。试求：，，。（）分析稳定性。 解：(1)n=0时， n=1时， n=2时， (2)对差分方程两边Z变换，得 特征方程： 解得：； 由于，所以系统稳定。 9.某离散系统差分方程为，其中：， 时，；时，。试求：，，。（2）分析
实验2 离散系统的差分方程、冲激响应和卷积分析 一、实验目的 加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。 二、实验原理 离散系统可表示为 其输入、输出关系可用以下差分方程描述： ∑∑==-=-M k m N k k m n x b k n y a 00][][ 输入信号分解为冲激信号， ∑∞ -∞=-= m m n m x n x ][][][δ。 记系统单位冲激响应 ][][n h n →δ， 则系统响应为如下的卷积计算式： ∑∞ -∞=-= *=m m n h m x n h n x n y ][][][][][ 当N k a k ,...2,1,0==时，h[n]是有限长度的（n ：[0，M]），称系统为FIR 系统；反之，称系统为IIR 系统。 在MATLAB 中，可以用函数y=filter(b,a,x)实现差分方程的仿真，也可以用函数 y=conv(x,h)计算卷积，用y=impz(b,a,N)求系统的冲激响应。 对于N 阶差分方程∑∑==-=-M k m N k k m n x b k n y a 00][][， 1) 当给定函数的系数和输入序列时，差分方程的递推过程在MA TLAB 中用函数y=filter(b,a,x)来实现，其中，b 为右端x 的系数,a 为左端y 的系数，a 0=1。求得的输出序列y 和输入序列x 的长度相等。若x 的长度太短，需要补零。用conv 函数计算能在输入序列后自动补零，而filter 函数不能。 2) MATLAB 中有一个求离散系统脉冲响应的专门函数y=impz(b,a,N)，其中，b 为右端x 的系数,a 为左端y 的系数，a 0=1。N 为要求的点数。键入impz(b,a)，程序将自动给出脉冲响应的曲线。 3) 当输入序列和脉冲响应序列都是以数值方式给出时，可以用MATLAB 中的卷积函数y=conv(x,h)来计算。
Chapter 10 基于有限体积法的控制方程离散 10.1 控制方程 稳态气相通用控制方程如下： () ()()+div ρU div Γgrad S t ρφφφ??=+ ()rv u w U U z r r r θU=ui+vj+wk div ???=??=++ ??? ()()()() z r r u r v w t z r r r r S z r r r ?ρφ??? ρφρφρφ????θ ??φ??φ??φ?????θ?θ+++=Γ?????? Γ+Γ++ ? ? ??????? 上式中φ表示气相场通用变量，在三维柱坐标下可以是轴向速度u 、径向速度v 、 周向速度w 以及焓h ，φ为1时代表连续方程，Γ表示输运系数，S 为气相源项，U 为速度矢量。密度由完全气体状态方程ρRT p =确定，温度由焓温关系式(h)T f =确定。每个方程中的φ、Γ、S 如表10.1所示。 表10.1 气相方程组
§10.2 方程的离散化 §10.2.1 有限体积法 有限差分法：原理比较简单，数学推演和程序编制的工作量较小，但其灵活性较差； 有限单元法：数学推演和程序编制工作要比有限差分法复杂，因而比较费时，但其灵活性较好； 有限体积法：作为有限差分法和有限单元法的中间产物，继承了两者的优点，克服了各自的缺点，因而在流场计算中占有重要的地位。 有限体积法的基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解微分方程对每一控制体积积分，得出离散方程，其中的未知数是网格点上的因变量φ。有限体积法在对控制体积积分时，必须假定φ值在网格点之间的分布。 由有限体积法导出的离散方程具有明确的物理意义，即因变量φ在有限大小的控制体积内的守恒原理，如同推导微分方程时因变量在无限小的控制体积中的守恒原理。有限体积法即使在粗网格条件下，也精确地满足
实验2 离散系统的差分方程、单位脉冲响应和卷积分析 一、 实验目的 1、 熟悉并掌握离散系统的差分方程表示法； 2、 加深对单位脉冲响应和卷积分析方法的理解。 二、 实验原理 (一), 1. 单位采样序列 ???=01 )(n δ 0 0≠=n n 在MATLAB 中可以利用zeros()函数实现。 ; 1)1();,1(==x N zeros x 如果)(n δ在时间轴上延迟了k 个单位，得到)(k n -δ即： ???=-01 )(k n δ 0≠=n k n 2．单位阶跃序列 1()=0 u n ??? 00<≥n n 在MATLAB 中可以利用ones()函数实现。 );,1(N ones x = 3．正弦序列 )/2sin()(?π+=Fs fn A n x 在MATLAB 中 ) /***2sin(*1:0fai Fs n f pi A x N n +=-=
4．复指数序列 n j e n x ?=)( 在MATLAB 中 ) **exp(1:0n w j x N n =-= 5．实指数序列 n a n x =)( 在MATLAB 中 n a x N n .^1:0=-= (二) 在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下: 其输入、输出关系可用以下差分方程描述： 00()()N M i i i i a y n i b x n i ==-=-∑∑ 输入信号分解为单位采样序列的移位加权和，即： ()()()m x n x m n m δ∞ =-∞= -∑ 记系统单位脉冲响应 ()()n h n δ→ 则系统响应为如下的卷积计算式：
6.3 离散化边界积分方程的建立 以二维边界离散化方程的建立为例，重点突出离散化方法的学习。 6.3.1建立Laplace 场的边界离散化方程 电磁场边界元法的通用积分方程 (4) 其中： ?????? ?∈∈∈=域外 光滑的边界上域内D D c i 0 211 设在Laplace 场中的二维边界上一点i 处，有方程： 在二维场的边界线l 上进行离散，将l 划分为许多小段，每段以直线段或曲线段逼近，作为一个单元。设l 点共被分为0N 个单元，其中在第一类边界1l 段上划分了1N 个单元，在第二类边界2l 段上划分了2N 个单元： 210N N N += 作为单元待求量的插值计算方式，可分为几种： ① 恒值单元 同一单元中的待求量u 和 n u ??都设为恒定值 (或称零次插值)，实际上是取单元中点的u 值(或 n u ??值)作为单元的u 值(或n u ??值)。这样，取单 元中点为节点，所以求解变量数等于节点数。 ② 线性单元
它也是直线单元，其u 值在单元两端点之间按线性变化(即线性插值)。单元两端点为单元的节点。 ③ 曲线单元 每单元上的节点数大于2，以多节点拟合的曲线逼近边界单元，以单元节点上的高阶插值函数作为待求位函数近似解。 取最简单的单元——恒值单元为例，介绍边界元离散方法。 按上面的方程对i 单元的“i ”节点离散化 ∑? ∑? ==??= ??+ o j o j N j l N j l i l n u F l n F u u 1 1 d d 2 1 ∑? ∑?=== ??+ 1 1 d d 2 1N j l N j l j j i j o j l F q l n F u u ，?= j l ij l F G d ，上式表示为： 设i 点为i 单元的中点（021N i 、、、 =），有 ()∑∑==== 1 01 21N j N j j ij j ij N i q G u H ,,, 式中： 于是上述0N 个方程写为矩阵形式 GQ HU = 由定解问题中的第一类边界1l ，对应有1N 个单元的位值u s 是已知的，2l 是第二类边界，对应有2N 个单元n u q s ??= 位是已知。所以上述矩阵方程中，有2N 个单元的u 值和 1N 个单元的q 值是未知的，即是说矩阵方程有021N N N =+个未知数。设单元排列顺序 在1l 边界上为1，2，……，1N ，在2l 边是上为11+N 、21+N 、…、0N ，则上述矩阵方
求解系统的状态方程 一、实验设备 PC计算机，MATLAB软件，控制理论实验台 二、实验目的 (1)掌握状态转移矩阵的概念。学会用MATLAB求解状态转移矩阵 (2)学习系统齐次、非齐次状态方程求解的方法，计算矩阵指数，求状态响应； (3)通过编程、上机调试，掌握求解系统状态方程的方法，学会绘制输出响应和状态响应曲线； (4)掌握利用MATLAB导出连续状态空间模型的离散化模型的方法。 三、实验原理及相关基础 (1)参考教材P99~101“3.8利用MATLAB求解系统的状态方程” (2)MATLAB现代控制理论仿真实验基础 （3）控制理论实验台使用指导 四、实验内容 (1)求下列系统矩阵A对应的状态转移矩阵 （a）
（b） 代码： syms lambda A=[lambda 0 0;0 lambda 0;0 0 lambda];syms t;f=expm(A*t) （c） 代码： syms t;syms lambda;A=[lambda 0 0 0;0 lambda 1 0;0 0 lambda 1;0 0 0 lambda];f=expm(A*t) (2) 已知系统
a) 用MATLAB求状态方程的解析解。选择时间向量t，绘制系统的状态响应曲线。观察并记录这些曲线。 （1） 代码： A=[0 1; -2 -3]; B=[3;0]; C=[1 1]; D=[0]; u=1; syms t; f=expm(A*t);%状态转移矩阵 x0=0; s1=f*B*u; s2=int(s1,t,0,t)%状态方程解析解 状态曲线： （2）A=[0 1;-2 -3]; syms t; f=expm(A*t); X0=[1;0]; t=[0:0.5:10]; for i=1:length(t); g(i)=double(subs(f(1),t(i))); end plot(t,g)
现设线性时变系统的离散状态方程和观测方程为： X(k) = F(k,k-1)·X(k-1)+T(k,k-1)·U(k-1) Y(k) = H(k)·X(k)+N(k) 其中 X(k)和Y(k)分别是k时刻的状态矢量和观测矢量 F(k,k-1)为状态转移矩阵 U(k)为k时刻动态噪声 T(k,k-1)为系统控制矩阵 H(k)为k时刻观测矩阵 N(k)为k时刻观测噪声 则卡尔曼滤波的算法流程为： 预估计X(k)^= F(k,k-1)·X(k-1) 计算预估计协方差矩阵 C(k)^=F(k,k-1)×C(k)×F(k,k-1)'+T(k,k-1)×Q(k)×T(k,k-1)' Q(k) = U(k)×U(k)' 计算卡尔曼增益矩阵 K(k) = C(k)^×H(k)'×[H(k)×C(k)^×H(k)'+R(k)]^(-1) R(k) = N(k)×N(k)' 更新估计 X(k)~=X(k)^+K(k)×[Y(k)-H(k)×X(k)^] 计算更新后估计协防差矩阵 C(k)~ = [I-K(k)×H(k)]×C(k)^×[I-K(k)×H(k)]'+K(k)×R(k)×K(k)' X(k+1) = X(k)~ C(k+1) = C(k)~ 重复以上步骤 ＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
Ｍatlab实现代码 ******************************************************************************* ************************************************** %%%% Constant Velocity Model Kalman Filter Simulation %%%% %========================================================================== clear all; close all; clc; %% Initial condition ts = 1; % Sampling time t = [0:ts:100]; T = length(t); %% Initial state x = [0 40 0 20]'; x_hat = [0 0 0 0]'; %% Process noise covariance q = 5 Q = q*eye(2); %% Measurement noise covariance r = 5 R = r*eye(2); %% Process and measurement noise w = sqrt(Q)*randn(2,T); % Process noise v = sqrt(R)*randn(2,T); % Measurement noise %% Estimate error covariance initialization p = 5; P(:,:,1) = p*eye(4); %========================================================================== %% Continuous-time state space model %{ x_dot(t) = Ax(t)+Bu(t) z(t) = Cx(t)+Dn(t) %} A = [0 1 0 0; 0 0 0 0;
求解离散化结构动力方程 问题描述： 图（1）离散化质量弹簧系统示意图 取三自由度系统研究。质量、弹簧、阻尼连接关系如图（1）所示。各质点所受力为 sin(),(1,2,3)i i F f t i ω== ，各质点离开静平衡位置的距离为i x ，初始条件为 ()(),0,00,0i i i i x x x x == ，上标表示对时间求导，下标i 表示第i 质点，第二下标表示为初 始时刻。 振型叠加法 为降低计算量，得出较简洁结果，特别取 123123412 34m m m m k k k k k c c c c c ===?? ====??====? (1) 由牛顿第二定律可得系统振动微分方程 MX +CX +KX =F (2) 其中 1 23000 00 m m m m ????=???????M =I ， 12 22 2333340 202002c c c c c c c c c c c c c c c c c +--???? ????-+-=--? ???????-+-???? C = ， 12 22 23 33 34020200 2k k k k k k k k k k k k k k k k k +--????????-+-=--??? ?????-+-???? K = []1 2 3x x x T X = ，[]123sin() sin() sin()f t f t f t ωωωT F = 其固有频率为 123Ω= Ω =Ω=正则振型矩阵为
11 11 ?? ?? ?? u= 令 ()() t t =? X uη(3) 代入(2)式，并在方程两边左乘T u，则正则坐标的振动微分方程为 ()()()() t t t t +?+?= T T T ηu Cuηu Kuηu F(4) 记 ( ( 11 * 22 33 2 00 00 2 0000 00 2 00 c m c c c m c c m ?? - ?? ???? ???? ?? ==?? ???? ?? ?? ?? ?? ?? T C u Cu=(5) ( ( 11 * 22 33 2 00 00 2 0000 00 2 00 k m k k k m k k m ?? - ?? ???? ???? ?? ==?? ???? ?? ?? + ?? ?? ?? T K u Ku=(6) ()() () * 1231 ** 132 * 3 123 sin sin f f f t t t t f f f f ω ?? ++?? ???? ==-= ??? ??? -+?? ? ?? T F u F(7) 方程(4)为 ()()()()() *sin,1,2,3 ii ii i t c t k t f t i ω +?+?== ηηη(8) (8)式已解耦，各式为常微分方程可独立求解 ( )()) ()()() () ,1,2 *2 2 222 exp exp 22 cos sin i i ii i ii i ii ii ii ii t t t C c C c f c t k t c k η ωωωω ωω ???? =+- ???? ???? ?? +- ?? +- (9)
第三章 控制方程及求解方法 本章介绍FDS 的控制方程并概述一个通用的求解方法。在后面的章节中将详细描述个别方程。控制方程是一组经过适当的简化和近似的偏微分方程。该数值计算方法基本上是由一个有限差分近似方程和一个及时更新这些方程的程序组成的。 3.1 控制方程 本节介绍了牛顿流体的基本的质量、动量和能量守恒方程。这些方程式可以在几乎任何关于流体动力学或CFD 的书籍中找到。Anderson 等人[13]的著作对方程的描述是一个特别有用的参考，方程所用的符号、各种近似都是从他们那引用的。请注意，这是一组由6个方程组成的偏微分方程组，共有6个未知数，三维空间的所有函数和时间：密度ρ，三个方向的速度分量u =[u ，v ，w] T ，温度T ，压力p 。 3.1.1 质量输运 用密度表示的的质量守恒定律： 或以单独的气态物质αY 表示的质量守恒定律： 这里∑'''='''αα，b b m m 表示的是物种通过蒸发液滴或颗粒的生产速率。由于1Y =∑α和 0m ='''∑α 和 b b m m '''='''∑ α， ，并且根据定义它是假定0Y D =?∑ααρ的，因此要对所有物种产生的原始的质量守恒方程进行总结。一般来说，这最后的说法是不正确的。然而， 输运方程是求解总质量而不是单个种类的质量，这意味着未知物体种类的扩散系数的选取，使所有的扩散通量的总和为零。 3.1.2 动量输运 动量守恒方程的形式为： uu 项表示一个双值张量。在矩阵表示法中，u =[u ，v ，w]T ，双值由张量乘以向量u 和u T 得到的。?·ρuu 项是通过应用矢量运算子?=? ?? ? ????????z y x ，，而形成的矢量。动量方程中的f b 项表示外力，如液滴所受到的阻力。应力张量τ ij 定义为： S ij 项为用传统的张量表示法表示的对称张量应变率，符号μ表示流体的动力粘度。 整体计算可以被视为直接数值模拟（DNS ）或作为一个大涡模拟（LES ），直接数值模
HEFEI UNIVERSITY 信号与系统项目设计报告 系别电子信息与电子工程系 题目项目第十题 专业电子信息工程 班级 11电子信息工程（2）班 小组成员钟文俊（1105012012）谢伟明（1105012041)授课老师纪平 完成时间 2014.01.02
用MATLAB仿真离散系统差分方程 一、设计题目 -n f - + - n n n y y n = y f (- ( )2 ) ( )1 ( 5.0 .0 ) 25 ( )1 + 系统输入序列为) f nε =。 n (n )5.0( ) ( 二、设计要求 1、试用MATLAB绘出输入序列的时域波形； 2、用MATLAB求出该系统0～20区间的样值； 3、用MATLAB画出系统的零状态响应波形。 三、功能分析 差分方程反映的是关于离散变量的取值与变化规律。通过建立一个或几个离散变量取值所满足的平衡关系，从而建立差分方程。差分方程就是针对要解决的目标，引路系统或过程中的离散变量，根据实际背景的规律、性质、平衡关系，建立离散变量所满足的平衡关系等式，从而建立差分方程。通过求出和分析方程的解，或者分析得到方程解的特别性质（平衡性、稳定性、渐近性、振动性周期性等），从而把握这个离散变量的变化过程的规律，进步再结合其他分析，得到原问题的解。 四、设计原理分析 1、差分方程定义 含有未知函数yt=f(t)以及yt的差分Dyt， D2yt，…的函数方程，称为常差分方程(简称差分方程)；出现在差分方程中的差分的最高阶数，称为差分方程的阶。n阶差分方程的一般形式为F(t，yt，Dyt，…， Dnyt)=0，其中F是t，yt, Dyt，…， Dnyt的已知函数，且Dnyt一定要在方程中出现。 含有两个或两个以上函数值yt，yt+1，…的函数方程，称为(常)差分方程，出现在差分方程中未知函数下标的最大差，称为差分方程的阶。n阶差分方程的一般形式为F(t，yt，yt+1，…，yt+n)=0,其中F为t，yt，yt+1，…，yt+n的已知函数，且yt和yt+n一定要在差分方程中出现。 2、差分方程的意义与应用 差分方程模型有着广泛的应用。实际上，连续变量可以用离散变量来近似逼
储运与建筑工程学院能源与动力工程系 计算传热学课程大作业报告 作业题目：五种对流离散项的对比研究 学生姓名：宋龙 学号：6030221 专业班级：能动－班 2017年 11 月 3 日
目录 1 计算题目 (1) 2 数学物理模型 (2) 3 计算区域及方程离散 (3) 3.1 区域离散 (3) 3.2 方程离散 (3) 3.2.1 中心差分格式 (3) 3.2.1 迎风差分格式 (4) 3.2.3 混合格式 (4) 3.2.4 指数格式 (5) 3.2.5 指数格式 (6) 3.2.6 五种格式格式系数aEDe的表达式 (6) 4 数值方法及程序流程 (7) 5 计算结果验证及网格独立性考核 (8) 5.1 计算结果验证 (8) 5.1 网格独立性考核 (11) 6 结果分析与讨论 (14) 7 参考资料 (15) 附录 (16) 附录A：计算环境及源程序 (16)
1计算题目对于有源项的一维稳态空气对流-扩散传热方程： d dx ρuT= d dx λ c dT dx +S 设源项S=0.5-100x，利用中心差分格式，一阶迎风格式，混合格式，指数格式，乘方格式求解在不同流速情况下，温度T的一维分布。
2数学物理模型物理模型：一维，稳态，有内热源，常物性 控制方程:d dx ρuT=d dx λ c dT dx +Sλ c =Γ 物理条件：ρ,Γ=const s=0.5-100x 边界条件：x=0 T=300K; x=1 T=500K 由于在常物性均分网格的情况下网格贝克勒数与流速成正比，所以在流速不同的情况下求温度场即在网格贝克勒数不同的情况下求温度场。
实验报告 实验课程：数字信号处理 实验内容：实验2离散系统的差分方程、冲激 响应和卷积分析 2013年6 月3 日
一、实验目的：加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。 二、实验原理： 离散系统 其输入、输出关系可用以下差分方程描述： ∑=∑=-= -M k k N k k k n x p k n y d 00][][ 输入信号分解为冲激信号，∑-= ∞-∞=m m n m x n x ][][][δ。 记系统单位冲激响应 ][][n h n →δ， 则系统响应为如下的卷积计算式： ∑∞-∞=-= *=m m n h m x n h n x n y ][][][][][ 当N k d k ,...2,1,0==时，h[n]是有限长度的（n ：[0，M]），称系统为 FIR 系统；反之，称系统为IIR 系统。 在MATLAB 中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv(x,h)计算卷积。 三、实验内容及步骤： 1、实验内容： 分别在x （n ）=δ(n)和x （n ）=cos(2π*0.47*n)的输入下，编制程序分别用上述两种 方法求解下列两个系统的响应，并得出系统零极点分布图，绘出其图形。 ]1[][]2[125.0]1[75.0][--=-+-+n x n x n y n y n y ⑴ ]}4[]3[]2[]1[{25.0][-+-+-+-=n x n x n x n x n y ⑵ 2、实验代码及结果 ]1[][]2[125.0]1[75.0][--=-+-+n x n x n y n y n y ⑴
实验二 差分方程的求解和离散系统频率响应的描述 一、 实验目的 1、掌握用MATLAB 求解差分方程的方法。 2、掌握绘制系统的零极点分布图和系统的频率响应特性曲线的方法。 3、 观察给定系统的冲激响应、阶跃相应以及系统的幅频特性和相频特性 二、 实验内容 1、已知描述离散新天地差分方程为：y(n+2)-0,25y(n+1)+0.5y(n)=x(n)+x(n-1)，且知该系统输入序列为)()2/1()(n u n x n =，试用MATLAB 实现下列分析过程：画出输入序列的时序波形；求出系统零状态响应在0～20区间的样值；画出系统的零状态响应波形图。 2、一离散时间系统的系统函数：5731053)(2323-+-+-=z z z z z z z H ，试用MA TLAB 求出系 统的零极点；绘出系统的零极点分布图；绘出响应的单位阶跃响应波形。 三、 实验报告要求 1、求出各部分的理论计算值， 并与实验结果相比较。 2、绘出实验结果波形（或曲线），并进行分析。 3、写出实验心得。 附录：本实验中所要用到的MATLAB 命令 1、系统函数H(z) 在MATLAB 中可调用函数zplane （），画出零极点分布图。调用格式为： zplane （b,a ） 其中a 为H （z ）分母的系数矩阵，b 为H(z)分子的系数矩阵。 例2－1：一个因果系统：y （n ）－0.8y(n －1)＝x(n) 由差分方程可求系统函数 8.0,8.011)(1>-= -z z z H 零极点分布图程序： b=[1,0]; a=[1,-0.8]; zplane(b,a)
2、求解差分方程 在MA TLAB中，已知差分方程的系数、输入、初始条件，调用filter（）函数解差分方程。 调用filter（）函数的格式为：y=filtier(b,a,x,xic)，参数x为输入向量（序列），b,a分别为（1-30）式中的差分方程系数，xic是等效初始状态输入数组（序列）。 确定等效初始状态输入数组xic（n），可使用Signal Processing toolbox中的filtic（）函数，调用格式为：y=filtic(b,a,y,x) 。其中y=[y(-1),y(-2),…,y(-N)]，x=[x(-1),x(-2),…,x(-M)] 。 例2－2：已知差分方程2y(n)-3y(n-1)+y(n-2)=2x(n) ，式中x(n)=(1/4)n u(n) ，y(-1)=4 ,y(-2)=10 ,求全响应y(n) 。 MATLAB程序如下： n=[0:7]; x=(1/4).^n; a=[2,-3,1]; b=[2]; y=[4,10]; xic=filtic(b,a,y) y1=filter(b,a,x,xic) y2=(1/3)*(1/4).^n+(1/2).^n+(2/3)*ones(1,8) %这是直接将差分方程Z变换后代入X(z)求出Y(z)，反变换后求出x(n)。 执行结果为： xic = 1 -2 y1 = 2.0000 1.2500 0.9375 0.7969 0.7305 0.6982 0.6824 0.6745 y2 = 2.0000 1.2500 0.9375 0.7969 0.7305 0.6982 0.6824 0.6745 3、求系统的冲激响应和阶跃响应 ⑴在MATLAB中，有专门求冲激响应并绘制其时域波形的函数impz( ) 格式： y=impz(b,a,n) %这是求数值解
02节airpak高级班48讲-《CFD流体流动控制方程》 各位同学，大家好，我是七师兄，今天我们来学习Airpak高级班的第二节课，《流体流动控制方程》。 室内空气流动数值计算基本有如下几个环节:建立流动的数学物理模型,得出流动的控制微分方程组:在计算域上对流动控制微分方程组进行离散,将离散结果整理为代数方程组:拟定特定的算法,然后据算法对离散所得代数方程组进行求解,从而得到计算域内空气流动的分布信息。 由于我们Airpak软件研究的一个重点方向是建筑室内空气流动的情况，那么我们首先来看下室内空气流动数值计算的物理模型。根据热工理论基础,可认为室内空气满足气体状态方程,即 p=pRT 式中,p为空气压力(Pa);p为空气密度(kg/m');R为空气常数,约287]/(kg.K);T为空气热力学温度(K)。 其次,通常可将室内空气流动的压力视为常数,于是可得: pr=常数 另外,常见的室内环境中的空气流动基本为低速流动,流速常在10~20m/s以下,因此叮将室内空气当作不可压缩流休看待,即 v.U=0 而且,空调通风房间内空气温度变化不大,也即密度变化不大,因此可认为室内空气流动符合Boussinesqu（布西内斯克）假设:密度变化并不显著改变流体性质,动量守桓中,密度的变化对惯性力项、压力差项和粘性力项的影响可忽略不计,而仅考虑对质量力项的影响"。综上所述,除了对质量力项的考虑之外,室内空气物性都可当作常物性看待。 最后,室内空气的粘性不可忽略,必须用粘性流体动力学的理论来研究。进一步而言,室内空气流动通常都是湍流流动,需要相应的湍流理论来模拟。而且.由于空调房间通常是通过送入与室内空气温度不同的冷风或者热风进行室内空气温度调节的,或者由散热器、冷辐射吊顶、电热膜等通过辐射换热方式室内空气温度,因此室内空气流动往往是自然对流和强迫对流并存的混合对流或者是辐射换
实验2离散系统的差分方程、冲激响应和卷积分析 一、实验目的 1 加深对离散系统的差分方程、冲激响应和卷积分析方法的理解。 二、实验原理 离散系统 其输入、输出关系可用以下差分方程描述： 输入信号分解为冲激信号 记系统单位冲激响应 则系统响应为如下的卷积计算式 当时，h[n]是有限长度的（n：[0，M]），称系统为FIR系统；反之，称系统为IIR系统。 在MATLAB中，可以用函数y=Filter(p,d,x) 求解差分方程，也可以用函数 y=Conv(x,h)计算卷积。 二、实验内容 编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。 y[n]+0.6y[n-1]+0.08y[n-2]=x[n]-x[n-1] y[n]=0.2{x[n-1]+x[n-2]+x[n-3]+x[n-4]+x[n-5]} 实验 程序1： A=[1,0.6,0.08];B=[1,-1]; x2n=ones(1,65); x1n=[1,zeros(1,30)]; y1n=filter(B,A,x1n); subplot(2,1,1);y='y1(n)'; stem(y1n,'g','.'); title('单位冲击响应') y2n=filter(B,A,x2n); subplot(2,1,2); y='y2(n)'; stem(y2n,'g','.'); title('阶跃响应')
程序2 A=[1];B=[0,0.2,0.2,0.2,0.2]; x2n=ones(1,25); x1n=[1,zeros(1,30)]; y1n=filter(B,A,x1n); subplot(2,1,1);y='y1(n)'; stem(y1n,'g','.'); title('单位冲击响应') y2n=filter(B,A,x2n); subplot(2,1,2); y='y2(n)'; stem(y2n,'g','.'); title('阶跃响应')
第三章 线性系统的运动分析 §3-1线性连续定常齐次方程求解 一、齐次方程和状态转移矩阵的定义 １、齐次方程 状态方程的齐次方程部分反映系统自由运动的状况（即没有输入作用的状况），设系统 的状态方程的齐次部分为：)()(t Ax t x = 线性定常连续系统：Ax x = ２、状态转移矩阵的定义 齐次状态方程Ax x = 有两种常见解法：（1）幂级数法；（2）拉氏变换法。其解为 )0()(x e t x At ?=。其中At e 称为状态转移矩阵（或矩阵指数函数、矩阵指数），记为：At e t =)(φ。 若初始条件为)(0t x ，则状态转移矩阵记为：)(00)(t t A e t t -=-Φ 对于线性时变系统，状态转移矩阵写为),(0t t φ，它是时刻t ，t 0的函数。但它一般不能写成指数形式。 （１）幂级数法 设Ax x = 的解是t 的向量幂级数 +++++=k k t b t b t b b t x 2210)( 式中 ，，， ，，k b b b b 210都是n 维向量，则 +++++=-1232132)(k k t kb t b t b b t x )(2210 +++++=k k t b t b t b b A 故而有： ????? ??????======00323021 201!1!31312121b A k b b A Ab b b A Ab b Ab b K K
且有0)0(b x =。 故 +++++=k k t b t b t b b t x 2210)( +++++=k k t b A k t b A t Ab b 020200!1 !21 )0()! 1!21(22x t A k t A At I k k +++++= 定义：∑∞==+++++=022! 1 !1!21K k k k k At t A k t A k t A At I e 则)0()(x e t x At ?=。 （２）拉氏变换解法 将Ax x = 两端取拉氏变换，有 )()0()(s Ax x s sx =- )0()()(x s x A sI =- )0()()(1 x A sI s x ?-=- 拉氏反变换，有 )0(])[()(11 x A sI L t x ?-=-- 则 ])[()(11---==A sI L e t At φ 【例3.1.1】 已知系统的状态方程为x x ?? ? ???=0010 ，初始条件为)0(x ，试求状态转移矩阵和状态方程的解。 解：（１）求状态转移矩阵 ++++ +==k k At t A k t A At I e t ! 1 !21)(22φ 此题中： ??????=0010A ， ?? ????====000032n A A A 所以