计算机仿真技术基础第3章-连续系统数值积分仿真方法学 课件

1、计算机仿真技术基础 第三章连续系统数值积分仿真方法学第二章 连续系统数值积分仿真方法学第一节 数值积分法的基本原理第二节 数值积分法的单步算法第三节 数值积分法的多步算法如何把已建立起来的数学模型转换成仿真运算模型（二次建模），以便为分析解决实际问题服务那是系统仿真学科的一个重要研究内容。对于复杂的数学模型来说，求其解析解是很烦琐和困难的，大多数情况下不求出解析解，或者根本不存在解析解，因此借助于数值解法对连续系统进行仿真研究。 用计算机不可能求出系统的解析解（连续）只能求出连续响应曲线上的有限个点，即数值解连续系统数值积分法：就是利用数值积分方法对常微分方程建立离散化形式的数学模型（差分方程

2、）并求出数值解。 最常用的数值解法有： 欧拉法、梯形法、Adams、RungeKutta法。 描述各类系统最基本的模型用微分方程或状态空间表达式，二次建模就是要求出适合用数字计算机求解的模型，就需要把微分运算转化成算术运算在用计算机求解。第一节 数值积分法的基本原理首先把需仿真研究的系统表示成一阶微分方程组或状态方程的形式。 以一阶连续系统为例，微分方程及初值如下： 设方程在处的连续解为：希望找到一个近似公式来表示方程的近似解：为精确值的近似值为准确积分值的近似值所谓微分方程初值问题的数值解法就是寻求真解在一系列离散点上的近似解（数值解）相邻两个时间离散点的间隔称为计算步距或步长，通常情况为定

3、值，也有变步长。可见，微分方程初值问题数值解法的主要问题归结为对为此，要先把连续的微分方程用数值积分法转化为离散的差分方程的初值问题，然后根据初始条件X0逐步递推计算出后续时刻的数值解：如何求出定积分的近似解数值解法的共同特点是步进式的递推算法，主要有：单步法、多步法和预估校正法，并有显式和隐式之分。一、欧拉法（Euler Method)Euler法是最简单的一种数值积分法的单步运算，虽然计算精度较差，但几何意义明显，便于理解，能说明构造数值积分算法的基本思想。下面采用三种方法推导出Euler法的数值近似公式，以便对数值积分器的基本思想能透彻了解。第二节 数值积分法的单步算法 、Taylor级

4、数展开x(t)为解析解，将x(t)展开成Taylor级数 以一阶连续系统为例，微分方程及初值如右：只取一次项，其余忽略写成差分方程为这就是解微分方程初值问题的欧拉算法。、矩形近似把积分区间 h 取得足够小，将 在 近似为常数用左矩形面积近似该区间的曲线面积对方程在上求积分也能得到：左矩形（也称为前向欧拉法）近似及误差将 在 近似为常数用右矩形面积近似该区间的曲线面积得到：这是右矩形欧拉公式，是一个隐式算法对积分右矩形（也称为后向欧拉法）近似及误差、切线近似（）取切线上处的值来近似在 的一个小邻域内，曲线x(t)可用 处的切线来表示， x(t) 在 处的斜率为：也能得到：前向欧拉法过点 以 fn

5、 为斜率的切线方程为：、切线近似（）过点 以 fn+1 为斜率的切线方程为：取切线上处的值来近似在 曲线 x(t) 可用 处的切线来表示， x(t) 在 处的斜率为：也能得到：后向欧拉法欧拉法（切线推导）的几何意义欧拉法实际计算时的几何意义例：设系统方程试用Euler法求其数值解，取步长 h=0.1,解:前向Euler法递推式：有初始条件：可进行递推：后向Euler法递推式：隐式算法，需先解此非线性方程：由此公式可进行递推。前向Euler法与精确值比较前向Euler法、后向Euler法与精确值比较前向Euler法在不同步长的结果比较二、梯形法 Euler法的计算精度较差，如果改用梯形面积代替每

6、个步距的曲线面积，就可提高精度。精确积分应为曲边梯形的面积：现用直边梯形的面积来近似：写成差分方程为梯形近似及其误差梯形法实质是采用了两点斜率平均值的结果，由于利用了两点的信息，从而提高了计算精度。和 这一思想被广泛地应用于许多算法中，实际计算时，如果在每个积分步矩中多取几个点，分别求出其斜率，然后取不同的权值为： 后面RungeKutta法就是采用这样的思想来进行计算的。梯形法的几何意义也可按折线理解， 梯形法大大提高了精度，但为隐式算法，每次递推计算时需解一次非线性方程，计算量较大由此考虑进行改进，先用Euler法计算出：的近似值代入导函数求出近似值再代入梯形公式求解。预估公式（Euler

7、法）校正公式（梯形法） 为预估校正法，也称为改进的Euler法。例：设系统方程试用梯形法求其数值解，取步长h=0.1,解:梯形法递推式：隐式算法，需先解此非线性方程：由此公式和初始条件可进行递推，见仿真结果。例：设系统方程用改进欧拉法求数值解，取步长h=0.1,解:改进欧拉法：由此公式和初始条件可进行递推，见仿真结果。前向Euler法、梯形法与精确值比较前向Euler法、改进欧拉法与精确值比较梯形法、改进Euler法与精确值比较三、Runge-Kutta法 、Taylor级数匹配原理由于输入函数是 t 的函数，则将记做 得微分方程：如果对变量t、x具有各阶导数，可推得x(t)的各阶导数。 设已

8、知 进行Taylor级数展开： 若已知 的值，则当 h 较小时，可用级数展开的前 p+1项作为近似,令 则 即 以上公式(1)就称为p阶的Taylor展开法递推公式 之间的误差为： 局部截断误差与hp+1是同阶无穷小量,记为O(hp+1) (1)欧拉法的Taylor级数展开只取一次项，其余忽略写成差分方程为这就是解微分方程初值问题的欧拉算法。所以欧拉法称为一阶的Taylor展开法递推公式 局部截断误差与h2是同阶无穷小量,记为O(h2)梯形法的Taylor级数展开取一次项和二次项，写成差分方程为所以梯形法称为二阶的Taylor展开法递推公式O(h3) 可见Taylor展开法需用 在 的高阶导数

9、计算 ，不便于数值计算。2、 Runge-Kutta法 于是Runge-Kutta法用 在一些点上的值表示 ，使局部截断误差的阶数与Taylor展开法相等。避免了求高阶导数，又保证高的精度。 对微分方程在区间的连续解为：在区间 取m 个点若已知则用它们的一次组合去近似 ，即？现在的问题是如何求设已知因为 是未知的，最简单的用欧拉法构造。由Euler法若按原有欧拉法计算下一步距的值得以二阶Runge-Kutta法为例说明。精度较差，为改进精度，由欧拉法以步长 ah得到另一个点设已知某步从点开始，沿斜率为方向移动一个步长h，得到从图中可见三点可得到比欧拉法精度较好的近似值。在一条直线上，可选取参数

10、 b1、 b构造：二阶R-K构造法tX(t)tmtm+1xEm+1x(tm+1)tm+ahxEm+ax(tm)欧拉法计算以步长ah得到另一个点从点开始，沿斜率为方向移动一个步长h，得到可选取参数 b1、 b构造：以得到比欧拉法精度较好的近似值。以上计算归结为：记则？如何选取参数a、b1、b2，可获得最高的精度将k2在处展开成Taylor级数见代入式子见若满足则与Taylor展开式前三项相同局部截断误差三个未知数，两个方程，有多组解：局部截断误差是h的三阶无穷小量,比欧拉法的精度高一阶若取a=1时，得到：改进的Euler公式若取时，得：修正的Euler公式或中点公式显式p阶Runge-Kutta

11、法的一般形式为：3、常用的Runge-Kutta法：） Kutta三阶法） Heun三阶法3） 经典显式四阶Runge-Kutta法四、微分方程数值积分的矩阵分析法 前述的各类数值积分公式都以一阶系统（单个的微分方程）进行讨论，而实际工程中大量的仿真对象是高阶系统，可用一阶微分方程组来描述。此时，数值积分公式有相应的矩阵形式。矩阵形式的数值积分公式：） 欧拉法公式前向欧拉法公式后向欧拉法公式） 梯形法公式）二阶龙格库塔法公式改进的欧拉法公式，是预估校正公式。） 四阶龙格库塔法公式（RK4）对于 n 阶系统，状态向量 x 为 n 维，计算中每前进一步 h ，要计算 4n 个 kij 值，对状态空

12、间表达式：此时，RK4公式的个 k 值：例：系统方程取步长 h=0.1,试用RK4法求t=0.1，0.2时的解解:将原系统方程化为状态方程形式：见仿真结果作业：P1493.2习题3-2：已知用前向欧拉法、梯形法求其数值解，取步长h=0.1解:前向欧拉法递推式：习题3-2：已知用前向欧拉法、梯形法求其数值解，取步长h=0.1解:梯形法递推式：隐式算法，需先解此非线性方程：单步法的特点：计算 n+1 时刻的值 yn+1 时，只用到第 n 时刻的 yn 和 fn 。如果能利用多步计算信息（历史时刻值），则可能既加快仿真速度又获得较高的仿真精度，这就是构造多步法的出发点。第三节 数值积分法的多步算法实

13、际在逐步递推过程中，计算 yn+1 时已经获得一系列的近似值： 以及 。多步法中以 Adams 法最具代表性，应用最为普遍。 对一阶连续系统：连续解为：现过三点按插值原理构造一个多项式 来逼近函数对函数,再对多项式积分近似积分 一、Adams算法得：多项式 中的系数由下决定： 拉格朗日插值公式令：同时考虑：因为有：进行变量替换：显然，对多项式的积分计算很容易。微分方程连续解为：写成差分方程：这就是显式两步二阶Adams递推式。显式 Adams 算法的系数值显式 Adams 算法的递推公式为： bibob1b2b3b40113/2-1/2223/12-16/125/12325/24-59/243

14、7/24-9/2441901/720-2774/7202616/720-1274/720-19/720k隐式 Adams 算法的系数值隐式 Adams 算法的递推公式为： bib-1b0b1b2b30111/21/225/128/12-1/1239/2419/24-5/241/244251/720646/720-264/720106/720-19/720k无论用显式或隐式 k 阶Adams 法求解微分方程初值问题数值，需要先知道 k+1 个初始值。例如三步Adams 法：于是初始值只能从初始条件得到 ，还需知道才能求出：需用单步法求出，才能使多步法的递推计算能够进行。为保证多步法的精度，注意选

15、择相应精度的单步法计算初始值。例：设系统方程用显式二阶Adams法求解，取步长h=0.1,解:显式二阶Adams法：起步初始值由梯形公式求出：下面就可以用Adams公式进行递推：有初始条件：可进行初值计算： 为计算 的值，用到 时刻以前的值来推导，可获得性能更好的算法。由 Adams 法得出更为一般的形式：即： 再令 二、线性多步法解出： 为计算 的值，用到和相应的导数值而且公式关于 x 、 f 是线性的，称为线性 k 步法来推导即： 当： 称为显式线性 k 步法当： 称为隐式线性 k 步法 时的情形。是线性 k 步法再看前面的 当： 显式 Adams 算法的递推公式为：隐式 Adams 算法

16、的递推公式为：？接下来的问题就是线性 k 步法如何确定其中的常数 应用 Taylor 级数匹配原理，使局部截断误差尽可能的小。将处展开成Taylor 级数，和再带入误差公式得：其中：对于若能选取 使：则此线性多步法为 p 阶 k 步算法。3、构造线性多步法 确定步数 k ，由Taylor级数匹配原理的 得到方程求出待定系数 ，得出尽可能高阶的算法。以二步法为例说明构造步骤：5个未知数，4个方程， 得： 得到一般的线性二步法形式：是三阶二步法； 又可得：是四阶二步法：称为Milne算法隐式算法为三阶隐式算法就是三阶隐式Adams法?如何得到二步显式算法得到二步显式算法但这是一个数值不稳定算法，计

17、算中出现的微小误差会迅速增长。因此，二步显式算法无法达到三阶，最多达到二阶。只用三个方程二阶二步显式算法的一般形式：二阶二步Adams显式算法为Adams法为向后微分公式显然是隐式公式若要求向后微分公式为K阶算法，则当K=1为向后Eular法K=1，2，.6的向后微分公式系数值 11-1122/31/3-4/3136/11-2/119/11-18/111412/253/25-16/2536/25-48/2515161k 计算 时，仅已知和相应的导数值 称为线性多步法的起步值，需用单步法来求出。 若线性多步法是P阶算法，计算起步值的算法不应低于P阶，否则影响计算精度。仿真步长 h 的选取，是否影

18、响仿真结果？先看一个例子：第四节 数值积分法稳定性分析系统稳定，精确解：用前向欧拉法求解： ch3_6.m由此可见，仿真步长 h 的选取，会影响仿真结果。 用前向欧拉法求解，当 h 0.2时不稳定，是由于步长太大，从而截断误差太大造成的。一、数值解法稳定性的含义 数值解的稳定性：在扰动（初始误差、舍入误差、截断误差）的影响下，计算过程中的累积误差不会随计算步数的增加而无限增长。 微分方程的数值积分方法，实质是微分方程的差分化，然后从初始条件递推迭代。不同数值解法对应着不同的差分方程，是否稳定取决于该差分方程的特征根是否满足稳定性要求。（处于Z平面上以原点为圆心的单位圆内）着重研究单步法的稳定性对步长的限制。二、数值解法稳定性分析这样做的根据是：1）试验模型简单，对其数值不稳定的方法，不可用；2）一般的初始问题在其解的存在