数值积分方法详解

连续系统仿真：数值积分法

2012-11-26

控制系统建模与仿真rjliu@

2

第二章连续系统仿真：数值积分方法

⏹

2.1 数值积分的基本概念⏹

2.2 龙格-库塔方法

⏹2.2.1 二阶RK 方法⏹2.2.2 四阶RK 方法⏹

2.2.3 多步法

⏹

2.3 数值计算稳定性分析

2012-11-26

控制系统建模与仿真rjliu@

3

问题描述

()(),(),()()

t t t t t f x A u B u x =+()()()()()()

t t t t t t =+⎧⎨

=+⎩x

Ax Bu y Cx Du 000,(),t t =x x 求:

x (t )= ?，t ≥t 0

(2-1)(2-2)对比上面两个问题,有:

y (t )是x (t ) u (t )的线性组合。

解出x (t )，就可计算出y (t )。x (t )由状态方程决定。数学上

的微分方程的初值问题。

已知: (1)状态方程(2)初值(3)输入u (t ) ，t ≥t 0

求: 系统的输出y (t )= ?, t ≥t 0已知: (1)微分方程

(2)初值()()(),d

x t f x t t dt

=000,(),t x t x =2012-11-26

控制系统建模与仿真rjliu@ 42.1 数值积分法的基本原理

微分方程的求解：

(1) 解析方法:只能求解一些特殊类型的方程.《高等数学》

(2) 数值解法:大部分的实际问题采用. 《计算方法》,《数值分析》差分方法是一类重要的数值解法.

⏹

1 数值积分法的原理

0()(),(0)x

t ax t x x == 0()at x t x e =解析解:

描点

t = [0, 1, 2, 3, 4, 5 ]x =[1, e -1,e -2,e -3,e -4,e -5]

1

2

3

4

5

00.10.20.30.40.50.60.70.80.91t

x (t )

e -t

为分析理解，作出函数的图形.

取a = -1, x 0=1, 则()t

x t e -=2012-11-26控制系统建模与仿真rjliu@ 52.1 数值积分法的基本原理

取点更密集一些，例如每隔0.25取一个点，

0.5

1

1.5

2

2.53

3.5

4

4.5

5

00.10.20.30.40.50.60.70.80.91t

x (t )

e

-t

00.51 1.52

2.53

3.54

4.55

0.1

0.20.30.40.50.60.70.80.91t

x (t )

e -t

t =[0, 0.25, 0.5, …]

2012-11-26控制系统建模与仿真rjliu@ 6

2.1 数值积分法的基本原理

⏹

1 数值积分法的原理

⏹

数值解法，就是寻求x=x (t )在一系列离散点t 1,t 2, …t k 上的近似解x 1, x 2, …,x k (此即数值解) .离散点满足：t 1

⏹

根据已知的初始条件x 0，采用“步进式”逐步递推的方法，依次计算出各时刻的数值x i .⏹

差分格式：

⏹由…x k -1, x k 计算x k +1的公式. x k +1＝？⏹

在仿真课程中，称为仿真模型

2012-11-26

控制系统建模与仿真rjliu@

7

2.1 数值积分的基本原理

001001121122

111

(1) (

)(2)()(3)()()()k k k k k x t x t t h x t x t t h x t x t t h x t x +++==+≈=+≈=+≈

已知初值 计算时刻的数值解 计算时刻的数值解 计算时刻的数值解 注1：步长:注2：一般计算中，h 0=h 1=h 2=…=h k =h ,等步长(固定步长)算法.注3：x (t k ) 表理论解,是精确的．x k 表示近似的数值解，含有误差

（舍入误差和截断误差）．

1k k k

h t t +=-数值积分求解的步骤：2012-11-26

控制系统建模与仿真rjliu@ 8

2.1 数值积分的基本原理

2 欧拉(Euler)法

微分方程(式2-2)的计算难点,在于其含有微分项. 因此要设法消

除微分项,这就是数值计算中的离散化. 方法:用差分代替微分.

取t = t k ，也即在t k 时刻, 微分方程变为：()()(),k k k d

x f x t d t t

t =2.1 显式格式: 一阶向前差分代替导数项1

()()()k k k x t x t d

x t dt h +-≈()

1()()(),k k k k x t x t h f x t t +≈+⋅()1,k k k k x x h f x t +=+⋅用x k 表示计算得到的近似值状态方程的仿真模型：

[]

1k k k k x x h Ax Bu +=+⋅+(2-3)

上式变成代数方程，2012-11-26控制系统建模与仿真rjliu@

10

gk=1;gT=0.5;tk = 0 ;% 系统时钟，初始时刻t 0= 0xk = 0;% 状态x k ≈x (t k ), 初始状态x 0= 0u0=100; % 输入信号u= 100 * 1(t)h = 0.05; % 仿真步长

timeend = 6*gT; % 仿真时间长度a = 1 -h/gT; b = gk*h/gT;

% 为了plot 作图, 状态x ,时种t 分别保存到数组t_eu1, x_eu1t_eu1(1)=tk; x_eu1(1)=xk;

for( ih = 2 : timeend/h +1 ) % ih 是循环变量, 也表示数组的下标

tk = tk + h;uk = u0;

xk = xk * a + b * uk;x_eu1(ih) = xk;t_eu1(ih) = tk;

end;

figure()；plot(t_eu1, x_eu1, ‘b’); % 蓝色线，数值解例2.1 采用欧拉方法的MATLAB 仿真程序2012-11-26

控制系统建模与仿真rjliu@

12

分析：

1) 红线表示理论解析解，蓝线表示欧拉法的仿真结果。两线没有重合，说明计算存在误差．误差范围[5%-0.06%]

例2-1 仿真结果分析

2) 如何减小误差？h 越小，差分越趋向于微分；减小h , 可以减小误差．