6.数值积分法仿真

试用梯形法求其数值解，取步长h=0.1, 0 ≤ t ≤ 1
! ( t ) = − y ( t ), 即f ( t , y ) = − y ( t ) 解: y
梯形法递推式： h y n = y n −1 + [f (t n −1, y n −1) + f (t n , y n )] 2 2 2 = y n −1 − 0.05( y n −1 + y n ) 隐式算法，需先解此非线性方程：
r
kw = ω1k1 + ω2k 2 + ! + ω r k r =
∑
i =1
ωi ki
后面Runge—Kutta法就是采用这样的思想来 进行计算的。
26
梯形法大大提高了精度，但为隐式算法，每次 递推计算时需解一次非线性方程，计算量较大 由此考虑进行改进，先用Euler法计算出： x (t n +1) 的近似值 x pn +1 代入导函数 fn +1 求出近似值 f
第6 章
数值积分法仿真
1
第6 章
数值积分法仿真
本章主要教学内容
本章主要介绍控制系统数学模型的相关知识， 通过本章的学习，应掌握以下内容：
Ø 求解常微分方方程数值解的一一般方方法 Ø 数值积分法的基本概念及其常用用方方法 Ø 以系统微分方方程或传递函数作为数学模型的仿真过程
及程序设计方方法 Ø 以系统动态结构图作为数学模型的仿真过程及程序设计方方法 Ø 仿真步⻓长的选择与系统仿真精度和稳定性的对应关系 Ø 快速仿真算法的概念、特点及其应用用
tn+1
11
t
３、切线近似（１）
t n 的一个小邻域内 [t n , t n +1] ，曲线x(t)可用 在 t n 处的切线来表示， x(t) 在 t n 处的斜率为：
dx = f (t n , x n ,un (t )) dt t = t n
x = x n + f (t n , x n ,un )(t − t n ) x (t n +1) t n +1 取切线上 处的值来近似 x (t n +1) ≈ x (t n ) + hf (t n , x n ,un )
前向Euler法递推式： y n = y n −1 + hf (t n −1, y n −1) 2 = y n −1 + 0.1× ( − y n −1)
= y n −1(1 − 0.1y n −1) 有初始条件： t 0 = 0, y 0 = 1
可进行递推：
t1 = 0.1, y1 = y 0 (1 − 0.1y 0 ) = 1 − 0.1 = 0.9
p n +1
再代入梯形公式求解。
= xn + hf (t n , xn ) 预估公式（Euler法） h c p x n+1 = xn + f (t n , xn ) + f (t n+1 , x n+1 ) 2
n +1
x
p
(
)
校正公式（梯形法）
27
为预估—校正法，也称为改进的Euler法。
! (t ) + y 2 (t ) = 0, y (0) = 1 例：设系统方方程 y
23
f(t) 误差 fn
fn+1
f(t) fn
误差
fn+1
左矩形近似
右矩形近似
tn
tn+1
t
tn
tn+1
t
f(t) 误差 fn
fn+1
梯形近似
24
tn
tn+1
t
梯形法实质是采用了 (t n , x (t n )) 和 (t n +1, x (t n +1)) 两点斜率平均值的结果，由于利用了两点的信息， 1 ⎧ 从而提高了计算精度。 x = x + h(k + k )
20
f(t) 误差 fn
fn+1
f(t)
左矩形近似
tn
tn+1
t
t
f(t)
f(t)
t
21
t
Ø 梯形法
Euler法的计算精度较差，如果改用梯形面积 代替每个步距的曲线面积，就可提高精度。 精确积分应为曲边梯形的面积： 现用直边梯形的面积来近似：
t n +1
t n +1
tn
∫
f (t , x ,u (t ))dt
3
6.1.2 离散化原理
• 在数字计算机上对连续系统进行仿真时，首先遇到的问题是，数字 计算机的数值及时间都是离散的（计算精度，指令执行时间），而 被仿真系统的数值和时间是连续的，后者如何用前者来实现？ ! = f ( y, u, t ) ,其中u(t)为输入变量，y(t)为系统状态变 y • 设系统模型为： 量。令仿真时间间隔为h,离散化后的输入变量为u’(tk), 其中tk表示 t=kh。如果 u’(tk) ≈ u(tk), y’(tk) ≈ y(tk)，则认为两模型等价，称为 相似原理。
2
6.1
6.1.1 概述
数值积分法
数字仿真模型、算法及仿真工工具
控制系统的数字仿真是利用数字计算机作为仿真工具，采用数学上的各 种数值算法求解控制系统运动的微分方程，得到被控物理量的运动规律。 通常，计算机模拟被控对象是用一定的仿真算法来实现被控对象的运动规 律，这是基于被控对象的数学模型来完成的。 控制系统的数学模型经过合理的近似及简化，大多数建立为常微分方程 的表达形式。由于数学计算的难度和实际系统的复杂程度，在实际中遇到 的大部分微分方程难以得到其解析解，通常都是通过数字计算机采用数值 计算的方法来求取其数值解。在高级仿真软件（例如 MATLAB ）环境下， 已提供了功能十分强大、且能保证相应精度的数值求解的功能函数或程序 段，使用者仅需要按规定的语言规格调用即可，而无需从数值算法的底层 考虑其编程实现过程。
用右矩形面积近似该区间的曲线面积
x (t n +1) ≈ x (t n ) + f (t n +1, x n +1,un +1)h
得到：
x n +1 = x n + hfn +1
这是右矩形欧拉公式，是一个隐式算法．
10
右矩形（也称为后向欧拉法）近似及误差
f(t) fn
误差
fn+1
右矩形近似
tn
""(t n ) 2 x " (t n )h + x (t n +1) = x (t n + h) = x (t n ) + x h +! 2!
= x (t n ) + f (t n , x n ,un )h
写成差分方程为 x n +1 = x n + hf (t n , x n ,un ) 这就是解微分方程初值问题的欧拉算法。
18
前向欧拉yn = yn−1 + hf (tn−1 , yn−1 )
后向欧拉yn = yn−1 + hf (tn , yn )
f(t) 误差 fn fn+1
f(t) fn
误差
fn+1
左矩形近似
右矩形近似
tn
tn+1
t
tn
tn+1
t
19
前向Euler法在不同步长的结果比较
y n = y n −1 + hf (t n −1, y n −1)
y(t ) =
y + ∫t
0
t
0
f (t , y)dt
右图所示曲线下的面积就 是y(t),由于难以得到积分 的数值表达式，所以采用 近似的方法，常用有三种形式： Ø 欧拉法 Ø 梯形法 Ø 龙格一库塔法
5
Ø 欧拉法（Euler Method)
Euler法是最简单的一种数值积分法的单步 运算，虽然计算精度较差，但几何意义明显， 便于理解，能说明构造数值积分算法的基本思 想。 下面采用三种方法推导出Euler法的数值近 似公式，以便对数值积分器的基本思想能透彻 了解。
6
１、Taylor级数展开
! (t ) = f (t , x ,u(t )) 以一阶连续系统为例， ⎧ x ⎨ x (t ) = x 微分方程及初值如右： 0 ⎩ 0 x(t)为解析解，将x(t)展开成Taylor级数
只取一次项，其余忽略 x (t ! ) ≈ x ( t ) + x (t n )h n +1 n
(t n , x n ) 以 fn 为斜率的切线方程为： 过点
也能得到： x n +1 = x n + hfn
前向欧拉法
12
３、切线近似（２）
[t n , 曲线 t n +1] t n +1 处的切线来表示， 在 x(t) 可用 x(t) 在 t n +1 处的斜率为：
dx = f (t n +1, x n +1,un +1(t )) dt t = t n +1
15
t 2 = 0.2, y 2 = y1(1 − 0.1y1) = 0.9 × 0.91 = 0.819 ! t10 = 1, y10 = y 9 (1 − 0.1y 9 ) = 0.4682 !
后向Euler法递推式：
y n = y n −1 + hf (t n , y n )
= y n −1 − 0.1y n 2
过点 (t n , x n )
以 fn+1 为斜率的切线方程为：
x = x n + f (t n +1, x n +1,un +1)(t − t n ) x (t n +1) t n +1 取切线上 处的值来近似 x (t n +1) ≈ x (t n ) + hf (t n +1, x n +1,un +1)