第2章 连续变量动态系统数字仿真算法

航空航天中的计算方法
授课教师：陈琪锋 中南大学航空航天学院
第2章 连续变量动态系统 数字仿真算法
内容提要 2.1 2.2 2.3 2.4 2.5 2.6 2.7 连续变量系统仿真概述 常微分方程数值算法原理 常微分方程数值求解的常用算法 常微分方程数值求解方法的稳定性 常微分方程数值求解方法的选择原则 刚性系统数值求解方法 含间断特性的系统仿真
t2
t
PageBiblioteka Baidu13
2.2 常微分方程数值算法原理 欧拉法的几何解释：折线近似，矩形求积 ∫t1t0f［t,y(t)］dt f(t0,y0)(t1-t0) f(t1,y1) f(t0,y0)
f(t)
以矩形面积 近似
t
t t0
2016/9/25
t1
航空航天中的计算方法 Page 14
2.2 常微分方程数值算法原理 欧拉法的特点： 方法简单，计算量小 计算精度较低 单步法，可自启动
如果微分方程中的控制参量全部给定为常量，该微分方 程就是定常系数微分方程，否则为变系数微分方程。
2016/9/25
航空航天中的计算方法
Page 7
2.1 连续变量系统仿真概述 随机的微分方程或确定的微分方程 在微分方程的全部控制参量和输入变量中，只要有一个 是随机变量，该微分方程就称为随机的微分方程，否则 为确定的微分方程。它们所描述的系统分别称为随机系 统或确定系统。
2016/9/25 航空航天中的计算方法 Page 6
2.1 连续变量系统仿真概述 定常系数微分方程或时变系数微分方程 微分方程中的变量：
状态变量刻划系统的演化规律，它在系统的状态空间变 化； 控制参量表示环境对系统的规定，它在系统的参量空间 变化； 输入变量表示环境对系统的作用，分为可控输入和不可 控输入，不可控输入又称为干扰输入。
f(t)
以中值点矩 形面积近似
t
t t0
2016/9/25
t1
航空航天中的计算方法
t2
Page 20
2.2 常微分方程数值算法原理 中值点法的特点： 方法简单，计算量小 计算精度较低 多步法，不可自启动，需借助其它算法启动
例1计算
2016/9/25
航空航天中的计算方法
Page 21
2.2 常微分方程数值算法原理 2.2.3 梯形法 常微分方程：
中值点法
y0
t
y(t)
dy/dt＝f(t,y)
y(t0)＝y0
y(tm1 ) y(tm 1 )
tm1 tm1
f (t , y)dt
t0
2016/9/25
t1
航空航天中的计算方法
t2
t
Page 19
2.2 常微分方程数值算法原理 中值点法的几何解释：按中值点矩形求积 ∫t2t0f［t,y(t)］dt f(t1,y1)(t2-t0) f(t1,y1) f(t0,y0)
导数： 小步长近似：
一般公式：
t m t 0 m t ,
ym y( t m ), m 1
2016/9/25
航空航天中的计算方法
Page 18
2.2 常微分方程数值算法原理 中值点法的几何解释：按中值点矩形求积 y1 y2 y2＝y0＋f(t1, y1 )2t 欧拉法
y(t1)
本课程主要解决第二个问题
2016/9/25
航空航天中的计算方法
Page 11
2.2 常微分方程数值算法原理 2.2 常微分方程数值算法原理 2.2.1 欧拉法 dy 常微分方程： f (t , y ) dt y ( t 0 ) y0
dy t O ( t 2 ) Taylor展开： y( t0 t ) y( t0 ) dt t0

v(t0 ) v0
MATLAB应用简介 环境界面 矩阵运算
语言
• • • • • • 循环（for、while） 条件if 开关switch 函数 画图plot、plot3、mesh 符号运算
2016/9/25
航空航天中的计算方法
Page 17
2.2 常微分方程数值算法原理 2.2.2 中值点法 常微分方程：
2016/9/25 航空航天中的计算方法 Page 22
2.2 常微分方程数值算法原理 梯形法的几何解释：按梯形求积 ∫t1t0f(t)dt 1/2*[f(t1,y1)+ f(t0,y0)]*(t1-t0) f(t1,y1) f(t0,y0) 以梯形面积 近似
f(t)
t
t t0
2016/9/25

v(t0 ) v0
航空航天中的计算方法 Page 5
2.1 连续变量系统仿真概述 线性微分方程或非线性微分方程 线性系统的基本特征——满足叠加原理（superposition principle）。 令f代表某种操作，如关系、变换、运算、函数、泛函、 方程等，x记操作对象（如变量），f(x)表示对x施加f操作 的结果。 如果f(x)满足条件： (1) 加性 f(x1+x2) = f(x1) + f(x2) (2) 奇性 f(k*x)= k*f(x)， k为常数 就称f为线性操作。加性与奇性合并表示为 f(a*x1+b*x2) = a*f(x1) +b*f(x2) ，a、b为常数 称为叠加原理。
利用多个点上的函数值 线性组合提高精度，同 时避免计算高阶导数
dy 1 d2y 2 3 t t O ( t ) Taylor展开： y( t0 t ) y( t0 ) 2 dt t0 2 dt t
0
多点斜率线性组合（两点为例）：
y( t 0 t ) y0 a1k1 a2 k2 t k1 f ( t 0 , y0 ) k2 f ( t 0 t1 , y0 y1 )
2.1 连续变量系统仿真概述 2.1 连续变量系统仿真概述 连续变量动态系统（CVDS）指系统状态随时间连续变化 的系统。 CVDS一般是遵守物理学定律或广义物理学定律的自然系 统或社会系统。通常以微分方程的形式来描述。 2.1.1 微分方程分类 常微分方程或偏微分方程 如果状态变量只是时间的函数，与空间分布无关，称为 集中参数系统，用常微分方程描述； 如果状态变量同时依赖于时间和空间分布（同一时刻不 同点的状态变量数值不同），称为分布参数系统，须用 偏微分方程来描述。
t1
航空航天中的计算方法 Page 23
2.2 常微分方程数值算法原理
1 f ( t0 , y0 ) f ( t1 , y1 ) t1 t 0 梯形求积近似： y1 y0 2
隐式方法
求y1: 预估-校正法： 解非线性方程，迭代？
P ym 1 ym f ( t m , y m )t
1.69 105 kg ( m s )
r 0.01m g 9.8 m s 2
24 6 C D ( Re ) 0.4 Re 1 Re Re 2 vr
初始条件： v0 0 m s 仿真设定： T 25 s h t 0.25 s
航空航天中的计算方法 Page 16
i 1
r
k i f ( t i , yi )
航空航天中的计算方法
2016/9/25
Page 25
2.3 常微分方程数值求解的常用算法 2.3 常微分方程数值求解的常用算法 2.3.1 龙格-库塔（Runge-Kutta）法 dy 常微分方程： f (t , y ) dt y ( t 0 ) y0
2.2 常微分方程数值算法原理 欧拉法的几何解释：折线近似，矩形求积 y1 y1＝y0＋f(t0, y0)t
y(t1)
y(t)
y0
t
dy/dt＝f(t,y) y(t0)＝y0
y(tm1 ) y(tm )
tm1 tm
f (t , y)dt
t0
2016/9/25
t1
航空航天中的计算方法
2.2 常微分方程数值算法原理
梯形法思想的推广： 多点斜率加权平均
1 1 ym 1 ym k1 k2 t 2 2 k1 f ( t m , ym )
P k 2 f ( t m 1 , ym 1 ) f ( t m 1 , ym 1 )
ym 1 ym wi ki t
t P ym 1 ym f ( t , y ) f ( t , y m m m 1 m 1 ) 2 tm t0 mt , ym y(t m ), m 0,1, 2,
欧拉法预估，梯形法校正，仅迭代一次
2016/9/25 航空航天中的计算方法 Page 24
参考资料：
[1] 黄柯棣，系统仿真技术，国防科技大学出版社，1998. [2] David L. Darmofal, Computational Methods in Aerospace Engineering (Lecture Notes), MIT, 2005.
2016/9/25 航空航天中的计算方法 Page 3
小步长近似： y(t0 t ) y0 f (t0 , y0 )t 一般公式：
ym 1 ym f ( t m , y m )t
t m t 0 m t ,
ym y( t m ), m 0,1, 2,
Page 12
2016/9/25
航空航天中的计算方法
2016/9/25
航空航天中的计算方法
Page 8
2.1 连续变量系统仿真概述 自由系统与强迫系统 若状态变量的变化中没有外作用项，称为自由系统： X’ = F(t, X; P) + u(t) 若演化方程中包含有外作用项u(t)（亦称强迫作用项，控 制作用项），如以下方程描述的系统，称为强迫系统： X’ = F(t, X; P) + u(t) 如果引入新的状态变量xn+1(t)，令x’n+1 = u(t)，强迫系统 (1.3-4)可以转化为n+1维自由系统来处理。
航空航天中的计算方法
Page 10
2.1 连续变量系统仿真概述 2.1.2 连续变量动态系统仿真的问题 如何对给定的CVDS建立数学模型——微分方程
如何把数学模型转换成计算机可接受的、等价的仿真模 型，并编制程序形成可执行的指令序列——运行模型
如何对所建模型进行校验、确认和对建模仿真应用进行 鉴定验收 如何在计算机环境下利用其运行模型作实验 如何分析评估仿真结果
dy f (t , y ) dt y ( t 0 ) y0
积分法：
y(t1 ) y0 f (t , y )dt
t0
t1
1 f ( t0 , y0 ) f ( t1 , y1 ) t1 t 0 梯形求积近似： y1 y0 2
线性近似，t0和t1处斜率的平均值 欧拉法？ 中值点法？
2016/9/25
航空航天中的计算方法
Page 9
2.1 连续变量系统仿真概述 自治系统与非自治系统 演化方程中不显含时间t的系统是自治系统： X’ = F(X; P) 否则，是非自治系统： X’ = F(t, X; P)
时变系数系统和含有作为时间函数的强迫作用u(t)的系统 都是非自治系统。
2016/9/25
dy f (t , y ) dt y ( t 0 ) y0
dy y ( t t ) y ( t t ) lim dt t 0 2 t
y(t t ) y(t t ) 2 f (t , y )t
y m 1 y m 1 2 f ( t m , y m ) t
2016/9/25
航空航天中的计算方法
Page 15
2.1 连续变量系统仿真概述 例1：常微分方程描述的系统举例——自由下落球体
dv 1 g D( v ) dt m
1 D(v ) r 2 v 2C D ( Re ) 2
参数： m 0.0038kg 0.9 kg m 3
2016/9/25 航空航天中的计算方法 Page 4
2.1 连续变量系统仿真概述 例1：常微分方程描述的系统举例——自由下落球体
dv m mg D(v ) dt
空气阻力
D( v )
1 r 2v 2C D ( Re ) 2
阻力系数
雷诺数 初始条件
2016/9/25
24 6 C D ( Re ) 0.4 Re 1 Re Re 2 vr