图2连续系统的冲激响应2
(3)impulse(num, den, t1:dt:t2) :绘出在t1~t2时间范围内,且以时间间隔dt均匀取样的冲激响应波形。对上例,若运行命令impulse(num, den, 1:0.1:2),则绘出1~2 秒内,每隔0.1秒取样的冲激响应的时域波形,如图 3 所示
图3 连续系统的冲激响应3
(4)y = impulse(num, den, t1:dt:t2):不绘出波形,而是求出系统冲激响应的数值解。对上例,若运行命令y=impulse(num, den, 0:0.2:2),则运行结果为:
y =
3.0000 1.1604 0.3110 -0.0477 -0.1726
-0.1928 -0.1716 -0.1383 -0.1054 -0.0777 -0.0559
2.step()函数:可绘出连续系统的阶跃响应g(t)在指定时间范围的时域波形并能求出其数值解,和impulse()函数一样,也有四种调用格式。
3、lsim()函数:前面介绍用impulse()和step()函数求取系统单位冲激响应和单位阶跃响应的方法,但如果输入信号由其他数学函数描述,则这两个函数就无能为力了,需要借助lsim()函数来绘制系统的时域响应曲线。lsim()函数的调用格式与impulse()函数的调用格式类似,所不同的是需要提供有关的输入信号向量,调用格式为:
lsim(num, den, u, t)
其中向量den和num表示连续系统常微分方程响应和激励的各阶导数的系数,u和t用于描述输入信号,u中的点对应t中各时间点处的输入值。调用这一函数可以绘制系统在给定输入作用下的零状态响应。
例:求r?(t)?4r(t)?4e(t)在矩形周期信号激励下的零状态响应。 t = 0: 0.01 : 10*pi; u = square(t); den = [1 4]; num = [4]; lsim(num, den, u, t); axis([0, t(end), -1.1, 1.1])
2.3.2 利用MATLAB 中的Simulink 进行系统时域特性仿真 1、利用系统传输函数进行时域分析
系统传输函数 H(s)可以表示为分子和分母多项式的形式,也可以表示为零极点形式,例如:
H(s)?s?1 2s?1.3s?0.8
在continuous的子库中选择传输函数(TransferFcn)子库,用鼠标把传输函数模块拖动untitled视窗中,置于激励信号源和示波器之间,然后用鼠标拖出的连线将信号源、传输函数模块和示波器等按照系统的要求连接起来即可。
如果需要对传输函数的参数进行设置,则可通过双击传输函数模块打开它的设置环境窗,可以分别设置分子、分母多项式的系数和阶数等。在untitled 窗口的菜单中选simulation 的start 功能则可执行仿真,待仿真完毕后,则可通过激活示波器观察到激励信号和响应信号的波形。
在source库中没有单位冲激信号模块,对于线性系统可以由系统的单位阶跃响应经过微分来单位冲击响应信号。
对于系统在任意激励下的零状态响应,可以采用的方法类似1中仿真系统的冲激响应,只是输入是任意波形例如周期的矩形波等。
2、根据系统方框图构建线性微分方程
除利用数学表达式描述系统模型之外,也可借助方框图表示系统模型。每个方框图反映某种数学运算功能,给出该方框图输出与输入信号的约束条件,若干个方框图组成一个完整的系统。对于线性微分方程描述的系统,它的基本运算单元是相加、倍乘(标量乘法)和积分(或微分)。虽然可以采用微分运算和其他运算来构成一个系统,弹在实际应用中考虑到抑制突发干扰(噪声)信号的影响,一般选用积分运算单元与其他运算来构成一个系统。
在Matlab的Simulink中相应的运算单元如下:
2.3.3 用符号运算工具箱求常微分方程的解析解
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解。如果有初始条件,则求出特解。MATLAB常微分方程符号解的语法是:
dsolve('equation', 'condition')
其中,equation代表常微分方程式,且以Dy代表一阶微分项y',D2y代表.一阶微分项y\ condition则为初始条件。
dsolve的调用格式如下; (1)dsolve('equation')
给出微分方程的解析解,表示为t的函数; (2)dsolve('equation', 'condition')
给出微分.方程初值问题的解,表示为t的函数; (3)dsolve('equation', 'v')
给出微分方程的解析解,表示为v的函数; (4)ddsolve('equation', 'condition', 'v')
给出微分方程初值问题的解,.表示为v的函数。
例:已知描述某二阶线性时不变连续时间系统的动态方程为:
y??(t)?6y?(t)?8y(t)?f(t),t?0
初始条件y(0)=1,y’(0)=2,输入信号f(t)=e-tu(t),求系统的零状态响应、零输入响应及全响应。
解:Matlab 零状态响应
dsolve('D2y+6*Dy+8*y = exp(-t)', 'y(0) = 0','Dy(0)=0') ans =
1/(3*exp(t)) - 1/(2*exp(2*t)) + 1/(6*exp(4*t)) 零输入响应
dsolve('D2y+6*Dy+8*y = 0', 'y(0) = 1','Dy(0)=2') ans =
3/exp(2*t) - 2/exp(4*t) 全响应
dsolve('D2y+6*Dy+8*y = exp(-t)', 'y(0) = 1','Dy(0)=2') ans =
1/(3*exp(t)) + 5/(2*exp(2*t)) - 11/(6*exp(4*t)) 2.3.4 用数值方法求常微分方程的数值解
在求常微分方程数值解方面,Matlab具有丰富的函数,将其统称为solver,其一般格式为:
[t,y] = solver(odefun, tspan, y0)
该函数表示在区间tspan = [t0,tf]上用初始条件声求解显式常微分方程y’=f(t,y)。
odefun为显式常微分方程y’=f(t,y)中的f(t,y)。tspan为求解区间,要获得问题在其他指定点t0,t1, t2, ……, tf的解,则令tspan = t0,t1, t2, ……, tf(要求ti单调),y0为初始条件。
solver为命令ode45, ode32, ode113, ode15s, ode23s, ode23t,ode23tb之一。其中ode45, ode32, ode113属于非刚性ODE类型,这些命令的特点如.表1所示。ode15s, ode23s, ode23t,ode23tb属于刚性ODE类型,这些命令的特点如表2所示。
函数ode45与ode23是常使用的求解方法,函数ode45的使用与ode32完全一样。两个函数的差别在于必须与所用的内部算法相关。两个函数都运用了基本的龙格-库塔(Runge-Kutta)数值积分法的变形。
ode32运用一个组合的2/3阶龙格-库塔-芬尔格(Runge-Kutta-Fehlerg)算法,而ode45运用组合的4/5l阶龙格-库塔-芬尔格算法。
一般地,ode45可取较多的时间步。所以,要保持与ode23相同误差时,在t0和tf之间可取较少的时间步。然而,在同一时间内ode23每时间步至少调用3次,而ode45每时间步至少调用6次。
正如使用高阶多项式内插常常得不到最好的结果一样,ode45也不总是比ode23好。如果ode45产生的结果,对作图间隔太大,则必须在更细的时间区间内,对数据进行内插,比如用函数interp1。这个附加时间点会使ode23更有效。作为一条普遍规则,在所计算的导数中,如有重复的不连续点,为保持精度致使高阶算法减少时间步长,这时低阶算法更有效。
采用求解器solver求解ODE的基本过程如下:
(1)根据问题所属学科中的规律、定律、公式,用微分方程与初始条件进行描述。F(y,y',y\(n), t)=0
y(0)=y0,y’(0)=y1,y\2, y(n-1)(0)=yn-1 写为向量的形式为