计算仿真课件 第三章1

济南大学控制科学与工程学院
计算机仿真技术
3
3.1.1 欧拉法和改进欧拉法
欧拉法是最简单的单步法，它是一阶的，精度较差. 但公式简单，有明显的几何意义， 适合初学者在直观上学习数值x(tn)是如何逼近微分方程的精确解x(t)的.
考虑初值问题：
x
dx dt
f t , x ,
x t 0 x0
计算机仿真技术
济南大学控制科学与工程学院 授课教师：李实
cse_lis@
1教1007室
第三章 数字仿真通用算法
在工程实际和科学研究中所遇到的问题往往很复杂， 很多情况下无法给出描述动态特性的微分方程解的解析 表达式，多数只能用近似的数值方法求解. 随着计算机软 硬件和数值理论的进展，微分方程的数值解方法已成为 当今研究、分析和设计系统的一种有力工具.
取
1 2!
h 2 t n x
(3.8)
Rn
1 2!
h 2 t n x
并截断，得到欧拉公式，Rn为截断误差，与h2成正比
Rn O h 2

(3.9)
解从t=0到t=tn所积累的误差为整体误差，比局部误差要大，欧拉法的整体误差与h成正比， 即为O1(h).
(3.1)
式 （ 3.1 ） 的 解 x(t) 是 一 连 续 变 量 x 的 函 数 ， 现 在 要 用 一 系 列 离 散 时 刻 的 近 似 值 x(t1),x(t2),…,x(tn)来代替，其中ti=t0+ih, h成为步长，是相邻两点的距离. 将式（3.1）在区 间（ti,ti+1）上积分，则可得：
(3.13)
因为λ为负，所以（3.13）必然成立，说明后向欧拉法是恒稳定的. 改进欧拉法（梯形公式）的稳定性： 改进欧拉法公式：
h h xn 1 xn xn xn 1 1 xn 1 xn 2 2 2 1 h / 2 xn 1 xn 1 h / 2
t 0 0,
y0 1
2 2
t1 0.1, t 2 0 . 2, t10 1.0,
济南大学控制科学与工程学院
y1 y 0 0.1 y 0 0.9 y 2 y1 0.1 y1 0.819 y10 y9 0.1 y9 0.4682
x t n h x t n hx t n
1 2!
h 2 t n x
1 3!
h 3t n x
(3.15)
欧拉法是截去h2以后各项所得到的一阶一步法，因此精度较低. 如果将泰勒展开式多选几 项以后截断，可以得到精度较高的高阶数值解，但是直接使用泰勒展开要计算函数的高阶 导数. 龙格-库塔法采用间接利用泰勒展开式的思路，用在n个点上的函数值f的线性组合来代替f的 导数，然后按泰勒展开式确定其中的系数，以提高算法的精度. 这样既避免计算函数的导数， 同时又保证了计算精度. 由于龙格-库塔法具有许多优点，因此在许多仿真程序包中，它是一个最基本的算法之一.
xn 1 xn xn 1h xn 1 f xn , t n
济南大学控制科学与工程学院
计算机仿真技术
8
改进欧拉法（预测-校正法）
xti 1 xti
x ti 1 x ti
写成递推差分格式：
t i 1 ti
f (t , x ) dt
济南大学控制科学与工程学院
计算机仿真技术
10
欧拉法：误差分析
舍入误差：由于计算机进行计算时，数的位数有限引起的，一般舍入误差与h-1成正比，可 以得到欧拉法总误差为：
n O1 h O2 h 1
(3.10)
由式（3.10）可以看出，步长h增加，截断误差增加，反之，步长h减小，舍入误差加大.
h
该差分方程稳定的条件：
1 h / 2 2 h / 2 2 1 h / 2 2 h / 2 2
计算机仿真技术
1
(3.14)
同理，改进欧拉法是恒稳定的.
济南大学控制科学与工程学院 13
3.1.2 龙格-库塔法
将x(tn+h)在t=tn点进行泰勒展开：
计算机仿真技术 4
欧拉法
因此，式（3.2）可近似为：
x ti 1 x ti hf (ti , x (ti ))
写成递推式：
x t n 1 x t n hf (t n , x (t n )),
n 0,1,2,..., N
(3.3)
已知x(0)=x0, 可以由上式求出x(t1), x(t2), …. 这种算法成为单步法. 可以直接由初始值递推 出其它各时间的值，因此单步法是一种自启动算法. 定义x(tn)=xn, 式子（3.3）还可以写成：
• 3.1 系统仿真中常用的数值积分法 • 3.2 刚性系统的特点及算法 • 3.3 实时半实物仿真 • 3.4 采样控制系统的仿真方法 • 3.5 分布参数系统的数字仿真
济南大学控制科学与工程学院 计算机仿真技术 2
3.1 系统仿真中常用的数值积分 法
• • • • 3.1.1 欧拉法和改进欧拉法 3.1.2 龙格-库塔法 3.1.3 线性多步法 3.1.4 MATLAB中常微分方程求解方法
济南大学控制科学与工程学院
y0 1
y1 0 .9161 y10 0 .5165
7
计算机仿真技术
欧拉法的几何意义
前向欧拉法 x(t) x(t) 后向欧拉法
xn+1
k xn
xn h tn tn+1 t
xn+1 xn
k xn 1
h tn tn+1 t
x n 1 x n x n h
xn 1 xn xn h,
该式称为前向欧拉公式，又称显式欧拉公式. 后向欧拉公式（隐式欧拉公式）：
n 0,1,2,..., N
(3.4)
xn 1 xn xn 1h,
济南大学控制科学与工程学院 计算机仿真技术
n 0,1,2,..., N
(3.5)
5
欧拉法
[例3.1] 设系统方程为
龙格-库塔法
当r=2时，二阶泰勒展开式为： 由式（3.17）得：
xn 1 xn hf
1 2
h 2 ( ft f x f )
x n 1 x n W1 K 1 W 2 K 2 K 1 hf t n , x n K hf t c h , x a K n 2 n 21 1 2
x
f t , x dt

f t

f dx x dt
ft f x f
f f d d d ( f t f x f ) df t df x df df t x x fx f fx f dt dt dt dt dt dt dt 2 f 2 f 2 f f f 2 f tt tx x f x ( t x x ) f xt xx f tt f tx f t f x ( f t f x f ) f ( f xt f xx f )
济南大学控制科学与工程学院 f tt f xx
x
15
f 2 f x2 f f x f t 2 f xt f 计算机仿真技术
龙格-库塔法
将各阶导数代入方程（3.15），得到（3.16）：
xn 1 xn hf
1 2!
h ( ft f x f )
2
1 3!
h 3 ( f tt f xx f 2 f x2 f f x f t 2 f xt f )
(3.16)
Wi: 待定的权因子， r：解公式的阶数，
假设问题的数值解公式为：
r x n 1 x n Wi K i i 1 (3.17) i 1 K hf t c h, x a K n ij j i i n j 1
济南大学控制科学与工程学院
计算机仿真技术
14
龙格-库塔法
将x(tn+h)在t=tn点进行泰勒展开：
x t n h x t n hx t n
根据偏导数关系：
1 2!
h 2 t n x
1 3!
h 3t n x
(3.15)
x f t , x f
xn 1 xn xn h xn h xn
1 h 1
或：
(3.11) (3.12)
h
2

选择的步长过大会使算法变得不稳定
济南大学控制科学与工程学院 计算机仿真技术 12
欧拉法：稳定性分析
后向欧拉法的稳定性：
后向欧拉法公式：
该差分方程稳定的条件：
xn 1 xn xn 1h xn h xn 1 1 xn 1 xn 1 h 1 h 1
济南大学控制科学与工程学院
计算机仿真技术
11
欧拉法：稳定性分析
求解微分方程的另一个重要问题是数值解是否稳定. 假设方程的特征根为λ. 方程可以写成：
x (t ) x (t )
若该方程稳定，需要满足 前向欧拉法的稳定性： 前向欧拉法公式： 该差分方程稳定的条件：
(3.10)
j 0
(3.2)
对积分公式（3.2）利用梯形面积公式计算积分项，得 f ti 1 , x(ti 1 )
2 f n 1
(3.6)
xn 1 xn
h 2
xn xn 1 xn h f n
(3.7)
可以先用欧拉法计算xn+1, 定义为预估值，写为xpn+1，再将该预估值代入原方程中计算函数， 最后利用式（3.7）得到修正后的xcn+1, 改进后的欧拉法描述如下：