《计算仿真课件第三章》由会员分享，可在线阅读，更多相关《计算仿真课件第三章（29页珍藏版）》请在金锄头文库上搜索。
1、计算机仿真技术授课教师：李实1教1007室1计算仿真课件第三章济南大学控制科学与工程学院计算机仿真技术第三章 数字仿真通用算法 在工程实际和科学研究中所遇到的问题往往很复杂，很多情况下无法给出描述动态特性的微分方程解的解析表达式，多数只能用近似的数值方法求解. 随着计算机软硬件和数值理论的进展，微分方程的数值解方法已成为当今研究、分析和设计系统的一种有力工具. 3.1 系统仿真中常用的数值积分法 3.2 刚性系统的特点及算法 3.3 实时半实物仿真 3.4 采样控制系统的仿真方法 3.5 分布参数系统的数字仿真2济南大学控制科学与工程学院计算机仿真技术3.1 系统仿真中常用的数值积分法3.1.
2、1 欧拉法和改进欧拉法3.1.2 龙格-库塔法3.1.3 线性多步法3.1.4 MATLAB中常微分方程求解方法3济南大学控制科学与工程学院计算机仿真技术3.1.1 欧拉法和改进欧拉法 欧拉法是最简单的单步法，它是一阶的，精度较差. 但公式简单，有明显的几何意义，适合初学者在直观上学习数值x(tn)是如何逼近微分方程的精确解x(t)的. 考虑初值问题：(3.1) 式 （ ） 的 解 x(t)是 一 连 续 变 量 x的 函 数 ， 现 在 要 用 一 系 列 离 散 时 刻 的 近 似 值x(t1),x(t2),x(tn)来代替，其中ti=t0+ih, h成为步长，是相邻两点的距离. 将式（）
3、在区间（ti,ti+1）上积分，则可得：(3.2) 积分项的几何意义是曲线f(x,y)在区间（ti,ti+1）上的面积，当（ti,ti+1）充分小时，可以用矩形面积来代替：4济南大学控制科学与工程学院计算机仿真技术欧拉法因此，式（）可近似为：写成递推式：(3.3) 已知x(0)=x0, 可以由上式求出x(t1), x(t2), . 这种算法成为单步法. 可以直接由初始值递推出其它各时间的值，因此单步法是一种自启动算法.定义x(tn)=xn, 式子（）还可以写成：(3.4)该式称为前向欧拉公式，又称显式欧拉公式.后向欧拉公式（隐式欧拉公式）：(3.5)5济南大学控制科学与工程学院计算机仿真技术欧
4、拉法例3.1 设系统方程为试用欧拉法求其在t=1时的数值解，仿真步长，要求分别使用前向欧拉法和后向欧拉法.方程可以写为：写出该系统向欧拉递推公式为：6济南大学控制科学与工程学院计算机仿真技术欧拉法该系统的后向欧拉递推公式为：整理后得到：将代入，求得方程的解为：7济南大学控制科学与工程学院计算机仿真技术欧拉法的几何意义htx(t)tnxntn+1xn+1前向欧拉法前向欧拉法tn+1htx(t)tnxnxn+1后向欧拉法后向欧拉法8济南大学控制科学与工程学院计算机仿真技术改进欧拉法（预测-校正法）对积分公式（）利用梯形面积公式计算积分项，得：(3.2)(3.6)写成递推差分格式：(3.7)可以先用
5、欧拉法计算xn+1, 定义为预估值，写为xpn+1，再将该预估值代入原方程中计算函数，最后利用式（）得到修正后的xcn+1, 改进后的欧拉法描述如下：(3.8)9济南大学控制科学与工程学院计算机仿真技术欧拉法：误差分析欧拉法以一定的步长来进行数值计算，所得到的解在tn点的近似值与微分方程的精确解存在误差.数值仿真的误差可分为截断误差和舍入误差两种. 截断误差与采用的计算方法有关，而舍入误差由计算机的字长决定.截断误差：将x(tn+h)在t=tn点进行泰勒展开：(3.8)取并截断，得到欧拉公式，Rn为截断误差，与h2成正比(3.9)解从t=0到t=tn所积累的误差为整体误差，比局部误差要大，欧拉
6、法的整体误差与h成正比，即为O1(h).10济南大学控制科学与工程学院计算机仿真技术欧拉法：误差分析舍入误差：由于计算机进行计算时，数的位数有限引起的，一般舍入误差与h-1成正比，可以得到欧拉法总误差为：(3.10)由式（）可以看出，步长h增加，截断误差增加，反之，步长h减小，舍入误差加大.11济南大学控制科学与工程学院计算机仿真技术欧拉法：稳定性分析求解微分方程的另一个重要问题是数值解是否稳定. 假设方程的特征根为. 方程可以写成：(3.10)若该方程稳定，需要满足前向欧拉法的稳定性：前向欧拉法公式：该差分方程稳定的条件：(3.11)或：(3.12)选择的步长过大会使算法变得不稳定12济南大
7、学控制科学与工程学院计算机仿真技术欧拉法：稳定性分析后向欧拉法的稳定性：后向欧拉法公式：该差分方程稳定的条件：(3.13)因为为负，所以（）必然成立，说明后向欧拉法是恒稳定的.改进欧拉法（梯形公式）的稳定性：改进欧拉法公式：该差分方程稳定的条件：(3.14)同理，改进欧拉法是恒稳定的.13济南大学控制科学与工程学院计算机仿真技术3.1.2 龙格-库塔法将x(tn+h)在t=tn点进行泰勒展开：(3.15)欧拉法是截去h2以后各项所得到的一阶一步法，因此精度较低. 如果将泰勒展开式多选几项以后截断，可以得到精度较高的高阶数值解，但是直接使用泰勒展开要计算函数的高阶导数. 龙格-库塔法采用间接利用
8、泰勒展开式的思路，用在n个点上的函数值f的线性组合来代替f的导数，然后按泰勒展开式确定其中的系数，以提高算法的精度. 这样既避免计算函数的导数，同时又保证了计算精度.由于龙格-库塔法具有许多优点，因此在许多仿真程序包中，它是一个最基本的算法之一.14济南大学控制科学与工程学院计算机仿真技术龙格-库塔法将x(tn+h)在t=tn点进行泰勒展开：(3.15)根据偏导数关系：15济南大学控制科学与工程学院计算机仿真技术龙格-库塔法将各阶导数代入方程（），得到（）：(3.16)假设问题的数值解公式为：(3.17)Wi: 待定的权因子，r：解公式的阶数，Ki：导数和步长的成绩，ci, aij: 待定系数
9、，c1=0, i=2,r当r=1时，一阶泰勒展开式为：由式（）得：16济南大学控制科学与工程学院计算机仿真技术龙格-库塔法由式（）得：当r=2时，二阶泰勒展开式为：将K2在点（tn,xn）附近泰勒展开为：将K1和K2代入xn+1中，整理得：与泰勒展开式比较，得到：取c2=1，得到二阶龙格-库塔计算公式：(3.17)17济南大学控制科学与工程学院计算机仿真技术龙格-库塔法由于泰勒展开式取h，h2项，h3及以上高阶项忽略，因此二阶龙格-库塔法的截断误差正比于h3.二阶龙格-库塔公式：三阶龙格-库塔公式：(3.18)18济南大学控制科学与工程学院计算机仿真技术四阶龙格-库塔法四阶龙格-库塔公式：(3
10、.19)四阶龙格-库塔公式的截断误差正比于h5，可以满足大部分实际工程问题，计算精度较高，编程容易，算法稳定，可以自启动，在系统仿真中应用最为广泛.19济南大学控制科学与工程学院计算机仿真技术四阶龙格-库塔法例3.2 对例的系统用显式四阶龙格-库塔法求t=1时的数值解. 系统方程为仿真步长h=0.1.方程可以写为：显式四阶龙格-库塔法递推公式为：20济南大学控制科学与工程学院计算机仿真技术四阶龙格-库塔法21济南大学控制科学与工程学院计算机仿真技术龙格-库塔法的稳定区域将x(tn+h)在t=tn点进行泰勒展开：同时方程的解为：龙格-库塔法取r阶，对于泰勒展开取r阶项，将方程解代入，可以得到:r
11、阶龙格-库塔法的稳定条件为：取 ：(3.20)22济南大学控制科学与工程学院计算机仿真技术龙格-库塔法的稳定区域该表格给出各阶龙格-库塔公式的稳定条件：r1稳定区域1(-2, 0)2(-2, 0)3(-2.51, 0)4(-2.78, 0)选择的步长应落在稳定区域内，从落在稳定区域可以看出，系统的特征根越大（系统时间常数越小），需要的积分步长就越小.23济南大学控制科学与工程学院计算机仿真技术3.1.3 线性多步法以上所述的数值解法均为单步法，只要知道xn, f(tn,xn)的值就可以递推出xn+1，也就是说，根据初始值可以递推出相继各时刻的x值，所以这种方法都可以自启动.多步法所不同的是，需
12、要x和f(t,x)在tn, tn-1, tn-2 各个时刻的值. 显然多步法计算公式不能自启动，并且在计算过程中占用的内存较大，但可以提高计算精度和速度.1. 亚当斯-贝希霍斯显式多步法简称亚当斯法（Adams）. 对于式(3.3)：利用插值多项式来近似代替f(t,x). 取tn和tn以前的k个节点，用这些节点的函数值fn, fn-1, , fn-k+1构造一个插值多项式Pk,n来近似fn, k为多项式系数. 根据牛顿后插公式，可以得到：(3.3)24济南大学控制科学与工程学院计算机仿真技术亚当斯-贝希霍斯显式多步法(3.21)设t-tn=sh，可以将多项式P推导成如下格式：(3.22)25济
13、南大学控制科学与工程学院计算机仿真技术亚当斯-贝希霍斯显式多步法(3.23)因此，可以得到亚当斯多步法的计算公式：当k=1时，得到一阶的欧拉公式：当k=2时，得到二阶的亚当斯计算公式：其中(3.24)26济南大学控制科学与工程学院计算机仿真技术亚当斯-贝希霍斯显式多步法当k=3时，得到三阶的亚当斯计算公式：其中(3.25)27济南大学控制科学与工程学院计算机仿真技术亚当斯-莫尔顿隐式多步法取n+1时刻的点来代替n时刻的点，相应的多项式P变为：(3.26)根据显式多步法的推导过程，可以得到隐式多步伐的计算公式：(3.27)28济南大学控制科学与工程学院计算机仿真技术亚当斯预测-校正法将亚当斯的显式公式与隐式公式联合使用，前者提供预测值，后者将预测值加以校正，这就是预测-校正法. 常用的四阶亚当斯预测-校正法计算公式为：(3.28)29
