计算仿真课件第三章

1、计算机仿真技术 授课教师：李实 1教1007室 1计算仿真课件第三章 济南大学控制科学与工程学院计算机仿真技术2 第三章 数字仿真通用算法 在工程实际和科学研究中所遇到的问题往往很复杂， 很多情况下无法给出描述动态特性的微分方程解的解析 表达式，多数只能用近似的数值方法求解. 随着计算机软 硬件和数值理论的进展，微分方程的数值解方法已成为 当今研究、分析和设计系统的一种有力工具. 3.1 系统仿真中常用的数值积分法 3.2 刚性系统的特点及算法 3.3 实时半实物仿真 3.4 采样控制系统的仿真方法 3.5 分布参数系统的数字仿真 济南大学控制科学与工程学院计算机仿真技术3 3.1 系统仿真中

2、常用的数值积分 法 3.1.1 欧拉法和改进欧拉法 3.1.2 龙格-库塔法 3.1.3 线性多步法 3.1.4 MATLAB中常微分方程求解方法 济南大学控制科学与工程学院计算机仿真技术4 3.1.1 欧拉法和改进欧拉法 欧拉法是最简单的单步法，它是一阶的，精度较差. 但公式简单，有明显的几何意义， 适合初学者在直观上学习数值x(tn)是如何逼近微分方程的精确解x(t)的. 考虑初值问题： 00 ,xtxxtf dt dx x (3.1) 式 （ ） 的 解 x ( t ) 是 一 连 续 变 量 x 的 函 数 ， 现 在 要 用 一 系 列 离 散 时 刻 的 近 似 值 x(t1),x

3、(t2),x(tn)来代替，其中ti=t0+ih, h成为步长，是相邻两点的距离. 将式（）在区间 （ti,ti+1）上积分，则可得： 1 ),( 1 i i t t ii dtxtftxtx(3.2) 积分项的几何意义是曲线f(x,y)在区间（ti,ti+1）上的面积，当（ti,ti+1）充分小时，可以用矩 形面积来代替： )(,),( 1 ii t t txthfdtxtf i i 济南大学控制科学与工程学院计算机仿真技术5 欧拉法 因此，式（）可近似为： )(,( 1iiii txthftxtx 写成递推式： Nntxthftxtx nnnn ,.,2 , 1 , 0),(,( 1 (3

4、.3) 已知x(0)=x0, 可以由上式求出x(t1), x(t2), . 这种算法成为单步法. 可以直接由初始值递推 出其它各时间的值，因此单步法是一种自启动算法. 定义x(tn)=xn, 式子（）还可以写成： Nnhxxx nnn ,.,2 , 1 , 0, 1 (3.4) 该式称为前向欧拉公式，又称显式欧拉公式. 后向欧拉公式（隐式欧拉公式）： Nnhxxx nnn ,.,2 , 1 , 0, 11 (3.5) 济南大学控制科学与工程学院计算机仿真技术6 欧拉法 例3.1 设系统方程为 试用欧拉法求其在t=1时的数值解，仿真步长，要求分别使用前向欧拉法和后向欧拉法. 100)( 2 yt

5、yt y 方程可以写为： )(,)( 2 tyytfty 写出该系统向欧拉递推公式为： 2 1nnn hyyy 1, 0 00 yt 4682. 01 . 0, 0 . 1 819. 01 . 0, 2 . 0 9 . 01 . 0, 1 . 0 2 991010 2 1122 2 0011 yyyt yyyt yyyt 济南大学控制科学与工程学院计算机仿真技术7 欧拉法 该系统的后向欧拉递推公式为： 2 11 nnn hyyy 1, 0 00 yt 5165. 0, 0 . 1 9161. 0, 1 . 0 1010 11 yt yt 0 1 2 1 nnn yyhy 整理后得到： 将代入，

6、求得方程的解为： 52510 1 nn yy 济南大学控制科学与工程学院计算机仿真技术8 欧拉法的几何意义 n xk h t x(t) tn xn tn+1 xn+1 hxxx nnn 1 前向欧拉法前向欧拉法 tn+1 1 n xk h t x(t) tn xn xn+1 hxxx nnn11 后向欧拉法后向欧拉法 nnn txfx, 1 济南大学控制科学与工程学院计算机仿真技术9 改进欧拉法（预测-校正法） 对积分公式（）利用梯形面积公式计算积分项，得： 1 ),( 1 i i t t ii dtxtftxtx(3.2) )(,)(, 2 111 iiiiii txtftxtf h txt

7、x (3.6) 写成递推差分格式： 111 22 nnnnnnn ff h xxx h xx (3.7) 可以先用欧拉法计算xn+1, 定义为预估值，写为xpn+1，再将该预估值代入原方程中计算函数， 最后利用式（）得到修正后的xcn+1, 改进后的欧拉法描述如下： p n p nn p n fxtfx 1111 , , 2 , 1 2 11 1 n xx h xx xhxx p nnn c n nn p n (3.8) 济南大学控制科学与工程学院计算机仿真技术10 欧拉法：误差分析 欧拉法以一定的步长来进行数值计算，所得到的解在tn点的近似值与微分方程的精确解存 在误差. 数值仿真的误差可分

8、为截断误差和舍入误差两种. 截断误差与采用的计算方法有关，而舍入 误差由计算机的字长决定. 截断误差：将x(tn+h)在t=tn点进行泰勒展开： nnnn txhtxhtxhtx 2 ! 2 1 (3.8) 取并截断，得到欧拉公式，Rn为截断误差，与h2成正比 nn txhR 2 ! 2 1 2 hORn (3.9) 解从t=0到t=tn所积累的误差为整体误差，比局部误差要大，欧拉法的整体误差与h成正比， 即为O1(h). 济南大学控制科学与工程学院计算机仿真技术11 欧拉法：误差分析 舍入误差：由于计算机进行计算时，数的位数有限引起的，一般舍入误差与h-1成正比，可 以得到欧拉法总误差为：

9、1 21 hOhO n (3.10) 由式（）可以看出，步长h增加，截断误差增加，反之，步长h减小，舍入误差加大. 济南大学控制科学与工程学院计算机仿真技术12 欧拉法：稳定性分析 求解微分方程的另一个重要问题是数值解是否稳定. 假设方程的特征根为. 方程可以写成： )()(txtx(3.10) 若该方程稳定，需要满足 0j 前向欧拉法的稳定性： 前向欧拉法公式： nnnnn xhxhxxx 1 该差分方程稳定的条件： 11 h (3.11) 或： 2 h (3.12) 选择的步长过大会使算法变得不稳定 济南大学控制科学与工程学院计算机仿真技术13 欧拉法：稳定性分析 后向欧拉法的稳定性： 后

10、向欧拉法公式： 111 nnnnn xhxhxxx 该差分方程稳定的条件： 11 h (3.13) nn x h x 1 1 1 因为为负，所以（）必然成立，说明后向欧拉法是恒稳定的. 改进欧拉法（梯形公式）的稳定性： 改进欧拉法公式： 111 22 1 2 nnnnnn x h x h xx h xx nn x h h x 2/1 2/1 1 该差分方程稳定的条件： 1 2/2/1 2/2/1 22 22 hh hh (3.14) 同理，改进欧拉法是恒稳定的. 济南大学控制科学与工程学院计算机仿真技术14 3.1.2 龙格-库塔法 将x(tn+h)在t=tn点进行泰勒展开： nnnnn tx

11、htxhtxhtxhtx 32 ! 3 1 ! 2 1 (3.15) 欧拉法是截去h2以后各项所得到的一阶一步法，因此精度较低. 如果将泰勒展开式多选几 项以后截断，可以得到精度较高的高阶数值解，但是直接使用泰勒展开要计算函数的高阶 导数. 龙格-库塔法采用间接利用泰勒展开式的思路，用在n个点上的函数值f的线性组合来代替f的 导数，然后按泰勒展开式确定其中的系数，以提高算法的精度. 这样既避免计算函数的导数， 同时又保证了计算精度. 由于龙格-库塔法具有许多优点，因此在许多仿真程序包中，它是一个最基本的算法之一. 济南大学控制科学与工程学院计算机仿真技术15 龙格-库塔法 将x(tn+h)在t

12、=tn点进行泰勒展开： nnnnn txhtxhtxhtxhtx 32 ! 3 1 ! 2 1 (3.15) 根据偏导数关系： fxtfx, fff dt dx x f t f dt xtf x xt , fffffffff fffffffffff x xx f tx f fx x f t f fx xt f tt f dt x f d f dt df f dt t f d dt df f dt df f dt df dt fffd x xttxxxxtt xxxtxtxttxtt x x x x txt 2 )()( )( )( 22 2222 济南大学控制科学与工程学院计算机仿真技术16

13、龙格-库塔法 将各阶导数代入方程（），得到（）： )2( ! 3 1 )( ! 2 1 2232 1 fffffffffhfffhhfxx xttxxxxttxtnn (3.16) 假设问题的数值解公式为： 1 1 1 1 , i j jijnini r i iinn KaxhcthfK KWxx (3.17) Wi: 待定的权因子， r：解公式的阶数， Ki：导数和步长的成绩， ci, aij: 待定系数，c1=0, i=2,r 当r=1时，一阶泰勒展开式为：),( 1nnnn xthfxx 由式（）得：),( 11nnnn xthfWxx 1 1 W 济南大学控制科学与工程学院计算机仿真技

14、术17 龙格-库塔法 12122 1 22111 , , KaxhcthfK xthfK KWKWxx nn nn nn 由式（）得： 当r=2时，二阶泰勒展开式为： )( 2 1 2 1 fffhhfxx xtnn ffWhafWhchfWWxx xtnn2 2 212 2 2211 将K2在点（tn,xn）附近泰勒展开为： xtxtxtnn fhKafhchfKa x f hhc t f hxthfK nnnn 121 2 2121),(2),(2 , 将K1和K2代入xn+1中，整理得： 与泰勒展开式比较，得到： 2 1 , 2 1 , 1 2212221 WaWcWW 取c2=1，得到

15、二阶龙格-库塔计算公式： 121 211 , 5 . 05 . 0 KxhthfKxthfK KKxx nnnn nn (3.17) 济南大学控制科学与工程学院计算机仿真技术18 龙格-库塔法 由于泰勒展开式取h，h2项，h3及以上高阶项忽略，因此二阶龙格-库塔法的截断误差正比 于h3. 121 211 , 5 . 05 . 0 KxhthfKxthfK KKxx nnnn nn二阶龙格-库塔公式： 23 12 1 3211 3 2 , 3 2 3 1 , 3 1 , )3( 4 1 KxhthfK KxhthfK xthfK KKKxx nn nn nn nn 三阶龙格-库塔公式： (3.1

16、8) 济南大学控制科学与工程学院计算机仿真技术19 四阶龙格-库塔法 34 23 12 1 43211 , 2 1 , 2 1 2 1 , 2 1 , )22( 6 1 KxhthfK KxhthfK KxhthfK xthfK KKKKxx nn nn nn nn nn 四阶龙格-库塔公式： (3.19) 四阶龙格-库塔公式的截断误差正比于h5，可以满足大部分实际工程问题，计算精度较高， 编程容易，算法稳定，可以自启动，在系统仿真中应用最为广泛. 济南大学控制科学与工程学院计算机仿真技术20 四阶龙格-库塔法 例3.2 对例的系统用显式四阶龙格-库塔法求t=1时的数值解. 系统方程为仿真步长

17、h=0.1. 100)( 2 ytyt y 方程可以写为： )(,)( 2 tyytfty 显式四阶龙格-库塔法递推公式为： 2 334 2 223 2 112 2 1 43211 , 2 1 2 1 , 2 1 2 1 2 1 , 2 1 , )22( 6 1 KyhKyhthfK KyhKyhthfK KyhKyhthfK hyythfK KKKKyy nnn nnn nnn nnn nn 济南大学控制科学与工程学院计算机仿真技术21 四阶龙格-库塔法 1, 0 00 yt 5000. 0)22( 6 1 , 0 . 1 9091. 0)22( 6 1 0826. 01 . 0,09118

18、. 05 . 01 . 0 09025. 05 . 01 . 0, 1 . 01 . 0, 1 . 0 432191010 432101 2 304 2 203 2 102 2 011 KKKKyyt KKKKyy KyKKyK KyKyKt 济南大学控制科学与工程学院计算机仿真技术22 龙格-库塔法的稳定区域 将x(tn+h)在t=tn点进行泰勒展开： nnnnn txhtxhtxhtxhtx 32 ! 3 1 ! 2 1 )()(txtx0j 同时方程的解为： 龙格-库塔法取r阶，对于泰勒展开取r阶项，将方程解代入，可以得到: n r nnnn xh r xhhxxx)( ! 1 )( !

19、 2 1 2 1 r阶龙格-库塔法的稳定条件为： 1)( ! 1 )( ! 2 1 1 2 r h r hh 取 ： hh 1 ! 1 ! 2 1 1 2 r h r hh (3.20) 济南大学控制科学与工程学院计算机仿真技术23 龙格-库塔法的稳定区域 该表格给出各阶龙格-库塔公式的稳定条件： r1稳定区域 1(-2, 0) 2(-2, 0) 3(-2.51, 0) 4(-2.78, 0) h1 2 2 1 1hh 32 6 1 2 1 1hhh 432 24 1 6 1 2 1 1hhhh 选择的步长应落在稳定区域内，从落在稳定区域可以看出，系 统的特征根越大（系统时间常数越小），需要的

20、积分步长就越小. h hh 济南大学控制科学与工程学院计算机仿真技术24 3.1.3 线性多步法 以上所述的数值解法均为单步法，只要知道xn, f(tn,xn)的值就可以递推出xn+1，也就是说， 根据初始值可以递推出相继各时刻的x值，所以这种方法都可以自启动. 多步法所不同的是，需要x和f(t,x)在tn, tn-1, tn-2 各个时刻的值. 显然多步法计算公式不能自 启动，并且在计算过程中占用的内存较大，但可以提高计算精度和速度. 1. 亚当斯-贝希霍斯显式多步法 简称亚当斯法（Adams）. 对于式(3.3)： 利用插值多项式来近似代替f(t,x). 取tn和tn以前的k个节点，用这些

21、节点的函数值fn, fn-1, , fn- k+1构造一个插值多项式Pk,n来近似fn, k为多项式系数. 根据牛顿后插公式，可以得到： (3.3) )(,( 1nnnn txthftxtx 济南大学控制科学与工程学院计算机仿真技术25 亚当斯-贝希霍斯显式多步法 n k k knnn n n nnk f kh tttttt f h tt ftP 1 1 1 , )!1( )()( )( (3.21) 1 22 1 11 11 2 1 0 n k n k nn k n k nnnnn nnn nn fffff fffff fff ff 设t-tn=sh，可以将多项式P推导成如下格式： n k knn k i n i ink fffftP 1 1 1 1 0 0 1 0 , )( 1 0 0 ) 1() 1( ! 1 1 dsisss i i (3.22) 济南大学控制科学与工程学院计算机仿真技术26 亚当斯-贝希霍斯显式多步法 (3.23) 因此，可以得到亚当斯多步法的计算公式： 1 0 1 k i n i inn fhxx 当k=1时，得到一阶的欧拉公式： nnnn hfxfhxx 0 01 当k=2时，得到二阶的亚当斯计算公式： ffhxx n