仿真_3_数值积分法.ppt

主讲教师 姜萍 第三章连续系统数值积分仿真方法学 第三章连续系统数值积分仿真方法学 第一节数值积分法的基本原理 第二节数值积分法的单步算法 第三节数值积分法的多步算法 如何把已建立起来的数学模型转换成仿真运算模型 二次建模 以便为分析解决实际问题服务那是系统仿真学科的一个重要研究内容 对于复杂的数学模型来说 求其解析解是很烦琐和困难的 大多数情况下不求出解析解 或者根本不存在解析解 因此借助于数值解法对连续系统进行仿真研究 用计算机不可能求出系统的解析解 连续 只能求出连续响应曲线上的有限个点 即数值解 连续系统数值积分法 就是利用数值积分方法对常微分方程建立离散化形式的数学模型 差分方程 并求出数值解 最常用的数值解法有 欧拉法 梯形法 Adams Runge Kutta法 描述各类系统最基本的模型用微分方程或状态空间表达式 二次建模就是要求出适合用数字计算机求解的模型 就需要把微分运算转化成算术运算在用计算机求解 第一节数值积分法的基本原理 首先把需仿真研究的系统表示成一阶微分方程组或状态方程的形式 以一阶连续系统为例 微分方程及初值如下 在时间点 处的解为 希望找到一个近似公式来表示方程的近似解 为精确值 的近似值 为准确积分值 的近似值 所谓微分方程初值问题的数值解法就是寻求真解在一系列离散点 上的近似解 数值解 相邻两个时间离散点的间隔 称为计算步距或步长 通常情况为定值 也有变步长 可见 微分方程初值问题数值解法的主要问题归结为对 为此 要先把连续的微分方程用数值积分法转化为离散的差分方程的初值问题 然后根据初始条件X0逐步递推计算出后续时刻 的数值解 如何求出定积分的近似解 数值解法的共同特点是步进式的递推算法 主要有 单步法 多步法和预估 校正法 并有显式和隐式之分 一 欧拉法 EulerMethod Euler法是最简单的一种数值积分法的单步运算 虽然计算精度较差 但几何意义明显 便于理解 能说明构造数值积分算法的基本思想 下面采用三种方法推导出Euler法的数值近似公式 以便对数值积分器的基本思想能透彻了解 第二节数值积分法的单步算法 Taylor级数展开 x t 为解析解 将x t 展开成Taylor级数 以一阶连续系统为例 微分方程及初值如右 只取一次项 其余忽略 写成差分方程为 这就是解微分方程初值问题的欧拉算法 矩形近似 把积分区间h取得足够小 将在 近似为常数 用左矩形面积近似该区间的曲线面积 对方程在上求积分 也能得到 左矩形 也称为前向欧拉法 近似及误差 将在 近似为常数 用右矩形面积近似该区间的曲线面积 得到 这是右矩形欧拉公式 是一个隐式算法 对积分 右矩形 也称为后向欧拉法 近似及误差 切线近似 取切线上处的值来近似 在的一个小邻域内 曲线x t 可用处的切线来表示 x t 在处的斜率为 也能得到 前向欧拉法 过点以fn为斜率的切线方程为 切线近似 过点以fn 1为斜率的切线方程为 取切线上处的值来近似 在曲线x t 可用处的切线来表示 x t 在处的斜率为 也能得到 后向欧拉法 欧拉法 切线推导 的几何意义 欧拉法实际计算时的几何意义 例 设系统方程 试用Euler法求其数值解 取步长h 0 1 解 前向Euler法递推式 有初始条件 可进行递推 后向Euler法递推式 隐式算法 需先解此非线性方程 由此公式可进行递推 前向Euler法与精确值比较 前向Euler法 后向Euler法与精确值比较 前向Euler法在不同步长的结果比较 二 梯形法 Euler法的计算精度较差 如果改用梯形面积代替每个步距的曲线面积 就可提高精度 精确积分应为曲边梯形的面积 现用直边梯形的面积来近似 写成差分方程为 梯形近似及其误差 梯形法实质是采用了 两点斜率平均值的结果 由于利用了两点的信息 从而提高了计算精度 和 这一思想被广泛地应用于许多算法中 实际计算时 如果在每个积分步矩中多取几个点 分别求出其斜率 然后取不同的权值为 后面Runge Kutta法就是采用这样的思想来进行计算的 梯形法的几何意义也可按折线理解 梯形法大大提高了精度 但为隐式算法 每次递推计算时需解一次非线性方程 计算量较大 由此考虑进行改进 先用Euler法计算出 的近似值 代入导函数 求出近似值 再代入梯形公式求解 预估公式 Euler法 校正公式 梯形法 为预估 校正法 也称为改进的Euler法 例 设系统方程 试用梯形法求其数值解 取步长h 0 1 解 梯形法递推式 隐式算法 需先解此非线性方程 由此公式和初始条件可进行递推 见仿真结果 例 设系统方程 用改进欧拉法求数值解 取步长h 0 1 解 改进欧拉法 由此公式和初始条件可进行递推 见仿真结果 前向Euler法 梯形法与精确值比较 前向Euler法 改进欧拉法与精确值比较 梯形法 改进Euler法与精确值比较 三 Runge Kutta法 Taylor级数匹配原理 由于输入函数是t的函数 则将记做得微分方程 如果对变量t x具有各阶导数 可推得x t 的各阶导数 设已知 进行Taylor级数展开 若已知 的值 则当h较小时 可用级数展开的前p 1项作为近似 令 则 即 以上公式 1 就称为p阶的Taylor展开法递推公式 之间的误差为 局部截断误差与hp 1是同阶无穷小量 记为O hp 1 1 欧拉法的Taylor级数展开 只取一次项 其余忽略 写成差分方程为 这就是解微分方程初值问题的欧拉算法 所以欧拉法称为一阶的Taylor展开法递推公式 局部截断误差与h2是同阶无穷小量 记为O h2 梯形法的Taylor级数展开 取一次项和二次项 写成差分方程为 所以梯形法称为二阶的Taylor展开法递推公式 O h3 可见Taylor展开法需用在的高阶导数计算 不便于数值计算 2 Runge Kutta法 于是Runge Kutta法用在一些点上的值表示 使局部截断误差的阶数与Taylor展开法相等 避免了求高阶导数 又保证高的精度 对微分方程在区间 的连续解为 在区间取m个点 若已知 则用它们的一次组合去近似 即 现在的问题是如何求 设已知 因为是未知的 最简单的用欧拉法构造 由Euler法 若按原有欧拉法计算下一步距的值得 以二阶Runge Kutta法为例说明 精度较差 为改进精度 由欧拉法以步长ah得到另一个点 设已知某步 从点开始 沿斜率为方向移动一个步长h 得到 从图中可见三点 可得到比欧拉法精度较好的近似值 在一条直线上 可选取参数b1 b 构造 二阶R K构造法 t X t tm tm 1 xEm 1 x tm 1 tm ah xEm a x tm 欧拉法计算 以步长ah得到另一个点 从点开始 沿斜率为方向移动一个步长h 得到 可选取参数b1 b 构造 以得到比欧拉法精度较好的近似值 以上计算归结为 记 则 如何选取参数a b1 b2 可获得最高的精度 将k2在处展开成Taylor级数 见 代入式子 见 若满足 则与Taylor展开式前三项相同 局部截断误差 三个未知数 两个方程 有多组解 局部截断误差是h的三阶无穷小量 比欧拉法的精度高一阶 若取a 1时 得到 改进的Euler公式 若取时 得 修正的Euler公式 或中点公式 显式p阶Runge Kutta法的一般形式为 3 常用的Runge Kutta法 Kutta三阶法 Heun三阶法 3 经典显式四阶Runge Kutta法 四 微分方程数值积分的矩阵分析法 前述的各类数值积分公式都以一阶系统 单个的微分方程 进行讨论 而实际工程中大量的仿真对象是高阶系统 可用一阶微分方程组来描述 此时 数值积分公式有相应的矩阵形式 矩阵形式的数值积分公式 欧拉法公式 前向欧拉法公式 后向欧拉法公式 梯形法公式 二阶龙格 库塔法公式 改进的欧拉法公式 是预估 校正公式 四阶龙格 库塔法公式 RK4 对于n阶系统 状态向量x为n维 计算中每前进一步h 要计算4n个kij值 对状态空间表达式 此时 RK4公式的 个k值 例 系统方程 取步长h 0 1 试用RK4法求t 0 1 0 2时的解 解 将原系统方程化为状态方程形式 见仿真结果 作业 P1493 2 单步法的特点 计算n 1时刻的值yn 1时 只用到第n时刻的yn和fn 如果能利用多步计算信息 历史时刻值 则可能既加快仿真速度又获得较高的仿真精度 这就是构造多步法的出发点 第三节数值积分法的多步算法 实际在逐步递推过程中 计算yn 1时已经获得一系列的近似值 以及 多步法中以Adams法最具代表性 应用最为普遍 对一阶连续系统 连续解为 现过三点 按插值原理构造一个多项式 来逼近函数 对函数 再对多项式 积分近似 积分 一 Adams算法 得 多项式中的系数由下决定 拉格朗日插值公式 令 同时考虑 因为有 进行变量替换 显然 对多项式的积分计算很容易 微分方程连续解为 写成差分方程 这就是显式两步二阶Adams递推式 显式Adams算法的系数值 显式Adams算法的递推公式为 k 隐式Adams算法的系数值 隐式Adams算法的递推公式为 k 无论用显式或隐式k阶Adams法求解微分方程初值问题数值 需要先知道k 1个初始值 例如三步Adams法 于是初始值 只能从初始条件得到 还需知道 才能求出 需用单步法求出 才能使 多步法的递推计算能够进行 为保证多步法的精度 注意选择相应精度的单步法计算初始值 例 设系统方程 用显式二阶Adams法求解 取步长h 0 1 解 显式二阶Adams法 起步初始值由梯形公式求出 下面就可以用Adams公式进行递推 有初始条件 可进行初值计算 为计算的值 用到时刻以前的值来推导 可获得性能更好的算法 由Adams法得出更为一般的形式 即 再令 二 线性多步法 解出 为计算的值 用到 和相应的导数值 而且公式关于x f是线性的 称为线性k步法 来推导 即 当 称为显式线性k步法 当 称为隐式线性k步法 时的情形 是线性k步法 再看前面的 当 显式Adams算法的递推公式为 隐式Adams算法的递推公式为 接下来的问题就是线性k步法 如何确定其中的常数 应用Taylor级数匹配原理 使局部截断误差尽可能的小 将 处展开成Taylor级数 和 再带入误差公式得 其中 对于 若能选取 使 则此线性多步法为p阶k步算法 3 构造线性多步法 确定步数k 由Taylor级数匹配原理的 得到方程求出待定系数 得出尽可能高阶的算法 以二步法为例说明构造步骤 5个未知数 4个方程 得 得到一般的线性二步法形式 是三阶二步法 又可得 是四阶二步法 称为Milne算法 隐式算法 为三阶隐式算法 就是三阶隐式Adams法 如何得到二步显式算法 得到二步显式算法 但这是一个数值不稳定算法 计算中出现的微小误差会迅速增长 因此 二步显式算法无法达到三阶 最多达到二阶 只用三个方程 二阶二步显式算法的一般形式 二阶二步Adams显式算法 为Adams法 为向后微分公式 显然是隐式公式 若要求向后微分公式为K阶算法 则 当K 1 为向后Eular法 K 1 2 6的向后微分公式系数值 k 计算时 仅已知 和相应的导数值 称为线性多步法的起步值 需用单步法来求出 若线性多步法是P阶算法 计算起步值的算法不应低于P阶 否则影响计算精度 仿真步长h的选取 是否影响仿真结果 先看一个例子 第四节数值积分法稳定性分析 系统稳定 精确解 用前向欧拉法求解 ch3 6 m 由此可见 仿真步长h的选取 会影响仿真结果 用前向欧拉法求解 当h 0 2时不稳定 是由于步长太大 从而截断误差太大造成的 一 数值解法稳定性的含义 数值解的稳定性 在扰动 初始误差 舍入误差 截断误差 的影响下 计算过程中的累积误差不会随计算步数的增加而无限增长 微分方程的数值积分方法 实质是微分方程的差分化 然后从初始条件递推迭代 不同数值解法对应着不同的差分方程 是否稳定取决于该差分方程的特征根是否满足稳定性要求 处于Z平面上以原点为圆心的单位圆内 着重研究单步法的稳定性对步长的限制 二 数值解法稳定性分析 这样做的根据是 1 试验模型简单 对其数值不稳定的方法 不可用 2 一般的初始问题在其解的存在区域内 可局部线性化转化为试验方程 只要原系统是稳定的 即 不等式成立 称后向欧拉法是恒稳定的算法 只要原系统是稳定的 即 不等式成立 与后向欧拉法一样