2013 第3章 数值积分算法

u0
1.计算所有变量的第一个RK系数向量
1 1 0 x10 0 k11 k1 A k21 x20 2 0.5 0 2
2. 计算所有变量的第二个RK系数向量 k 2
x10 h k11 0 1 1 0.1 0 0.1 k12 k2 A x k k 2 0.5 0 2 1.95 2 2 22 21 20
对于线性定常系统:
(t ) Ax (t ) Bu(t ) x
f (x, t ) Ax Bu(t )
RK4的4个 k 向量可表示为
k1 f (tn , x n ) Ax n Bu(tn ) k f (t h , x h k ) A x h k Bu (t h ) 2 n n 1 1 n n 2 2 2 2 k f (t h , x h k ) A x h k Bu (t h ) n n 2 2 n n 3 2 2 2 2 k 4 f (tn h, x n hk 3 ) A x n hk 3 Bu (tn h)
t10 1,

x10 x9 0.1x92 0.4682

3.1 欧拉法和改进的欧拉法
(例1 续) 后向欧拉法的递推公式为
xn1 xn hf (tn1, xn1 ) xn 0.1x
解此方程可得
2 n1
xn1 10 xn 25 5
3.1 欧拉法和改进的欧拉法
f (t, x) x2
h f (tn , xn ) f (tn1 , xn1 ) 2 2 2 xn 0.05( xn xn 1 )
2 xn 1 20 xn xn 100 10
解此代数方程可得：
t0 0, x0 1
2 t1 0.1, x1 20 x0 x0 100 10 20 1 12 100 10 0.9087
显式算法：计算 xn+1 时，没有用到 tn 时刻以后的状态或输入。
隐式算法：计算xn+1 时，用到 tn 时刻以后的状态或输入。
单步法：计算 xn+1 时只与 tn 时刻的状态或输入有关。 多步法：计算 xn+1 时要用到 tn 时刻以前的状态或输入。
3.1 欧拉法和改进的欧拉法
一、前向欧拉(Euler)算法计算公式：
梯形法的几何意义
3.1 欧拉法和改进的欧拉法
2 x(t ) x (t ) 0 ，x(0) 1, x 0 0 t 1 试用梯形法求其数值解（取仿真步距 h 0.1)，
例2 设系统方程为
(t ) x2 (t ), 解：原方程为 x
梯形法的递推公式为
xn 1 xn
3.1.3 线性多步法
一、显式多步积分算法
四阶亚当姆斯积分算法
h xn 1 xn [9 f (t n 3 , xn 3 ) 37 f (tn 2 , xn 2 ) 59 f (t n 1 , xn 1 ) 55 f (t n , xn )] 24
3.1.2 龙格-库塔法
t0 0, x0 1
2 k1 f (0, x0 ) x0 1
h h h 2 k 2 f ( , x0 k1 ) ( x0 ) 0.9025 2 2 2 h h h 2 k3 f ( , x0 k 2 ) ( x0 0.9025 ) 0.9118 2 2 2
二阶、单步、显式
3.1.2 龙格-库塔法
一、龙格-库塔（Runge-Kutta）积分算法思路 间接利用泰勒展开式。用在若干个点上函数值
f(t,y) 的线性组合来代替高阶导数项，既可以避免计
算高阶导数，又可以提高数值计算精度。
3.1.2 龙格-库塔法
二、二阶Rungeቤተ መጻሕፍቲ ባይዱKutta法
xn 1 xn hk2 k1 f (tn , xn ) h h k2 f (tn , xn k1 ) 2 2
xn1 xn hf (tn , xn ), n 0, 1, , M 1
xn 1
截断误差：
R(tn1 ) O(h2 )
x(t )
xn
x(t )
一阶、单步、显式
tn
h
tn 1
t
3.1 欧拉法和改进的欧拉法
二、后向欧拉算法计算公式：
xn1 xn hf (tn1, xn1 ), n 0, 1, , M 1
k4 f (h, x0 hk3 ) ( x0 0.9118hk3 )2 0.8260
h t1 0.1 , x1 x0 (k1 2k 2 2k3 k 4 ) 0.9091 6
3.1.2 龙格-库塔法
k1 f (h, x1 ) x12 0.8264
3h h h 2 k 2 f ( , x1 k1 ) ( x1 0.8264 ) 0.7530 2 2 2 3h h h 2 k3 f ( , x1 k 2 ) ( x1 0.7530 ) 0.7594 2 2 2
k4 f (2h, x1 hk3 ) ( x1 0.7594h)2 0.6942
三、梯形公式（改进的欧拉公式）
h xn 1 xn [ f (t n , xn ) f (t n 1 , xn 1 )], 2 n 0, 1, , M 1
A
x(t )
xn 1 xn
x(t )
截断误差：
R(tn1 ) O(h3 )
t
二阶、单步、隐式
tn
h 2
h 2
tn 1
4. 计算所有变量的第四个RK系数向量
x10 1 1 k13 0 k14 0.0975 0.196125 k4 A h 0 0.1 1.96125 1.921438 x k k 2 0.5 24 23 20
二阶、单步、显式
一种改进的Euler 公式
3.1.2 龙格-库塔法
三、经典显式四阶Runge-Kutta法 (RK4)
h xn 1 xn 6 (k1 2k2 2k3 k4 ) k f (t , x ) n n 1 h h k2 f (tn , xn k1 ) 2 2 h h k3 f (tn 2 , xn 2 k2 ) k f (t h, x hk ) 4 n n 3
第3章 连续时间系统 数字仿真算法
章节概览
3.1 欧拉法和改进的欧拉法
3.2 龙格－库塔法 3.3 线性多步算法
3.4 数值积分算法的选用
数值积分算法
(t ) f (t , x(t ), u (t )), t t0 x 系统状态方程的一般形式： x(t0 ) x0
数值积分算法的类型
2 x(t ) x (t ) ， f (t, x) x2
0 t 1
前向欧拉法的递推公式为:
xn1 xn hf (tn , xn ) xn 0.1x n2
t0 0, x0 1 2 t1 0.1, x1 x0 0.1x0 0.9
t2 0.2, x2 x1 0.1x12 0.819
y(t )
yn yn1
截断误差：
R(tn1 ) O(h2 )
y(t )
一阶、单步、隐式
tn
h
tn 1
t
3.1 欧拉法和改进的欧拉法
例1 设系统方程为
(t ) x (t ) 0 , x(0) 1, x 0 x
2
试用Euler法求其数值解（取仿真步距 h 0.1 )，
解：原方程为
h xn 1 xn (k1 2k 2 2k 3 k 4 ) 6
3.1.2 龙格-库塔法
例4 已知系统方程为
(t ) 0.5 y (t ) 2 y(t ) 0 , y
t 0.1 , 0.2 时
(0) 0 , y
y(0) 1
取积分步长 h 0.1，试用四阶龙格库塔法计算
t2 0.2, x2 20 x1 x12 100 10 20 0.9087 0.9087 2 100 10 0.8328
3.1.1 欧拉法和改进的欧拉法
四、预测-校正法（改进的欧拉公式）
h x x ( k k ) n 1 n 1 2 2 k1 f (tn , xn ) k f (t h, x hk ) n n 1 2
h t2 0.2 , x2 x1 (k1 2k2 2k3 k4 ) 0.8333 6
h t10 1 , x10 x9 (k1 2k 2 2k3 k 4 ) 0.5000 6
3.1.2 龙格-库塔法
四、矩阵分析法（RK4解状态方程）
5. 计算第一步的近似值。
x11 x10 h x x 6 k1 2k 2 2k 3 k 4 21 20 1 0.1 0 0.1 0.0975 0.196125 1.0099 2 2 0 6 2 1.95 1.96125 1.921438 0.1957
3.1.2 龙格-库塔法
3. 计算所有变量的第三个RK系数向量
k3
k4
x h k 0 1 1 0.1 0.1 0.0975 k k 3 13 A 10 12 k23 x20 2 k22 2 0.5 0 2 1.95 1.96125
y 的值。
解： 原系统的状态方程为
x1 y, x2 y
1 0 1 x1 x x x 2 0.5 2 2
x1 (0) 1 x2 (0) 0
3.1.2 龙格-库塔法
1 0 , B 0, 对应于系数矩阵为 A 2 0.5
四阶、单步、显式
3.1.2 龙格-库塔法
例3 设系统方程为
(t ) x2 (t ) 0, x(0) 1 x
试用四阶Runge—Kutta法求其数值解
（取仿真步距h 0.1,0 t 1)
(t ) x (t ), 解：原方程为 x
2
f (t, x) x2
h 递推公式可写为 xn 1 xn (k1 2k2 2k3 k4 ) k1 k2 k3 k4 6 2 f (tn , xn ) xn h h h 2 f (tn , xn k1 ) ( xn k1 ) 2 2 2 h h h 2 f (tn , xn k2 ) ( xn k2 ) 2 2 2 f (tn h, xn hk3 ) ( xn hk3 ) 2