仿真数值积分法

仿真数值积分法第1页，课件共104页，创作于2023年2月第三章连续系统数值积分仿真方法学第一节数值积分法的基本原理第二节数值积分法的单步算法第三节数值积分法的多步算法第2页，课件共104页，创作于2023年2月

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

欧拉法、梯形法、Adams、Runge—Kutta法。描述各类系统最基本的模型用微分方程或状态空间表达式，二次建模就是要求出适合用数字计算机求解的模型，就需要把微分运算转化成算术运算在用计算机求解。第4页，课件共104页，创作于2023年2月

第一节数值积分法的基本原理 首先把需仿真研究的系统表示成一阶微分方程组或状态方程的形式。以一阶连续系统为例，微分方程及初值如下：在时间点处的解为：第5页，课件共104页，创作于2023年2月希望找到一个近似公式来表示方程的近似解：为精确值的近似值为准确积分值的近似值所谓微分方程初值问题的数值解法就是寻求真解在一系列离散点上的近似解（数值解）相邻两个时间离散点的间隔称为计算步距或步长，通常情况为定值，也有变步长。第6页，课件共104页，创作于2023年2月第7页，课件共104页，创作于2023年2月

可见，微分方程初值问题数值解法的主要问题归结为对 为此，要先把连续的微分方程用数值积分法转化为离散的差分方程的初值问题，然后根据初始条件X0逐步递推计算出后续时刻的数值解：如何求出定积分的近似解 数值解法的共同特点是步进式的递推算法，主要有：单步法、多步法和预估－校正法，并有显式和隐式之分。第8页，课件共104页，创作于2023年2月

一、欧拉法（EulerMethod) Euler法是最简单的一种数值积分法的单步运算，虽然计算精度较差，但几何意义明显，便于理解，能说明构造数值积分算法的基本思想。 下面采用三种方法推导出Euler法的数值近似公式，以便对数值积分器的基本思想能透彻了解。

第二节数值积分法的单步算法第9页，课件共104页，创作于2023年2月１、Taylor级数展开x(t)为解析解，将x(t)展开成Taylor级数以一阶连续系统为例，微分方程及初值如右：只取一次项，其余忽略写成差分方程为这就是解微分方程初值问题的欧拉算法。第10页，课件共104页，创作于2023年2月

２、矩形近似把积分区间h取得足够小，将在近似为常数用左矩形面积近似该区间的曲线面积对方程在上求积分也能得到：第11页，课件共104页，创作于2023年2月左矩形（也称为前向欧拉法）近似及误差第12页，课件共104页，创作于2023年2月将在近似为常数用右矩形面积近似该区间的曲线面积得到：这是右矩形欧拉公式，是一个隐式算法．对积分第13页，课件共104页，创作于2023年2月右矩形（也称为后向欧拉法）近似及误差第14页，课件共104页，创作于2023年2月

３、切线近似（１）取切线上处的值来近似在的一个小邻域内，曲线x(t)可用处的切线来表示，x(t)在处的斜率为：也能得到：前向欧拉法过点以fn

为斜率的切线方程为：第15页，课件共104页，创作于2023年2月

３、切线近似（２）过点以fn+1为斜率的切线方程为：取切线上处的值来近似在曲线x(t)可用处的切线来表示，x(t)在处的斜率为：也能得到：后向欧拉法第16页，课件共104页，创作于2023年2月欧拉法（切线推导）的几何意义第17页，课件共104页，创作于2023年2月欧拉法实际计算时的几何意义第18页，课件共104页，创作于2023年2月第19页，课件共104页，创作于2023年2月

例：设系统方程试用Euler法求其数值解，取步长h=0.1,解:前向Euler法递推式：有初始条件：可进行递推：第20页，课件共104页，创作于2023年2月后向Euler法递推式：隐式算法，需先解此非线性方程：由此公式可进行递推。第21页，课件共104页，创作于2023年2月前向Euler法与精确值比较第22页，课件共104页，创作于2023年2月前向Euler法、后向Euler法与精确值比较第23页，课件共104页，创作于2023年2月第24页，课件共104页，创作于2023年2月前向Euler法在不同步长的结果比较第25页，课件共104页，创作于2023年2月第26页，课件共104页，创作于2023年2月

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

例：设系统方程试用梯形法求其数值解，取步长h=0.1,解:梯形法递推式：隐式算法，需先解此非线性方程：由此公式和初始条件可进行递推，见仿真结果。第34页，课件共104页，创作于2023年2月

例：设系统方程用改进欧拉法求数值解，取步长h=0.1,解:改进欧拉法：由此公式和初始条件可进行递推，见仿真结果。第35页，课件共104页，创作于2023年2月前向Euler法、梯形法与精确值比较第36页，课件共104页，创作于2023年2月前向Euler法、改进欧拉法与精确值比较第37页，课件共104页，创作于2023年2月梯形法、改进Euler法与精确值比较第38页，课件共104页，创作于2023年2月

三、Runge-Kutta法１、Taylor级数匹配原理由于输入函数是t的函数，则将记做得微分方程：如果对变量t、x具有各阶导数，可推得x(t)的各阶导数。第39页，课件共104页，创作于2023年2月

设已知

进行Taylor级数展开：

若已知的值，则当h较小时，可用级数展开的前p+1项作为近似,令第40页，课件共104页，创作于2023年2月

则

即

以上公式(1)就称为p阶的Taylor展开法递推公式

之间的误差为：

局部截断误差与hp+1是同阶无穷小量,记为O(hp+1)

(1)第41页，课件共104页，创作于2023年2月欧拉法的Taylor级数展开只取一次项，其余忽略写成差分方程为这就是解微分方程初值问题的欧拉算法。所以欧拉法称为一阶的Taylor展开法递推公式

局部截断误差与h2是同阶无穷小量,记为O(h2)第42页，课件共104页，创作于2023年2月梯形法的Taylor级数展开取一次项和二次项，写成差分方程为所以梯形法称为二阶的Taylor展开法递推公式O(h3)第43页，课件共104页，创作于2023年2月可见Taylor展开法需用在的高阶导数计算，不便于数值计算。2、Runge-Kutta法于是Runge-Kutta法用在一些点上的值表示，使局部截断误差的阶数与Taylor展开法相等。避免了求高阶导数，又保证高的精度。第44页，课件共104页，创作于2023年2月对微分方程在区间的连续解为：在区间取m个点若已知则用它们的一次组合去近似，即第45页，课件共104页，创作于2023年2月？现在的问题是如何求设已知因为是未知的，最简单的用欧拉法构造。由Euler法第46页，课件共104页，创作于2023年2月若按原有欧拉法计算下一步距的值得以二阶Runge-Kutta法为例说明。精度较差，为改进精度，由欧拉法以步长ah得到另一个点设已知某步从点开始，沿斜率为方向移动一个步长h，得到第47页，课件共104页，创作于2023年2月从图中可见三点可得到比欧拉法精度较好的近似值。在一条直线上，可选取参数b1、b２构造：第48页，课件共104页，创作于2023年2月二阶R-K构造法tX(t)tmtm+1xEm+1x(tm+1)tm+ahxEm+ax(tm)第49页，课件共104页，创作于2023年2月欧拉法计算以步长ah得到另一个点

从点开始，沿斜率为方向移动一个步长h，得到可选取参数b1、b２构造：以得到比欧拉法精度较好的近似值。第50页，课件共104页，创作于2023年2月以上计算归结为：记则第51页，课件共104页，创作于2023年2月？如何选取参数a、b1、b2，可获得最高的精度将k2在处展开成Taylor级数见第52页，课件共104页，创作于2023年2月代入式子见若满足则与Taylor展开式前三项相同第53页，课件共104页，创作于2023年2月局部截断误差三个未知数，两个方程，有多组解：局部截断误差是h的三阶无穷小量,比欧拉法的精度高一阶若取a=1时，得到：改进的Euler公式第54页，课件共104页，创作于2023年2月若取时，得：修正的Euler公式或中点公式显式p阶Runge-Kutta法的一般形式为：第55页，课件共104页，创作于2023年2月3、常用的Runge-Kutta法：１）Kutta三阶法第56页，课件共104页，创作于2023年2月２）Heun三阶法第57页，课件共104页，创作于2023年2月3）经典显式四阶Runge-Kutta法第58页，课件共104页，创作于2023年2月

四、微分方程数值积分的矩阵分析法

前述的各类数值积分公式都以一阶系统（单个的微分方程）进行讨论，而实际工程中大量的仿真对象是高阶系统，可用一阶微分方程组来描述。此时，数值积分公式有相应的矩阵形式。第59页，课件共104页，创作于2023年2月矩阵形式的数值积分公式：１）欧拉法公式前向欧拉法公式后向欧拉法公式第60页，课件共104页，创作于2023年2月２）梯形法公式３）二阶龙格－库塔法公式改进的欧拉法公式，是预估－校正公式。第61页，课件共104页，创作于2023年2月４）四阶龙格－库塔法公式（RK4）第62页，课件共104页，创作于2023年2月对于n阶系统，状态向量x为n维，计算中每前进一步h，要计算4n个kij值，对状态空间表达式：此时，RK4公式的４个k值：第63页，课件共104页，创作于2023年2月例：系统方程取步长h=0.1,试用RK4法求t=0.1，0.2时的解解:将原系统方程化为状态方程形式：第64页，课件共104页，创作于2023年2月第65页，课件共104页，创作于2023年2月见仿真结果作业：P1493.2第66页，课件共104页，创作于2023年2月第67页，课件共104页，创作于2023年2月 单步法的特点：计算n+1时刻的值yn+1时，只用到第n时刻的yn和fn。 如果能利用多步计算信息（历史时刻值），则可能既加快仿真速度又获得较高的仿真精度，这就是构造多步法的出发点。

第三节数值积分法的多步算法 实际在逐步递推过程中，计算yn+1时已经获得一系列的近似值：以及。 多步法中以Adams法最具代表性，应用最为普遍。第68页，课件共104页，创作于2023年2月对一阶连续系统：连续解为：现过三点按插值原理构造一个多项式来逼近函数对函数,再对多项式积分近似积分

一、Adams算法第69页，课件共104页，创作于2023年2月第70页，课件共104页，创作于2023年2月得：多项式中的系数由下决定：拉格朗日插值公式令：同时考虑：第71页，课件共104页，创作于2023年2月因为有：进行变量替换：第72页，课件共104页，创作于2023年2月显然，对多项式的积分计算很容易。第73页，课件共104页，创作于2023年2月微分方程连续解为：写成差分方程：这就是显式两步二阶Adams递推式。第74页，课件共104页，创作于2023年2月显式Adams算法的系数值显式Adams算法的递推公式为：bibob1b2b3b40113/2-1/2223/12-16/125/12325/24-59/2437/24-9/2441901/720-2774/7202616/720-1274/720-19/720k第75页，课件共104页，创作于2023年2月隐式Adams算法的系数值隐式Adams算法的递推公式为：bib-1b0b1b2b30111/21/225/128/12-1/1239/2419/24-5/241/244251/720646/720-264/720106/720-19/720k第76页，课件共104页，创作于2023年2月

无论用显式或隐式k阶Adams法求解微分方程初值问题数值，需要先知道k+1个初始值。例如三步Adams法：于是初始值只能从初始条件得到，还需知道才能求出：需用单步法求出，才能使多步法的递推计算能够进行。为保证多步法的精度，注意选择相应精度的单步法计算初始值。第77页，课件共104页，创作于2023年2月

例：设系统方程用显式二阶Adams法求解，取步长h=0.1,解:显式二阶Adams法：起步初始值由梯形公式求出：第78页，课件共104页，创作于2023年2月下面就可以用Adams公式进行递推：有初始条件：可进行初值计算：第79页，课件共104页，创作于2023年2月为计算的值，用到时刻以前的值来推导，可获得性能更好的算法。由Adams法得出更为一般的形式：即：再令

二、线性多步法解出：第80页，课件共104页，创作于2023年2月为计算的值，用到和相应的导数值而且公式关于x、f是线性的，称为线性k步法来推导即：当：称为显式线性k步法当：称为隐式线性k步法第81页，课件共104页，创作于2023年2月时的情形。是线性k步法再看前面的当：显式Adams算法的递推公式为：隐式Adams算法的递推公式为：第82页，课件共104页，创作于2023年2月？接下来的问题就是线性k步法如何确定其中的常数

应用Taylor级数匹配原理，使局部截断误差尽可能的小。第83页，课件共104页，创作于2023年2月将处展开成Taylor级数，和再带入误差公式得：其中：第84页，课件共104页，创作于2023年2月对于若能选取使：则此线性多步法为p阶k步算法。第85页，课件共104页，创作于2023年2月3、构造线性多步法

确定步数k，由Taylor级数匹配原理的

得到方程求出待定系数，得出尽可能高阶的算法。以二步法为例说明构造步骤：5个未知数，4个方程，第86页，课件共104页，创作于2023年2月

得：

得到一般的线性二步法形式：是三阶二步法；

又可得：是四阶二步法：称为Milne算法隐式算法第87页，课件共104页，创作于2023年2月为三阶隐式算法就是三阶隐式Adams法?如何得到二步显式算法得到二步显式算法但这是一个数值不稳定算法，计算中出现的微小误差会迅速增长。因此，二步显式算法无法达到三阶，最多达到二阶。第88页，课件共104页，创作于2023年2月只用三个方程二阶二步显式算法的一般形式：二阶二步Adams显式算法第89页，课件共104页，创作于2023年2月为Adams法为向后微分公式显然是隐式公式若要求向后微分公式为K阶算法，则第90页，课件共104页，创作于2023年2月当K=1为向后Eular法第91页，课件共104页，创作于2023年2月K=1，2，….6的向后微分公式系数值

11-1122/31/3-4/3136/11-2/119/11-18/111412/253/25-16/2536/25-48/2515161k第92页，课件共104页，创作于2023年2月

计算时，仅已知和相应的导数值称为线性多步法的起步值，需用单步法来求出。若线性多步法是P阶算法，计算起步值的算法不应低于P阶，否则影响计算精度。第93页，课件共104页，创作于2023年2月仿真步长h的选取，是否影响仿真结果？先看一个例子：

第四节数值积分法稳定性分析系统稳定，精确解：用前向欧拉法求解：ch3_6.m第94页，课件共104页，创作于20