《火电厂仿真课件--第5章 数值积分法仿真》由会员分享，可在线阅读，更多相关《火电厂仿真课件--第5章 数值积分法仿真（39页珍藏版）》请在金锄头文库上搜索。
1、第第 5 5 章章 面向微分方程的数值积分法仿真面向微分方程的数值积分法仿真 数值积分是数值分析的一个基本问题 也是复杂计算问题中的一个基本组成部分 数值积分往往用极简单的方法就能较好地得 出对所求解的具体数值问题的解答 但数值积分的难点在于计算时间有时会过长 有时会出现数值不稳定现象 另外 数值积分的理论性较强 其理论和方 法都已经比较成熟 计算精度也比较高 5 1 仿真中研究数值积分法的意义 数值解的一种近似方法 对于连续系统的高阶微分方程 可化为若干个一阶微分方程组成的方程组 数值积分法是求解微分方程 例如 下式所示的状态方程 可以化为一个一阶微分方程组 5 15 1 5 25 2 所以
2、 连续系统的仿真就是从给定的初始条件出发 对描述系统动态特性的常微分方程或常微分方程组进行求 解 从而得到系统在一定输入作用下的变化过程 在求解这些微分方程时 最常用 也是最有效的一种 方法就是数值积分法 5 35 3 5 2 数值积分法仿真的基本原理 对微分方程 5 1 两端同时取积分 可得 当 时 上式变为 5 45 4 5 55 5 将积分项拆成两项 5 65 6 故上式可写为 5 75 7 此式是方程 5 1 在tn 1时刻的精确解 在数值解法中 希望用近似解 代替准确解 其中 5 85 8 为 为 为 的近似值 令 称为计算步长或步距 式 式 5 85 8 是从常微分方程 5 1 出
3、发建立的离散数学 模型 差分方程 由此可见 数值积分法就是在已知微分方程初值的情况下 求解该方程在一系列离散点 处的近似值 其 特点是步进式 根据初始值逐步递推地计算出以后各时 刻的值 从式 5 8 可知 数值积分法的主要问题归结为如何对f t y 进行数值积分 求出f t y 在区间 tn tn 1 上 定积分的近似值Qn 采用不同的方法求Qn 就出现了各种各样的数值积分方法 不同的数值积分将对求解的精度 速度和数值稳定性会 产生不同的影响 这将在下述内容中具体介绍 数值积分法种类繁多 在此从实用角度介绍几种基本的方法 5 3 欧拉 Euler 法 5 3 1 简单欧拉法 欧拉法是一种最简单
4、的数值积分法 对于方程 在区间 tn tn 1 上求积分 得到 若区间 tn tn 1 足够小 则 tn tn 1 上的f t y 可近似地 看成常数f tn yn 故可用矩形面积近似代替 即 tntn 1 f tn yn 于是有 5 95 9 将此式写成差分方程为 5 105 10 著名的欧拉公式 5 3 2 改进的欧拉法 如果用梯形面积而不是矩形面积来代替每一个小区间上的 曲线积分 就可以提高计算精度 梯形法的计算公式为 5 11 式中的右端含有待求量yn 1 因而它是隐函数形式 这种 方法不能自行启动运算 需要依赖其它算法的帮助 每次计算都用欧拉法算出y t n 1 的近似值 以 此计算
5、近似值 然后利用梯形公式 5 11 求出修正后的 即有 帮助方法 预估式 校正式 5 12 改进的欧拉公式 5 3 3 几个基本概念 简单的欧拉法是用前一时刻tn的yn求出后一时刻的yn 1 这种方法称为单步法 它是一种自行启动的算法 如果 求yn 1时需要tn tn 1 tn 2 时刻yn yn 1 yn 2 的值 这 种方法为多步法 改进的欧拉法为两步法 它是一种不 能自行启动的算法 1 单步法与多步法 简单的欧拉法表达式的右端 计算 所用的数据均已 求出 这种公式称为显式公式 改进的欧拉法表达式的右端 有待求量 这种公式 称为隐式公式 隐式公式不能自行启动 需要用预估 校正法 单步法和显
6、式在计算上比多步法和隐式方便 但有时 为了满足精度 稳定性等方面的要求 需要采用隐式算 法 2 显式与隐式 3 截断误差 这里用泰勒级数为工具来分析数值积分公式的精度 假定yn是精确的 用泰勒级数表示 处的精确解 即 显然 简单的欧拉法是从以上精确解中取前两项之和来 近似计算 每一步由这种方法引入的误差称为局部截断 误差 简称截断误差 简单的欧拉法的截断误差为 不同的数值方法有不同的截断误差 一般若截断误 差为 则方法为r阶的 所以方法的阶数可以作为 衡量方法精确度的一个重要标志 5 13 由此可见 简单的欧拉法是1阶精度 改进的欧拉法 由于采用了平均斜率 相当于取了泰勒级数的前3项 因 此为
7、2阶精度 分析欧拉法截断误差的思想 同样也适用于其它数值 积分方法 4 舍入误差 由于计算机的字长有限 数字不能表示得完全精确 在计算过程中 不可避免地会产生舍入误差 舍入误差与计算步长成反比 如果计算步长小 计 算次数多 则舍入误差就大 舍入误差除了与计算机字长有关以外 还与计算机 所使用的数字系统 数的运算次序以及计算所用的子程 序的精度等因素有关 5 数值解的稳定性问题 采用数值积分法求解稳定的常微分方程 应该保持 原系统的稳定特征 但是 1 在计算机逐步计算时 初始数据的误差及计算 过程的舍入误差对后面的计算结果将产生影响 2 如果计算步长取的不合格 有可能使仿真出现 不稳定的结果 下
8、面我们简单讨论一下这个问题 差分方程的解与微分方程的解类似 可分为特解和通 解两部分 与稳定性有关的是方程的通解 它取决于差 分方程的特征根是否满足稳定性条件 例如 为了考查欧拉法的稳定性 我们研究检验方程 Test Equation 其中 为方程的特征根 对此有 显然 要使该差分方程是稳定的 必须使下式成立 即 表明 为使数字仿真稳定 对计算步长应有所限制 另外 稳定性还与系统的特性以及数值积分方法有关 上述分析欧拉法稳定性的思想 同样适用于其它数值积 分方法 5 14 5 15 5 4 龙格 库塔 Runge Kutta 法 由前面的分析可知 将泰勒展开式多取几项以后截 断 就能提高截断误
9、差的阶数和计算精度 然而 直接 采用泰勒展开方法要计算函数的高阶导数 运用起来不 便 龙格 库塔方法的基本思想是 用几个点上的函数值 的线性组合代替函数的各阶导数 然后按泰勒级数展开 确定其中的系数 这样既可避免计算高阶导数 又可提 高积分的精度 龙格 库塔法有多种形式 以下从实用角度直接给出 公式的形式和相应的精度 5 4 1 龙格 库塔方法的基本思想 5 4 2 二阶龙格 库塔方法 2阶龙格 库塔方法的公式为 5 16 上式表示的数值解是用的泰勒级数在2阶导数以后截断 所求得的 因此称为2阶方法 故 2阶龙格 库塔法与式改进的欧拉法相比 实质完全相同 所以改进的欧拉法实质上是2阶龙格 库塔
10、法 截断误差为 实时仿真的2阶龙格 库塔方法 5 17 其截断误差为 5 4 3 四阶龙格 库塔方法 4阶龙格 库塔法是一种最常用的方法 其经典表达式为 5 18 其截断误差为 显然 4阶龙格 库塔法的计算量较大 但计算精度 较高 在比较不同算法的计算精度时 常以它的计算结 果作为标准 实时仿真的4阶龙格 库塔方法 5 19 以上公式都是标量形式 如果要换成向量形式 只 要把式中的标量y f k换成向量Y F K即可 从理论上讲 可以构造任意阶数的龙格 库塔方法 但是 精度的阶数与计算函数值f 的次数之间并非等量增 加的关系 见下表所列 由此可见 4阶经典龙格 库塔方法有其优越性 而4 阶以上
11、的龙格 库塔方法计算f所需的次数比阶数多 增加 了计算量 从而限制了更高龙格 库塔方法的应用 对于 大量的实际问题 4阶方法已可满足精度要求 所以得到 了广泛的应用 我们仍采用检验方程 进行讨论 对它利用泰勒级 数展开得 5 4 4 龙格 库塔公式的稳定区域 对于 有 将它代入上式得 5 20 5 21 令 将其代入上式 得到该式的稳定条件为 5 22 由此稳定条件 下表给出了各阶龙格 库塔公式的稳定区域 表5 2 龙格 库塔公式的稳定区域 在使用龙格 库塔公式时 选取的步长应使落在稳定 区域内 否则 在计算时会产生很大的误差 从而得不到 稳定的数值解 例如 用4阶龙格 库塔方法解 1 用解析
12、法求解 本例是稳定的 2 用数值法求解 当h 0 1 时 数值解是稳定的 当h 0 2时 数值解就不稳定了 这时因为 此数值在稳定区间以外 所以数值解不 能收敛 这种对步长有限制的数值积分法称为条件稳定积分法 从4阶龙格 库塔方法的稳定条件 中可以看出 系统的特征根越大 需要的积分步长 就越小 这一点可作为选择步长的依据 步长的大小 除了与数值积分方法的阶数有关外 还与方程本身的性质有关 除以上介绍的欧拉法 龙格 库塔法外 数值积分 方法还有许多种 如亚当姆斯方法 吉尔方法等等 此 处不一一介绍 5 5 计算方法的选择 数值积分方法很多 在实际使用时存在一个选择问 题 对于一个具体问题 如何选
13、择具有一定的难度 至 今尚无一种确定的方法 一般来说 数值积分方法的选择应考虑的因素有 1 精度要求 数值积分法的精度受 截断误差 舍入误差 积累误差 的影响 而这些 误差与积分方法 步长 计算时间 计算机精度等有关 一般 1 积分方法的阶数越高 截断误差越小 精度越高 2 步长越小 精度越高 3 多步法的精度高于单步法 4 隐式算法的精度高于显式算法 因此 当需要高精度时 可采用高阶 多步 较小 步长 隐式算法 但是 步长的减小往往会增加迭代次数并增大舍入 误差和积累误差 与此相反 在精度要求不高时 最好 使用低阶方法 计算速度取决于计算步数 每步积分所需的时间 而每步的计算时间又与积分方法
14、有关 它主要取决于 计算导数的次数 4阶龙格 库塔方法每步要计算4次导数 导数的计算是最费时的 为了加快计算速度 在积分方法已定的条件下 应要 保证精度的前提下 尽量选择较大步长 以缩短积分时 间 2 计算速度 数值解的稳定性必须保证 否则 计算结果将失 去真实意义 从稳定性来看 不同的数值积分算法有不同的稳 定性 应用时应控制步长h 使数值积分算法在稳定域 内 3 数值稳定性 步长的选择很重要 步长过大会增大截断误差 甚 至出现数值不稳定现象 过小了又因增加了步数 而使 舍入误差增大 所以 仿真的总误差与步长的关系不是 单调函数关系 而是一个具有极值的函数 如图所示 5 6 计算步长的选择
15、误差 总误差 截断误差 舍入误差 步长 由图可知 存在一个最佳步长 当积分方法确定以后 在选取步长时 需要考虑的 一个重要因素就是仿真系统的动态响应特性 如果系统 的动态响应快 导数变化激烈 则应取高阶方法和小步 长进行计算 为了保证数值稳定性 步长应限制在最小时间常数 相当于最大特征值的倒数 的数量级上 用经验方法选取步长的两种方法 1 由系统方程中最小时间常数Tmin来决定 2 由系统开环频率特性的剪切频率 来决定 5 7 面向微分方程的仿真程序构成 一般的仿真程序的组成可用下图表示 输入或设置输入或设置 参数子程序参数子程序 主程序主程序 运行管理运行管理 子程序子程序 解答存储及解答存
16、储及 打印子程序打印子程序 在线绘图在线绘图 子程序子程序 仿真系统计算仿真系统计算 子程序子程序 图中各方块的功能为 主程序 实现对整个仿真计算的逻辑控制 输入或预置参数模块 输入系统的参数初值 计算步长 计算时间等参数 运行管理模块 这是数字仿真程序的核心 对仿真计算 进行时间控制 以保证计算机按要求进行计算及输出 计算模块 根据被仿真的系统及所选的仿真方法编写的计 算程序 输出及显示模块 将仿真结果以数据或图表的形式输出 给用户 4阶龙格 库塔法仿真程序设计 目前 应用于系统仿真的商用软件包诸多 但无论 是开发者还是应用者 了解程序的设计思想 分析程序 的结构与功能 完善和编写仿真程序都是必要的 因此 此处通过对经典龙格 库塔方法 即式 7 22 仿真 程序的介绍 使读者了解和掌握编写仿真程序的基本技 术 为了便于说明程序的实质 我们用具体的例子来讨 论 假设有一个2阶系统 令 则 以下是用C语言编写的原理性程序 程序中变量说明 i n 分别为循环控制和积分器个数变量 h t tmax 分别为步长 时间和仿真总时间变量 yn 10 存放tn时刻y值的临时变量数组 y 10 存放y
