摘要: 求解微分方程是连续系统半实物仿真中的一项关键技术 应用于复杂系统半实物仿真的实时算法要求具有较大的稳定区域较高的计算精度和较小的运算量，以便满足实时性要求 常用的 法基于离散相似原理，计算量小但稳定区域小 常用的实时 法基于泰勒级数匹配原理，稳定区域较大但计算量大 为了满足复杂系统半实物仿真的需要，综合了泰勒级数匹配原理和离散相似原理的优点，构造了一类实时的四阶混合数字仿真算法，并且分析了该类算法的稳定性实时性和计算精度 最后，利用标准算例对该类算法进行了检验，仿真结果表明，该类算法计算量小稳定区域大精度高，在半实物仿真中具有一定的优越性

关键词: 数字仿真算法; 混合法; 精度; 稳定性; 实时

引言

连续系统的数学模型一般采用微分方程表示，根 据 统计，在一帧的仿真计算中，求解微分方程占据了大部分的计算量 因而，在连续系统的数字仿真中，求解微分方程具有重要的地位 实时仿真算法是半实物仿真的基础 在半实物仿真中，由于有实物进入仿真回路，故而，对仿真算法的实时性稳定性和精度有较高的要求 常用的 法基于离散相似原理，计算量小但稳定区域小 常用的实时 法基于泰勒级数匹配原理，稳定区域较大但计算量大 为了满足复杂系统半实物仿真的需要，本文综合了泰勒级数匹配原理和离散相似原理的优点，构造了一类实时的四阶混合数字仿真算法，综上所述，构造优良的仿真算法求解微分方程，在满足精度的前提下，使得仿真算法具有较少的计算量和较大的稳定区域是追求的目标

本文在研究了基于泰勒级数匹配原理的仿真算法和基于离散相似原理的仿真算法的基础上，构造了一类实时的四阶混合实时算 法，并 对 其 精 度稳 定 区 域计 算 量 进 行 了 分析，最后，通过算例对算法进行了检验

构造混合算法

构造混合算法的思想对于一般的非线性微分方程:( ( ) ， ( ) ) ，( ) ( )其中( ) 为状态向量，( ) 为外部输入向量计算机仿真输入为 ( ) ，它来自实物系统或物理环境

泰勒级数匹配原理是先计算出点( ，( ) ) 邻域中的若干个点( ，) 上的值( ，) ，利用这些值构造( ) 的近似值 的计算公式，并且要求该公式在( ，( ) ) 处的泰勒展开式与真值( ) 在 处的泰勒展开式的前若干项一致例如: 算法其算法的一般形式为:， ( )， ( )，，，离散相似原理是将连续系统通过虚拟采样开关和信号保持器转换成离散系统，由此建立数字仿真模型例 如:算法其算法的一般形式为:[ ] ( )

算法的构造思想是在一个时间步长的计算中，安排多个子步长计算相应的右函数值，最后，利用这些右函数值和当前帧的被积变量值计算出下一帧的被积变量值而 算法的思想是利用一系列各时间点右函数的历史值和当前帧的被积变量值计算出下一帧的被积变量值 据此，本文将两者结合起来，采用混合法

构造混合算法的过程

首先在一个时间步长的计算中，利用一系列各时间点右函数的历史值预估出一个中间时刻的被积变量值，并计算出相应时刻的右函数值然后，再利用这个右函数值一系列各时间点右函数的历史值和当前帧的被积变量值计算出下一帧的被积变量值构造的算法如式( ) 所示， ( ){( )其中对于式( ) 中的第一个公式( 以下简称预估公式) ，进行截断误差的分析注意 到 的 精 确 值 为 ( )， 的 精 确 值 为( ) 将预估公式进行泰勒展开，则有( )( )其中:( ) ， ( )( )为了保证预估公式的精度达到 阶，要求( )的多项式中，零阶项一阶导数项二阶导数项三阶导数项的系数均为零，则有( )对于式( ) 中的第三个公式( 下面简称校正公式) ，进行截断误差的分析注意到 的精确值为( ) 将校正公式进行泰勒展开，则有( )( )( ) ( )( )( )其中:( ) ， ( ) ( )为了保证校正公式的精度达到 阶，要求( )的多项式中，零阶项一阶导数项二阶导数项三阶导数项四阶导数项的系数均为零，则有( )将作为参数，可以求出相应的各个系数的值本文选取了 ，，三组值，构造了 三种混合算法各系数如表 所示表 四阶混合算法的系数算法四阶混合法的稳定性实时性精度分析稳定性分析下面对本文所设计的 等三种算法的稳定区域进行分析，并与常用实时算法 () 的稳定区域进行比较首先分析四阶混合法的稳定区域将线性试验方程 代入式( ) ，记 ，则有( )( )( )( )( )( )因而可得四阶混合法的稳定多项式( )( )取 ，可得 的稳定多项式( )的稳定区域如图 所示图 算法的稳定区域取图 中曲线所包围的区域中的一点( ，) ，即，代入式( ) ，则有( )利用复函数方程根的下山法［］，用 编程，求解出式( ) 的三个复数根为{( )显然， ， ，因此，图中曲线及 直线所包围的内部及边界区域为稳定区域在实轴上的稳定区域为( ，)取 ，可得 的稳定多项式( )同理，可以确 定 的稳定区域为图 中 曲 线 及直线所包围的内部及边界区域图 算法的稳定区域在实轴上的稳定区域为( ，)取 ，可得 的稳定多项式( )同理，可以确 定 的稳定区域为图 中 曲 线 及直线所包围的内部及边界区域图 算法的稳定区域在实轴上的稳定区域为( ，)将线性试验方程 代入 的 算 法 公 式，可 得的稳定多项式为( )可以确定 的稳定区域为图 中曲线及 直线所包围的内部及边界区域图 算法的稳定区域在实轴上的稳定区域为( ，)将线性试验方程 代入 的 算 法 公 式，可 得的稳定多项式为( )可以确定 的稳定区域为图 中曲线及 直线所包围的内部及边界区域图 算法的稳定区域在实轴上的稳定区域为( ，)将线性试验方程 代入 的算法公式，可得的稳定多项式为( )可以确定 的稳定区域为图 中曲线及 直线所包围的内部及边界区域在实轴上的稳定区域为( ，)以为可变参数，以 混合算法的稳定区域在负实轴上交点 的绝对值最大为目标函数，编制程序进行寻优结果表明: 当从逐步增大到，的绝对值逐步增大，到时其最大值为 但是，考虑到 算法首先需要对进行预估，然后，再计算 的值为了保证在当前帧完成所有的计算，因此，不能取得太大，以便给校正公式的计图 算法的稳定区域算留出充足的时间为此，本文只选择了 ，，三种情况进行研究本文中的混合法及常用的实时算法的数值稳定性区域( 实轴上部分) 如表 表 表 所示表 算法的数值稳定性区域( 实轴上部分)算法 数值稳定性区域表 算法的数值稳定性区域( 实轴上部分)算法 数值稳定性区域表 算法的数值稳定性区域( 实轴上部分)算法 数值稳定性区域从表 中可以看出，本文所设计的 算法，其稳定区域明显高于 法，而与 法接近但是， 的算法精度明显高于 ，适用于高精度仿真的情况实时性分析算法的实时性要求算法中所用到的输入信号都应该是数字处理过程中已经从实物系统或其它过程获取的下面分析式( ) 中的三个公式第一个公式中，是前一帧积分的结果，， 是前一帧前两帧的右函数值均为已知，而是在当前帧可以计算的因此，满足实时性要求第二个公式中， 已经由预 估 公 式 计 算 得 到， 可以从外界采样得到，因此，满足实时性要求第三个公式中，，，，，均已计算得到，因此，满足实时性要求值得一提的是，在第一帧计算时，， 是没有值的，在第二帧计算时， 是没有值的 因此，本算法是不能自起步的 本文在前两帧的计算中，采用 算法进行起步计算，从第三帧开始采用式( ) 算法进行计算 同样，也是不能自起步的，也采用 算法进行起步计算从第三帧开始采用 算法进行计算 从第四帧开始采用 算法进行计算精度及计算量分析均达到 阶精度，同样 也达到了 阶 精 度，而 只 达 到 阶 精 度 显 然，具有精度优势在保证阶精度的前提下， 均只需要计算两次右函数， 只需要计算一次右函数，而只 达 到 阶精度却需要计算 次 右 函 数 显 然，的计算量 最 少， 的 计 算 量次之， 的计算量最大综上所述， 虽然计算量最少，但是其稳定区域比其它两组算法小很多; 而 的稳定区域最大，但是，其精度只达到 阶，且计算量较大;的稳定区域接近 ，计算量不大，而精度较高，因此，在半实物仿真中，具有一定的综合优势混合法仿真实例将本文所构造的实时算法 用下面的衰减振荡型微分方程组的算例进行检验，( )，( )( )该方程组的精确解为( ) ( )( ) ( )( )检验过程采用双精度计算，步长取为 ，积分区间为［，］采用 编程 表 为求解方程组后的数值解与精确解的误差 ( )表 为 求解方程组后的数值解与精确解的误 差 ( ) 表 为改变积分步长的值时，求解方程组后的数值解与精确解的误差 ( )表 混合法的精度比较表 常用实时算法的精度比较表 步长变化时各算法的误差结果从表 表 表 可以看出，本文所构造的四阶混合法的精度高于常用的实时算法

结论

本文综合了泰勒级数匹配原理和离散相似原理的优点，构造了一类实时的四阶混合数字仿真算法，并且分析了该类算法的稳定性实时性和计算精度 最后，利用标准算例对该类算法进行了检验，仿真结果表明，该类算法计算量小稳定区域大精度高，在半实物仿真中具有一定的优越性

    参考文献:［］熊光楞数字仿真算法与软件［］北京: 宇航出版社，［］黄柯棣，等系统仿真技术［］长沙: 国防科技大学出版社，［］陈进，等科学与工程计算常用算法程序库［］上海: 上海交通大学出版社，［］尹彦芝 语言常用算法与子程序［］北京: 清华大学出版社，［］黄振全，等低阶实时最优 算法［］数值计算与计算机应用， ，( ) :