3第三章 数值积分法在系统仿真中的应用

对式(3.1)两端由 t o到t1进行积分 , 得到
y (t1 ) = y 0 + ∫ f ( y , t )dt
t0 t1
（3-3） ）
0
t0
t1
t
图3.1 矩形近似及其误差
第三章 数值积分法在系统仿真中的应用
欧拉法的几何意义
欧拉法的几何意义 十分清楚。
y
y0
y1
y2
图3.2 通过 ( t 0 , y 0 )点作为积分曲线的切线 ，其 斜率为 f ( y 0 , t 0 )，此切线与 t1处平行于 y轴直 线的交点即为 y1，再过 ( t1 , y1 )点作积分曲线 的切线 , 它与过 t 2 平行于 y轴直线的交点为 y 2。 这样过 ( t 0 , y 0 )，1 , y1 )， 2 , y 2 )， ，得到一 (t (t ⋯
其中∇f n 为向前插分算子，定义为
∇ 0 f n +1 = f n +1 ∇1 f n +1 = f n +1 − f n ⋯⋯ ∇ i f n +1 = ∇(∇ i −1 f n +1 ) = ∇ i −1 f n +1 − ∇ i −1 f n

（3-32） ）
第三章 数值积分法在系统仿真中的应用
t0
t1 t 2
t
称为欧拉折线法。 条折线，
图3.2 欧拉折线
第三章 数值积分法在系统仿真中的应用
改进的欧拉法
在推导时用图中的阴影面积来近似 式(3.3)时，由于梯形公式中隐含有待求 时 通常可用欧拉法启动初值， 量，通常可用欧拉法启动初值，算出近 似值，然后带如微分方程，最后利用梯 似值，然后带如微分方程， 形公式求出修正。为提高精度， 形公式求出修正。为提高精度，简化计 只迭代一次。这样可得改进的欧拉 算，只迭代一次。这样可得 公式：
∇f n = f n − f n −1 ∇ 2 f n = ∇( f n − 2 f n −1 ) + f n − 2 ⋯⋯ ∇ k f n = ∇ k −1 f n − ∇ k −1 f n −1
∇0 fn = fn
（3-26） ）
第三章 数值积分法在系统仿真中的应用
f (t )
ynp+1 = yn + h ⋅ f ( yn , tn ) （3-8） ） yn+1 = yn + h [ f ( yn , tn ) + f p ( ynp+1, tn+1)] 2
0
t0
t1
t
第一式称为预估公式， 第一式称为预估公式， 第二式称为校正公式。 第二式称为校正公式。
i =0 i = −1
m
m
（3-2） ）
第三章 数值积分法在系统仿真中的应用
几种常用的积分法
欧拉法 欧拉法的几何意义 改进的欧拉法 龙格-库塔法 亚当斯法(显式) 亚当斯法(隐式)
第三章 数值积分法在系统仿真中的应用
欧拉法
f (t )
误差
欧拉法虽然计算精度较低，实际中很少采用， 虽然计算精度较低，实际中很少采用， 虽然计算精度较低 但器推倒简单， 但器推倒简单，能说明旧够数值解法一般计算公 式的基本思想。 式的基本思想。
yn +1 = yn + hk 2 k1 = f (t n , yn ) k 2 = f (t n + h , yn + h k1) 2 2
图3.3 梯形近似及其误差
第三章 数值积分法在系统仿真中的应用
龙格-库塔法
将式 (3 .3)在 t n点展开 泰勒级数
h2 ɺ ɺɺ(t n ) + o(h 3 ) y (t n + h) = y (t n ) + hy (t n ) + y 2
龙格-库塔（RK）法的一般形式为 （3-9） ）
yn +1 = yn + h∑ wi ki
第三章 数值积分法在系统仿真中的应用
3.3 实时仿真法
假设仿真的连续动力学由非线形常微分方程描述为： 假设仿真的连续动力学由非线形常微分方程描述为：
dy = f ( y, u, t ) dt
y (t0 ) = y0
（3-37） ）
对(3-37)式采用二阶龙格 库塔公式求解，其递推方程可写为 式采用二阶龙格 库塔公式求解 式采用二阶龙格-库塔公式求解，
0
（k=1时可得欧拉公式）
1
第三章 数值积分法在系统仿真中的应用
当k=2时，得到亚当斯多步法的计算公式，(3-28)式各系数为
v0 = 1
v1 = ∫ f n ds = 1 0 2
1
将v 0 , v1代入 (3 − 27 )式得
y (t n +1 ) = y (t n ) + 1 h(3 f n − f n −1 ) 2
i =1
r
（3-10） ）
式中
k1 = f ( y n , t n )
k i = f ( y n + h∑ β ij k j , t n + α i h)
j =1 i −1
(i = 2, 3, ⋯ r )
第三章 数值积分法在系统仿真中的应用
而
α
i
=
∑
i −1 j =1
β ij
ai , β ij , wi为待定系数 , r为使用 k的个数 (即阶数 )。在给定 r值后，通 过把式 (3 .10 )展开成 h 的幂级数，然后和台劳 展开式的系数进行对比 ， 以确定 β ij , wi的值。
yn +1 = yn + 1 (k1 + k 2 ) 2 k1 = hf (t n , yn ) k 2 = hf (t n + h, yn + k1 )

（3-38） ）
F为函数，外部输入为u(t) 。 为函数，外部输入为 为函数
第三章 数值积分法在系统仿真中的应用
采入un +1
第三章 数值积分法在系统仿真中的应用
3.1 连续系统仿真中常用的数值积分法
1. 数值积分法
如果已知某一系统的一阶向量微分方程为
ɺ y = f ( y, t ), y (t0 ) = y0
对式子(3.1),数值积分可写成统一的公式
（3-1） ）
yn +1 = ∑ ai yn −i + h ∑ β i f n −i
（3-34） ）
Baidu Nhomakorabea
c 校正： y n+1 = y n + (9 fn+1 + 19 fn − 5 fn−1 + fn−2 ) （3-35） ）
第三章 数值积分法在系统仿真中的应用
3.2 刚性系统的特点及算法
一个刚性系统可以这样描述，对于n阶微分方程组 可以这样描述，对于 阶微分方程组 一个 可以这样描述
整理上式得
y (t n +1 ) = y (t n ) + 1 h(23 f n − 16 f n −1 + 5 f n − 2 ) 12
（3-30） ）
第三章 数值积分法在系统仿真中的应用
亚当斯法(隐式)
牛顿前插公式为
Pk , n ( t ) = f n +1 +
( t − t n +1 ) ∇ f n +1 + ⋯ + h （3-31） ） ( t − t n −1 )( t − t n ) ⋯ ( t − t n + 2 − k ) k ∇ f n +1 k h k
dy = f (t , y ) dt
y (a) = η
（3-36） ）
其中 y = [ y1 ⋯ yn ]T , f = [ f1 ⋯ f n ]T ,η = [η1 ⋯ η n ]T 为n维向量 其雅可比矩阵 ∂f 的特征值实部表示系统的衰减速度，其最大特 ∂y
征值和最小特征值实部之比，即
ρ = max Re λk / min Re λk
亚当斯多步法的计算公式是
y (t n +1 ) = y (t n ) + h∑ vi ∇ i f n
i =0
k −1
（3-27） ）
其中
1 1 s ( s + 1) ⋯ ( s + i − 1)ds i ≥ 1 v1 = ∫ it 0
v0 = 1
（3-28） ）
y (t n +1 ) = y (t n ) + h ∫ f n ds = y (t n ) + hf n
1≤ k ≤ n 1≤ k ≤ n
作为系统刚性程序的度量。
第三章 数值积分法在系统仿真中的应用
当 ρ >> 1 时，系统为刚性系统，或称为stiff系统。对与这样的系统作做 数字仿真，其最大的困惑是：积分步长由最大的特征值来确定，最小的
特征值决定数值求解总的时间。
刚性系统在时间中的普遍性和重要性已得到广泛的重视，这种方程的数 值解已成为常微分方程的数值研究的重点。
目前解刚性方程的数值方法基本分为： 显式公式 隐式公式 预测校正
第三章 数值积分法在系统仿真中的应用
显式公式常用雷纳尔法。其中着眼点是，在保证稳定的前提下，尽 可能地扩大稳定区域。这一方法的优点是，它是显式的，所以便于程序 设计。对一般好的方程设计。对一般条件好的方程，它就还原为四阶龙 格-库塔方法，而对刚性方程它又有增加稳定性的好处。 众所周知，隐式公式都是稳定的，故都大于解描述刚性系统的方程 组，如隐式的龙格-库塔法。但这种方法每计算一步都要进行迭代，故计 算量大。在工程上使用又一定捆年。因此在解刚性方程时，常Rosenbrock 提出的半隐式龙格-库塔法。 预测-校正型中常用的解刚性方程的方法式Gear算法。Gear首先应 引进刚性稳定性的概念，它可以满足稳定型，而减低对h的要求。Gear 方法是一格通用的方法，它不但使用于解刚性方程组，而且也适用于解 非刚性方程组。
4阶龙格-库塔法式使用较多的一种方法，其公式如下
y n +1 = y n + h ( k1 + 2 k 2 + 2 k 3 + k 4 ) 6
k1 = f ( y n , t n )
k 2 = f ( yn + h h k1 , t n + ) 2 2
（3-11） ）
h h k3 = f ( yn + k2 , tn + ) 2 2 k 4 = f ( y n + hk 3 , t n + h)
第三章 数值积分法在系统仿真中的应用
亚当斯法(显式)
在解决积分问题时，采用亚当斯 贝喜霍斯显示多步法 简称亚当斯法。 贝喜霍斯显示多步法， 在解决积分问题时，采用亚当斯-贝喜霍斯显示多步法，简称 根据牛顿后插公式
Pi , j (t ) = fn +
(t − t n )(t − t n −1 ) ⋯ (t − t n +1− k ) k （3-25） (t − tn) ） ∇fn + ⋯ + ∇ fn k h h k!
第三章 数值积分法在系统仿真中的应用
第三章 数值积分法在系统仿真中的应 用
3.1 连续系统仿真中常用的数值积分法……………. 3.2 刚性系统的特点及算法…………………………. 3.3 实时仿真法………………………………………. ………………………………………. 3.4 分布参数系统的数字仿真………………………. 3.5 面向微分方程的仿真程序设计…………………. 本章小结……………………………………………….
第三章 数值积分法在系统仿真中的应用
为了适用于实时仿真计算，一般经常采用以下方法： 为了适用于实时仿真计算，一般经常采用以下方法：
（1）选择Adams多步法。 （2）合理地选择龙格-库塔法计算公式中的系数，使之适用于 实时仿真。
在方程 (3 − 18)式中，令 W1 = 0可得 W2 = 1, c2 = 1 , 此时， 2 (3 − 18)式化为
仿照显式多法的推倒过程， 仿照显式多法的推倒过程，得亚当斯-摩尔顿隐式多步法 的计算公式
yn +1 = yn + h∑ β k ,i f n −i +1
i =0
k −1
（3-33） ）
常用的四阶亚当斯预测-校正法的计算公式为 的计算公式为 常用的
预测：
y n +1 = y n +
p
h ( 55 f n − 59 f n −1 ) 24
当 k = 3时，3 − 28 ) 式的系数 v i 为 (
1 1 s ( s + 1) ds = 5 v2 = ∫ 2! 0 12
3-29） （3-29）
第三章 数值积分法在系统仿真中的应用
故可得三阶亚当斯公式
y (t n +1 ) = y (t n ) + hf n + 1 h∇f n + 5 h∇ 2 f n 2 12
采入u (t n + h ) 2
计算yn +1 , 并输入
计算k 2
计算下 一个k1
tn
tn + h 2
t n +1
图3.6
t n +1 + h 2
t n +1 + h
RK-2 的计算流程
由于当 t n + h = t n +1时极端即才能真正采入 外部输入 u n +1 , 才具备 计算 k 2的条件。所以 y n +1要到 t n +1 + h 时才能计算出来，并输 入到外 2 部设备，也就是说，计 算机输出要迟后半个计 算步距。