连续系统仿真：数值积分法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 =+⎧⎨=+⎩xAx 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)初值()()(),dx t f x t t dt=000,(),t x t x =2012-11-26控制系统建模与仿真rjliu@ 42.1 数值积分法的基本原理微分方程的求解：(1) 解析方法:只能求解一些特殊类型的方程.《高等数学》(2) 数值解法:大部分的实际问题采用. 《计算方法》,《数值分析》差分方法是一类重要的数值解法.⏹1 数值积分法的原理0()(),(0)xt 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]1234500.10.20.30.40.50.60.70.80.91tx (t )e -t为分析理解，作出函数的图形.取a = -1, x 0=1, 则()tx t e -=2012-11-26控制系统建模与仿真rjliu@ 52.1 数值积分法的基本原理取点更密集一些，例如每隔0.25取一个点，0.511.522.533.544.5500.10.20.30.40.50.60.70.80.91tx (t )e-t00.51 1.522.533.544.550.10.20.30.40.50.60.70.80.91tx (t )e -tt =[0, 0.25, 0.5, …]2012-11-26控制系统建模与仿真rjliu@ 62.1 数值积分法的基本原理⏹1 数值积分法的原理⏹数值解法，就是寻求x=x (t )在一系列离散点t 1,t 2, …t k 上的近似解x 1, x 2, …,x k (此即数值解) .离散点满足：t 1 <t 2 < t 3 < …< t k < …⏹根据已知的初始条件x 0，采用“步进式”逐步递推的方法，依次计算出各时刻的数值x i .⏹差分格式：⏹由…x k -1, x k 计算x k +1的公式. x k +1＝？⏹在仿真课程中，称为仿真模型2012-11-26控制系统建模与仿真rjliu@72.1 数值积分的基本原理001001121122111(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 kh t t +=-数值积分求解的步骤：2012-11-26控制系统建模与仿真rjliu@ 82.1 数值积分的基本原理2 欧拉(Euler)法微分方程(式2-2)的计算难点,在于其含有微分项. 因此要设法消除微分项,这就是数值计算中的离散化. 方法:用差分代替微分.取t = t k ，也即在t k 时刻, 微分方程变为：()()(),k k k dx f x t d t tt =2.1 显式格式: 一阶向前差分代替导数项1()()()k k k x t x t dx 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@10gk=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 , 可以减小误差．2012-11-26控制系统建模与仿真rjliu@13分析：1) h=0.02, 红线与蓝色线重合度高。

表明仿真解的误差减小．误差范围[2%-0.005%]例2-1 仿真结果分析2)还有没有其他减小误差的方法呢？研究其他数值积分方法2012-11-26控制系统建模与仿真rjliu@142.1 数值积分的基本概念欧拉公式的几何理解(之1) :从导数的概念理解：欧拉公式是一种折线近似. 看下图：t x0x 1x 2x 1t 2t kt 1k t +3t 0A1A 2A 3A kx 1n x +1k A +()()x x t =精确值未知的k x x =数值计算()010000011011100000110001111()(),(,).-(,)..[,)](, A A x t xt f t x A A t t h A x x x hf A A x t t t A t x t A t x t x ==+ 做直线,斜率为在的一阶导数即，直线与横坐标交于点，也即：欧拉公式是以直线来近似纵坐标就是纵轴方向的增量:＝对比发现上代替精确的在已知点()，求式就是欧时刻的点区间的式：拉公曲线段。

11211122122(,).A A A f t x A A t t h A x =+再在点做线段，斜率为在时刻时的位置，求得。

12n A A A 依此类推得到，，，表示2012-11-26控制系统建模与仿真rjliu@17局部截断误差:在x k 假定为准确的前提下(即x k =x (t k )),计算得到的x k +1(记作), 与精确值x (t k +1)的差.(只考虑第k 步的误差)2.4 欧拉格式的精度总体截断误差:它描述了从x 0, x 1, x 2,…计算到x k ，k 步误差的积累。

1ˆk x+(),0,1,2,k k ke x t x k =-= 111ˆ()k k k T x t x +++=-如果其局部截断误差为O (h r +1)，其总体截断误差为O (h r )。

收敛性: 对于如果表明当h 足够小时，x k 与x (t k )充分接近，数值解x k 收敛于精确值x (t k )则称该计算方法是收敛的。

(),1r k e O h r =≥0,0k h e →→当有r 阶方法：一个方法，它的总体截断误差若为O (h r )，则它的精度是r 阶的，称其为r 阶方法。

2012-11-26控制系统建模与仿真rjliu@182.4 欧拉格式的精度1. 欧拉公式(2-3):()()1ˆ,()(),()()k k k k k k k k k x x h f x t x t h f x t t x t h x t +'=+⋅=+⋅=+⋅21()()()()()2k k k k h x t x t h x t h x t x ξ+'''=+=++⋅⋅根据泰勒公式,精确值x (t k +1)精确值为:()0221122ˆ()():max ()2nk k t t h x t xx h M M x ξξξ++≤≤''''-=≤=⋅⋅注单步显式的欧拉公式的精度是一阶的.因为前向差分与后向差分的精度相同,所以隐式的欧拉公式也是一阶精度的.局部截断误差()1,k k k k x x h f x t +=+⋅2012-11-26控制系统建模与仿真rjliu@202.1 数值积分的基本概念000011000102211121100001()()()()(,),,()()(,)(),()()(,)(),()()(,(,)(,))(()tk t k k k k x t x t x t x t f x t t t x f x t dtt t t t x t x t f x t t t t t x t x t f x t t t t t x t x t f x t t f x t t t ++=+≈+=≈+-=≈+-=≈+--⎰当的时候，在几何意义上，上式右侧的红色部分，是由，,轴所围成的面积。

有有 很小 有 1)k k t +-欧拉公式采用矩形面积表示定积分。

如果采用梯形面积，则更准确一些。

欧拉公式的几何理解(之2) 数值积分求解角度，如何计算x (t )?考虑高等数学中的定积分公式：2012-11-26控制系统建模与仿真rjliu@212.1 数值积分的基本概念⏹3 梯形公式: 欧拉法用矩形的面积表示积分项。

更精确一些，可以用梯形的面积来表示积分项。

梯形差分格式是:[]1111(,)2(,)k k k k k k x x h f x t f x t +++=++()()1111()()(),2(),k k k k k k f x t x t h x t x t f t t +++⎡⎤≈++⎣⎦(2-6)状态空间模型:()()1111()()2k k k k k k Ax Bu x x h Ax Bu t t +++⎡⎤=++++⎣⎦现在是t k 时刻, x k , f (x k , t k )已知，求x k +1.而x k +1出现在等式的两端，故梯形格式是隐式的,计算比较复杂.2012-11-26控制系统建模与仿真rjliu@222.1 数值积分的基本概念⏹3 改进的欧拉公式: 先用欧拉公式预报出x k +1的近似值,再把这个近似值代入梯形公式,再利用梯形公式计算x k +1的校正值.()()1111(,),,2k k k k k k k k k k x h f x t hx x f x t f t xx ++++=+⋅⎧⎪⎨⎡⎤=++⎪⎣⎦⎩ 预报公式校正公式()()1211111211(,) ,2k k k k k k k k k k K f x t x h K f t t h hx K x x K K x t +++++=⎧⎪=+⎧⎨=⎨⎪=+⎩⎩=++ 其中因为f (·)是x 对t 的一阶导数,它的几何意义表示斜率.因此，有(2-7)(2-8)总体截断误差O (h 2)，精度是二阶的．在实时仿真的应用中较多.3.1 表示形式1：本式与RK4表示形式，是统一的。