连续系统的数字仿真1_08解析
合集下载
第09讲 连续系统的数字仿真
设置完以上仿真控制参数后，则可选择 Simulation→Start 命令来启动仿真过程 , 在仿真结 束时会自动发出一声鸣叫。在仿真过程中，还允许 采用Simulation菜单下的Pause和Continue命令来暂 停或继续仿真过程，若选择Stop命令,则人为中止仿 真过程。结果分析有助于模型的改进和完善，同时 结果分析也是仿真的主要目的。
9
ode45：四/五阶龙格－库塔法，它采用的是单步法， 也就是说它在计算y(tn)时，它仅需要最近处理时刻的 结果y(tn-1)。该算法对于大多数系统有效，最常用， 但不适用于刚性（stiff）系统。 discrtet：当Simulink检测到模型中没有连续状态 时所选择的一种算法。
10
定步长(Fixed-step)算法
15
单任务模式(Single
Tasking)
该模式不检查模块间的采样速率转换。该模式 对于建造单任务系统模型非常有用，在此类系统中， 任务同步不是问题。
自动模式(Auto)
当选用此模式时，如果模型中所有模块运行于 同样的采样率下， Simulink 使用单任务模式；如果 模型包含有不同采样速率运行的模块，则使用多任 务模式。
25
(1) 利用示波器模块(Scope)得到输出结果
当利用示波器模块作输出时，它不仅会自动地将 仿真的结果从示波器上实时地显示出来。而且也可同 时把示波器缓冲区存储的数据，送到MATLAB工作 空间指定的变量中保存起来，以便利用绘图命令在 MATLAB命令窗口里绘制出图形。
26
在示波器模块的窗口中，利用快捷按钮“ ”， 可打开如图6-35所示的示波器模块参数(parameters) 对话框。示波器参数对话框中有两个页面，图6-35 （a）为一般参数设置（General），图6-35（b）为 数据存储参数设置（Data history）。
9
ode45：四/五阶龙格－库塔法，它采用的是单步法， 也就是说它在计算y(tn)时，它仅需要最近处理时刻的 结果y(tn-1)。该算法对于大多数系统有效，最常用， 但不适用于刚性（stiff）系统。 discrtet：当Simulink检测到模型中没有连续状态 时所选择的一种算法。
10
定步长(Fixed-step)算法
15
单任务模式(Single
Tasking)
该模式不检查模块间的采样速率转换。该模式 对于建造单任务系统模型非常有用，在此类系统中， 任务同步不是问题。
自动模式(Auto)
当选用此模式时，如果模型中所有模块运行于 同样的采样率下， Simulink 使用单任务模式；如果 模型包含有不同采样速率运行的模块，则使用多任 务模式。
25
(1) 利用示波器模块(Scope)得到输出结果
当利用示波器模块作输出时，它不仅会自动地将 仿真的结果从示波器上实时地显示出来。而且也可同 时把示波器缓冲区存储的数据，送到MATLAB工作 空间指定的变量中保存起来，以便利用绘图命令在 MATLAB命令窗口里绘制出图形。
26
在示波器模块的窗口中，利用快捷按钮“ ”， 可打开如图6-35所示的示波器模块参数(parameters) 对话框。示波器参数对话框中有两个页面，图6-35 （a）为一般参数设置（General），图6-35（b）为 数据存储参数设置（Data history）。
第六讲 连续系统仿真概述祥解
四个一阶微分方程，有4个状态变量y(1) y(2)…….。 计算机仿真技术
26
今希望用直角坐标输出，故还引入两个定义变量 x y(21)，y)* cos[ y(2)] y(22) y(1)* sin[ y(2)]
根据卫星的发射速度，可以建立起方程组的初始条件： y(1)0=6400km(卫星到地心的距离，即为地球之半径)，
任务失败。
计算机仿真技术
6
Xb＝Vb*t
Yb=100
Xf’=Vf*cosθ
Yf’=Vf*sinθ
tanθ=(100-Yf)/(Xb-Xf) 描述方程－》龙格库塔法求解位置－》求解我机角度
及位置－》判断距离时间－》end
计算机仿真技术
7
二阶系统：汽车轮子悬置系统
y’’(t)+p1y’(t)+y(t)=1.0
Y ( y1 , y2 , yn ) P ( p1, p2 , pk )
U (u1, u2 ,um )
G ( g1 , g2 , gn )
是 n 维状态向量 是 k 维参数向量 是 m 维控制向量 是n 维函数向量 计算机仿真技术
24
Y (0) [ y1 (0), y2 (0), yn (0)] 是给定的初始状态
计算机仿真技术
4
微分方程数学模型
一般形式
dny d n1 y dy d n1u d n2 u a1 n1 a n1 a n y c1 n1 c2 n1 cn u n dt dt dt dt dt
n为系统的阶次，ai为系统的结构参数，ci为输入函数的
16
连续系统数字仿真中离散化最基本的算法是数值积分 f ( y ,u ,t ) 的系统，已知系统变量 算法。对于形如 y 的初始条件 y (t 0 ) y 0 ，现在要求y随时间变化的过 程y(t)。计算过程可以这样考虑：首先求出初始点 y(t 0 ) y0 的 f (t 0 ，y0 ) ，微分方程可以写作：
26
今希望用直角坐标输出，故还引入两个定义变量 x y(21)，y)* cos[ y(2)] y(22) y(1)* sin[ y(2)]
根据卫星的发射速度，可以建立起方程组的初始条件： y(1)0=6400km(卫星到地心的距离，即为地球之半径)，
任务失败。
计算机仿真技术
6
Xb＝Vb*t
Yb=100
Xf’=Vf*cosθ
Yf’=Vf*sinθ
tanθ=(100-Yf)/(Xb-Xf) 描述方程－》龙格库塔法求解位置－》求解我机角度
及位置－》判断距离时间－》end
计算机仿真技术
7
二阶系统：汽车轮子悬置系统
y’’(t)+p1y’(t)+y(t)=1.0
Y ( y1 , y2 , yn ) P ( p1, p2 , pk )
U (u1, u2 ,um )
G ( g1 , g2 , gn )
是 n 维状态向量 是 k 维参数向量 是 m 维控制向量 是n 维函数向量 计算机仿真技术
24
Y (0) [ y1 (0), y2 (0), yn (0)] 是给定的初始状态
计算机仿真技术
4
微分方程数学模型
一般形式
dny d n1 y dy d n1u d n2 u a1 n1 a n1 a n y c1 n1 c2 n1 cn u n dt dt dt dt dt
n为系统的阶次，ai为系统的结构参数，ci为输入函数的
16
连续系统数字仿真中离散化最基本的算法是数值积分 f ( y ,u ,t ) 的系统，已知系统变量 算法。对于形如 y 的初始条件 y (t 0 ) y 0 ，现在要求y随时间变化的过 程y(t)。计算过程可以这样考虑：首先求出初始点 y(t 0 ) y0 的 f (t 0 ，y0 ) ，微分方程可以写作：
系统工程之连续系统仿真
b1s n 1 b2 s n 2 bn Y (s) G ( s) n n 1 U ( s ) s a1s an 1s an Y (s) X (s) X (s) U (s) (b1s
n 1
b2 s
n2
1 bn ) n s a1s n 1 an 1s an
Part 1 连续系统模型概述
系统的概念：所谓系统，是由相互联系、相互作 用的若干部分，以一定的结构组成的具有特定功 能的整体。 模型的概念：对系统的特征与变化规律的一种定 量抽象，是人们用以认识事物的一种手段。
•模型是现实系统的一种抽象，是在一定假设条件下对系统 的转化。 •模型中必须包含系统中的主要因素。 •模型中必须反映出主要因素之间的逻辑关系和数学关系。
系统仿真的应用
•对已有系统进行分析时，采用系统仿真时仿真成 为系统分析器。 •对尚未有的系统进行设计时，采用仿真技术考察 其性能是否满足预定要求，这时仿真是系统设计 器。 •在系统运行时，利用仿真模型作为观测器。 •在系统运行前，利用仿真模型作为预测器。 •利用仿真模型作为训练器，训练系统操作人员， 仿真成为训练仿真器。
系统
系统建模
模型
仿真试验
计算机
系统是研究的对象 模型是系统的抽象 仿真是对模型实验
仿真建模
计算机仿真的三个活动 • 系统建模是通过对实际系统的观测和检测，在忽略次要因 素及不可检测变量的基础上，用物理或数学的方法进行描 述，从而获得实际系统的简化近似模型。 • 仿真建模是将系统模型转化为仿真模型的过程，仿真模型 反映了系统模型同计算机之间的关系，它能为计算机所接 受并在其上运行。 • 仿真实验是将系统的仿真模型置于计算机上运行的过程。
系统工程导论 第五章 系统建模与仿真 第五节连续系统仿真与离散系统仿真
每天接待顾客人 次
发生天数
30～39 5
40～49 25
50～59 40
60～69 28
70～79以上 2
每位顾客购货金额 （元）
发生天数
10～19 20～29 30～39 40～49
40
30
15
10
50～59以 上
5
据此可以列出相应的概率分布，今若以随机数01、02、…、98、 99、00来表示上述概率分布
若没有事件发生，则仿真时钟继续等距推进；若有事件发生，则 记录此事件区间，从而可以得到有关事件的时间参数。
若有若干事件同时发生，除了记录该事件的时间参数外，还需事 先规定这种情况下对各类事件处理的优先顺序。
5.5 连续系统仿真与离散系统仿真
2）仿真建模策略
如何建立仿真模型，也就是采用何种方法推进仿真时钟，建立起系 统中各类实体之间的逻辑联系，一般有三种策略，这就是事件调度法、 活动扫描法和进程交互法。
随机数的生成方法大致有如下三种： 1）随机数表(random number table)法 即由人们在事先人为地产生
出一批均匀随机数，并制成表格形式备用。当需要使用它时，直接调 用这张随机数表就可以了。
2）随机数发生器 即在计算机上附加一个能产生随机数的装置，如 附加一个某种放射粒子的发射源装置，由于发射源在单位时间内发射 的粒子数量是随机的，所以计数器记录下来的数值就是随机数了。
① 事件调度法（Event Scheduling）。用事件的观点来分析真实系统， 通过定义事件及每个事件发生时系统状态的变化，按时间顺序确定并执 行每个事件发生时有关的逻辑关系，这就是事件调度法的基本思想。它 直接对事件加以调度。
所有事件连同其发生时间均放在事件表中。模型中有一个时间控制 模块，不断地从事件表中选择具有最早发生时间的事件，推进仿真时钟 到该事件发生的时间，并调用与该事件类型相应的事件处理模块，处理 完后在返回事件控制模块。
连续系统的数字仿真1_08解析
图3-2 义
梯形法的几何意
xk 1
xk
h[ f 2
(tk , xk )
f
(tk1, xk1)]
8
由于上式右边包含未知量xk+1,所以每一步都 必须通过迭代求解，每一步迭代的初值xk+1(0 通常采用欧拉公式来计算，因此梯形法每一步
迭代公式为
xk
(0) 1
xk
hf
(tk , xk )
xk
( 1
一般控制系统的输出动态响应在开始段变化较快， 到最后变化将会很缓慢。这时，计算可以采用变步长 的方法，即在开始阶段步长取得小一些，在最后阶段 取得大一些，这样即可以保证计算的精度，也可以加 快计算的速度。
23
对于一般工程计算，计算精度要求并不 太高，故常用定步长的方法。作为经验 数据，当采用四阶龙格－库塔法作数值 积分计算时，取计算步长
j 1
（3-8）
c1 0, m 1,2,, v
即：
k1 f (tk , xk )
k2 f (tk c2h, xk b21k1h)
k3 f (tk c3h, xk b31k1h b32k2h) 将式（3-7）展 开成h的幂级数并与微分方程式 （3-1）精确解式（3-6）逐项比较,便可求得式
t k 1
x(tk1) x(tk ) f (t, x(t))dt
tk
通常假设离散点t0,t1,…,tn是等距离的，即tk+1tk=h, 称h为计算步长或步距。
5
当t>t0时，x(t)是未知的，因此式（3-2） 右端的积分是求不出的。为了解决这个问题， 把积分间隔取得足够小，使得在tk与tk+1之间 的 f(t,x(t)) 可 以 近 似 看 作 常 数 f(tk,x(tk)) ， 这样便得到用矩形公式积分的近似公式
连续系统仿真概述
第一篇 连续系统仿真篇第一章 连续系统仿真概论按系统模型的特征分类，可以有连续系统仿真及离散事件系统仿真两大类。
本篇讨论连续系统仿真问题。
过程控制系统、调速系统、随动系统等这类系统称作连续系统，它们共同之处是系统状态变化在时间上是连续的，可以用方程式（常微分方程、偏微分方程、差分方程）描述系统模型。
1.1 连续系统模型描述连续系统仿真中的数学模型有很多种，但基本上可分为三类：连续时间模型、离散时间模型及连续－离散混合模型。
本节将对它们的线性定常形式作一介绍，1.2节将介绍几种模型结构变换的方法。
1.1.1 连续时间模型如果一个系统的输入量u(t)，输出量y(t)，系统的内部状态变量x(t)都是时间的连续函数，那么我们可以用连续时间模型来描述它。
系统的连续时间模型通常可以有以下几种表示方式：常微分方程，传递函数，权函数，状态空间描述．本节仅对其一般描述形式作一简要介绍，偏微分方程将在第8章讨论。
1.常微分方程常微分方程可用(1.1)式表示：u c dt ud c dt u d c y a dt dy a dt y d a dt y d n n n n n n n n n nn +++=++++------- 1221111111(1.1) 其中n 为系统的阶次，a i n i (,,,,)=012 为系统的结构参数，),,2,1(n j c j =为输入函数的结构参数，它们均为实常数。
2.传递函数若系统的初始条件为零，即系统在 t=0时已处于一个稳定状态，那么对（1.1）式两边取拉氏变换后可得：)()()()()()()(2211111s U c s U sc s U sc s Y a s sY a s Y s a s Y s n n n n n n n +++=++++----稍加整理，并记：jn j jn n j jjn s a s cs U s Y s G ∑∑=--=-==10)()()( (1.2)(1.2)式称为系统的传递函数。
实验二 面向系统结构图的连续系统数字仿真实验 matlab程序
4.离散相似法编程
时域离散编程求解 程序：
clear num=[8 10]; den=[0.1 1 0 0]; [A,B,C,D]=tf2ss(num,den) sysc=ss(A,B,C,D); T=0.05; sysd=c2d(sysc,T); Ad=sysd.a; Bd=sysd.b;
0 0.5 1
程序：
clear all A = [-1000.25 999.75 0.5;999.75 -1000.25 0.5;0 0 0]; t = 0; x = [1 -1 1]; h = 0.05; M = round(10/h); for k = 1:M tt = t(k) + h; x1 = -exp(-0.5*tt)+exp(-2000*tt)+1; x2 = -exp(-0.5*tt)-exp(-2000*tt)+1; x3 = 0; temp = [x1 x2 x3]; x = [x; temp]; t = [t; tt]; end T = 0.003125; m = 4; t_frog = 0; x_frog = [1 -1 1]; start_time_frog = clock; F = expm(A*T); for k = 0:m tt = 2^(k) * T; temp = F*x_frog(1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; if k<m F=F*F; end end qT = 2^(m)*T; for k = 1:M-1 tt = t_frog(m+k+1) + qT; temp = F*x_frog(m+k+1,:)'; temp = temp'; x_frog = [x_frog; temp]; t_frog = [t_frog; tt]; end pass_time_frog = etime(clock,start_time_frog); hh = 0.001; end pass_time_rk4=etime(clock,start_time_rk4); t_display = t(2:31,:); x1_display = x(2:31,1); x1_frog_display = x_frog(6:35,1); x1_rk4_display = x_rk4(51:50:1501,1); disp(' 时间 解析解 蛙跳法 RK4法 ') t_rk4 = 0; x_rk4 = [1 -1 1]; start_time_rk4=clock; for k = 1:50*M tt = k*hh; K1 = A*x_rk4(k,:)'; K2 = A*(x_rk4(k,:)'+hh*K1/2); K3 = A*(x_rk4(k,:)'+hh*K2/2); K4 = A*(x_rk4(k,:)'+hh*K3); temp = x_rk4(k,:)'+hh*(K1/2+2*K2+2*K3+K4)/6; temp = temp'; x_rk4 = [x_rk4; temp]; t_rk4 = [t_rk4; tt];
连续系统数字仿真数值积分法
main() { float x11,x21,x31,x10,x20,x30,R; int i; x10=0;x20=0;x30=0;outputy[0][0]=0; ST=10;DT=0.005; LP=ST/DT;VN=1;R=20; for (i=1;i<=LP;i++) { x11=x10-DT*x30+DT*R; x21=DT/0.5*x10+(1-DT/0.5)*x20; x31=DT/0.1*x20+(1-DT/0.1)*x30; outputy[0][i]=x31; x10=x11;x20=x21;x30=x31; } dispcurve(); } 程序为llx56.c
4．2 梯形法 为了提高仿真精度，离散-再现环节采取图4．3 的形式。 (4-11)
X (k 1) X (k )
T [e(k ) e(k 1)] 2
(4-12)
式(4-12)称为梯形公式，其几何解释如图4．4所示。
[例4．1] 已知一多变量系统的结构框图如图4．5所示，请 用梯形公式对此系统进行仿真，并输出y1、y2的仿真结果。
假设线性系统的状态空间描述为
AX BU X Y CX DU
式中：X为n×1维状态向量；U为r×l维输入 向量；A为n×n维系统矩阵；B为n×r维输入 矩阵；Y为m×1维输出向量；C为m×n维输出 矩阵；D为m×r维传递矩阵。
(4-3) (4-4) (4-5) (4-6)
4.3 龙格-库塔(Runge-Kutta)法
在非实时仿真中，有时需要更高的精度。 4．3节和4．4节中将再介绍两种更精确的方 法及其离散-再现环节。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第3章 连续系统的数字仿真
本章内容
(1) 熟悉在数字计算机仿真技术中常用的几种数值积分法 ，
特别是四阶龙格-库塔法； (2) 典型环节及其系数矩阵的确定； (3) 各连接矩阵的确定； (4) 利用MATLAB在四阶龙格-库塔法的基础上，对以状 态
空间表达式和方框图描述的连续系统进行仿真； (5) 了解以增广矩阵法为基础的连续系统的快速仿真方法1 。
xk 1
xk
hxk
1 2!
h2
xk
hp p!
xk( p)
0(h p31-)5）
可以看出，提高截断误差的阶次，便可提高其精
度，但是由于计算各阶导数相当麻烦，所以直接采
用泰勒级数公式是不适用的，为了解决提高精度问
Hale Waihona Puke 题，龙格和库塔两人先后提出了间接使用泰勒级数
公式的方法，即用函数值f (t,x)的线性组合来代替f
j 1
（3-8）
c1 0, m 1,2,, v
即：
k1 f (tk , xk )
k2 f (tk c2h, xk b21k1h)
k3 f (tk c3h, xk b31k1h b32k2h) 将式（3-7）展 开成h的幂级数并与微分方程式 （3-1）精确解式（3-6）逐项比较,便可求得式
t k 1
x(tk1) x(tk ) f (t, x(t))dt
tk
通常假设离散点t0,t1,…,tn是等距离的，即tk+1tk=h, 称h为计算步长或步距。
5
当t>t0时，x(t)是未知的，因此式（3-2） 右端的积分是求不出的。为了解决这个问题， 把积分间隔取得足够小，使得在tk与tk+1之间 的 f(t,x(t)) 可 以 近 似 看 作 常 数 f(tk,x(tk)) ， 这样便得到用矩形公式积分的近似公式
干个离散点a≤ t0 < t1<… <tn≤b处的近似值
x(t1),x(t2),…,x(tn)。
4
3.1.1 欧拉法
欧拉法又称折线法或矩形法，是最简单也是最 早的一种数值方法
将式（3-1）中的微分方程两边进行积分，得
tk1 dx(t)
t k 1
dt f (t, x(t))dt
tk dt
tk
即
图3-2 义
梯形法的几何意
xk 1
xk
h[ f 2
(tk , xk )
f
(tk1, xk1)]
8
由于上式右边包含未知量xk+1,所以每一步都 必须通过迭代求解，每一步迭代的初值xk+1(0 通常采用欧拉公式来计算，因此梯形法每一步
迭代公式为
xk
(0) 1
xk
hf
(tk , xk )
xk
( 1
或简化为 x(tk1) x(tk ) f (tk , x(tk ))h
这就是欧拉公x式k。1 xk f (tk , xk )h
6
以x(t0)=x0作为初始值，应用欧拉公式， 就可以一步步求出每一时刻tk的xk值，
即 k =0,x1≈x0+f(t0,x0)h
k=1, x2≈x1+f(t1,x1)h x2
┇
x1
k=n-1,
x
(t2 ,x2)
(t1 ,x1)
xn≈xn-1+f(tn-1,xn-1)h 这样式(3-1)的解x(t) 就求出来了。欧拉法的计 算虽然比较简单，但精度 较低。图3-1为欧拉法的 几何解释。
x0
t
t0 t1 t2 图3-1 欧拉法的几何解释
7
3.1.2 梯形法
由上可知欧拉公式中的积分是用矩形面积
xk 1
xk
hfk
h2 2!
(
f
' k
f' xk
fk ) 0(h p1（) 3-6）
为了避免计算式（3-6）中的各阶导数项，可令xk+1由
以下多项式表示。
v
xk1 xk h amkm m1
（3-7）
12
式中 am为待定因子，v为使用f函数值的个数，
km满足下列方程
m1
km f (tk cmh, xk bmjk jh)
hf
(tk , xk )
xk
(1) 1
xk
h[ f 2
(tk , xk )
f
(tk
1
,
xk
(0) 1
)]
通常称这类方法为预估－校正方法。它首
先根据欧拉公式计算出xk+1的预估值
xk+1(0) ，然后再对它进行校正，以得到更
准确的近似值xk+1(1)。
10
3.1.4 龙格－库塔法
根据泰勒级数将式（3-1）在tk+1=tk+h时刻的解 xk+1=x(tk+h) 在tk附近展开，有
（3-7）和式（3-8）中的系数am ,bmj和cm等。
13
现以v=2为例，来说明这些参数的确定方法。
用数字计算机来仿真或模拟一个连续控 制系统的目的就是求解系统的数学模型。 由控制理论知，一个n阶连续系统可以被描 述成由n个积分器组成的模拟结构图。因此 利用数字计算机来进行连续系统的仿真， 从本质上讲就是要在数字计算机上构造出n 个数字积分器，也就是让数字计算机进行n 次数值积分运算。可见，连续系统数字仿 真中的最基本的算法是数值积分算法。
2
3
3.1 数值积分法
连续系统通常把数学模型化为状态空间表达式，为
了对n阶连续系统在数字计算机上仿真及求解，就要 采用数值积分法来求解系统数学模型中的n个一阶微 分方程。
设n阶连续系统由以下n个一阶微分方程组成
dx(t)
dt
f (t, x(t))
（3-1）
x(t0 ) x0
所谓数值积分法，就是要逐个求出区间［a ,b］内若
(t,x)的导数，然后按泰勒公式确定其中的系数, 这
样既能避免计算f (t,x)的导数，又可以提高数值计算
精度，其方法如下。
11
因
xk f (tk , xk ) fk
xk
df
(t, dt
x)
t tk
f f dx
( t
x
) dt
t tk
fk f xk fk
x xk
x xk
故式（3-5）可写成
R1)
xk
h[ 2
f
(tk , xk )
f
(tk
1
,
x(R) k 1
)]
(3-3）
式中 迭代次数R=0,1,2,…
9
3.1.3 预估－校正法
虽然梯形法比欧拉法精确，但是由于每 一步都要进行多次叠代，计算量大，为了 简化计算，有时只对式（3-3）进行一次叠 代就可以了，因此可得
xk
(0) 1
xk
f(tk,xk)h 来近似的。 由图3-2知，用矩形面积
f dx
dt
c
tkabtk+1代替积分，其误差就
是图中阴影部分。为了提高精
a
b
度现用梯形面积tkactk+1来代
替积分，即
tk1
h
tk
t k 1
t
tk
f
(t, x(t))dt
[ 2
f
(tk , xk )
f
(tk1, xk1)]
于是可得梯形法的计算公式为
本章内容
(1) 熟悉在数字计算机仿真技术中常用的几种数值积分法 ，
特别是四阶龙格-库塔法； (2) 典型环节及其系数矩阵的确定； (3) 各连接矩阵的确定； (4) 利用MATLAB在四阶龙格-库塔法的基础上，对以状 态
空间表达式和方框图描述的连续系统进行仿真； (5) 了解以增广矩阵法为基础的连续系统的快速仿真方法1 。
xk 1
xk
hxk
1 2!
h2
xk
hp p!
xk( p)
0(h p31-)5）
可以看出，提高截断误差的阶次，便可提高其精
度，但是由于计算各阶导数相当麻烦，所以直接采
用泰勒级数公式是不适用的，为了解决提高精度问
Hale Waihona Puke 题，龙格和库塔两人先后提出了间接使用泰勒级数
公式的方法，即用函数值f (t,x)的线性组合来代替f
j 1
（3-8）
c1 0, m 1,2,, v
即：
k1 f (tk , xk )
k2 f (tk c2h, xk b21k1h)
k3 f (tk c3h, xk b31k1h b32k2h) 将式（3-7）展 开成h的幂级数并与微分方程式 （3-1）精确解式（3-6）逐项比较,便可求得式
t k 1
x(tk1) x(tk ) f (t, x(t))dt
tk
通常假设离散点t0,t1,…,tn是等距离的，即tk+1tk=h, 称h为计算步长或步距。
5
当t>t0时，x(t)是未知的，因此式（3-2） 右端的积分是求不出的。为了解决这个问题， 把积分间隔取得足够小，使得在tk与tk+1之间 的 f(t,x(t)) 可 以 近 似 看 作 常 数 f(tk,x(tk)) ， 这样便得到用矩形公式积分的近似公式
干个离散点a≤ t0 < t1<… <tn≤b处的近似值
x(t1),x(t2),…,x(tn)。
4
3.1.1 欧拉法
欧拉法又称折线法或矩形法，是最简单也是最 早的一种数值方法
将式（3-1）中的微分方程两边进行积分，得
tk1 dx(t)
t k 1
dt f (t, x(t))dt
tk dt
tk
即
图3-2 义
梯形法的几何意
xk 1
xk
h[ f 2
(tk , xk )
f
(tk1, xk1)]
8
由于上式右边包含未知量xk+1,所以每一步都 必须通过迭代求解，每一步迭代的初值xk+1(0 通常采用欧拉公式来计算，因此梯形法每一步
迭代公式为
xk
(0) 1
xk
hf
(tk , xk )
xk
( 1
或简化为 x(tk1) x(tk ) f (tk , x(tk ))h
这就是欧拉公x式k。1 xk f (tk , xk )h
6
以x(t0)=x0作为初始值，应用欧拉公式， 就可以一步步求出每一时刻tk的xk值，
即 k =0,x1≈x0+f(t0,x0)h
k=1, x2≈x1+f(t1,x1)h x2
┇
x1
k=n-1,
x
(t2 ,x2)
(t1 ,x1)
xn≈xn-1+f(tn-1,xn-1)h 这样式(3-1)的解x(t) 就求出来了。欧拉法的计 算虽然比较简单，但精度 较低。图3-1为欧拉法的 几何解释。
x0
t
t0 t1 t2 图3-1 欧拉法的几何解释
7
3.1.2 梯形法
由上可知欧拉公式中的积分是用矩形面积
xk 1
xk
hfk
h2 2!
(
f
' k
f' xk
fk ) 0(h p1（) 3-6）
为了避免计算式（3-6）中的各阶导数项，可令xk+1由
以下多项式表示。
v
xk1 xk h amkm m1
（3-7）
12
式中 am为待定因子，v为使用f函数值的个数，
km满足下列方程
m1
km f (tk cmh, xk bmjk jh)
hf
(tk , xk )
xk
(1) 1
xk
h[ f 2
(tk , xk )
f
(tk
1
,
xk
(0) 1
)]
通常称这类方法为预估－校正方法。它首
先根据欧拉公式计算出xk+1的预估值
xk+1(0) ，然后再对它进行校正，以得到更
准确的近似值xk+1(1)。
10
3.1.4 龙格－库塔法
根据泰勒级数将式（3-1）在tk+1=tk+h时刻的解 xk+1=x(tk+h) 在tk附近展开，有
（3-7）和式（3-8）中的系数am ,bmj和cm等。
13
现以v=2为例，来说明这些参数的确定方法。
用数字计算机来仿真或模拟一个连续控 制系统的目的就是求解系统的数学模型。 由控制理论知，一个n阶连续系统可以被描 述成由n个积分器组成的模拟结构图。因此 利用数字计算机来进行连续系统的仿真， 从本质上讲就是要在数字计算机上构造出n 个数字积分器，也就是让数字计算机进行n 次数值积分运算。可见，连续系统数字仿 真中的最基本的算法是数值积分算法。
2
3
3.1 数值积分法
连续系统通常把数学模型化为状态空间表达式，为
了对n阶连续系统在数字计算机上仿真及求解，就要 采用数值积分法来求解系统数学模型中的n个一阶微 分方程。
设n阶连续系统由以下n个一阶微分方程组成
dx(t)
dt
f (t, x(t))
（3-1）
x(t0 ) x0
所谓数值积分法，就是要逐个求出区间［a ,b］内若
(t,x)的导数，然后按泰勒公式确定其中的系数, 这
样既能避免计算f (t,x)的导数，又可以提高数值计算
精度，其方法如下。
11
因
xk f (tk , xk ) fk
xk
df
(t, dt
x)
t tk
f f dx
( t
x
) dt
t tk
fk f xk fk
x xk
x xk
故式（3-5）可写成
R1)
xk
h[ 2
f
(tk , xk )
f
(tk
1
,
x(R) k 1
)]
(3-3）
式中 迭代次数R=0,1,2,…
9
3.1.3 预估－校正法
虽然梯形法比欧拉法精确，但是由于每 一步都要进行多次叠代，计算量大，为了 简化计算，有时只对式（3-3）进行一次叠 代就可以了，因此可得
xk
(0) 1
xk
f(tk,xk)h 来近似的。 由图3-2知，用矩形面积
f dx
dt
c
tkabtk+1代替积分，其误差就
是图中阴影部分。为了提高精
a
b
度现用梯形面积tkactk+1来代
替积分，即
tk1
h
tk
t k 1
t
tk
f
(t, x(t))dt
[ 2
f
(tk , xk )
f
(tk1, xk1)]
于是可得梯形法的计算公式为