一种仿真中刚性常微分方程初值问题的数值积分求解方法与流程

1.本发明属于仿真技术领域，涉及一种仿真中刚性常微分方程初值问题的数值积分求解方法。


背景技术：

2.工程技术中的很多现象，如控制系统的运行、电力系统的运行、飞行器的运动、化学反应的过程等，其数学模型都是常微分方程(ode,ordinary differential equation)的初值问题(ivp,initialvalue problem)，很多偏微分方程也可以近似简化成为常微分方程的初值问题进行求解。一般来说，常微分方程的数值解法是比较成熟的，理论也很完善，有许多成熟的求解算法可供选择。
3.然而，有一类常微分方程在数值求解时会遇到相当大的困难，这类常微分方程解的分量有的变化很快，有的变化很慢。从数值求解的观点来看，解的变化快时应采用较小的积分步长，当变化快的分量趋于稳定时，就可以采用较大的积分步长。但是大量的理论和实践都证明，无论解的变化是否趋稳，积分步长都不能放大，否则就会出现数值求解不稳定进而发散。当进行系统动力学分析研究时，这种现象会给数值求解带来很大的难度和困扰，尤其是当系统规模较大或涉及多个学科时这种现象非常普遍。
4.常微分方程的这种性质叫做刚性(stiff)，这一类问题叫做刚性问题。刚性问题在有的文献中也被叫做病态方程或者坏条件方程，这种现象在实际工程中经常会遇到，如控制系统中被控对象惯性较大，具有较大的时间尺度，而控制器则反应灵敏，时间尺度较小；对于飞行动力学来说，其姿态运动较快而质心运动相对较慢；在热流体系统中，压力的波动很快，在流体中以音速传播，而传热过程较慢，其变化尺度显然不在同一量级上；对于复杂的电力系统，受寄生电容的影响，不同元件之间的时间尺度也会有非常大的差异。
5.伴随着20世纪五十年代开始的以原子能、航空航天、计算机等高新技术为代表的第三次工业革命浪潮，系统动力学问题数值计算也逐渐从数学理论走向工程实践，进而诞生了众多开源程序和商业仿真工具。
6.应用比较广泛的开源求解器包括：采用bdf算法且至今依然经常更新版本的cvode；采用隐式runge-kutta法的radau5；采用显式runge-kutta法的dop853和dopri5；采用对角隐式runge-kutta法的sdirk；基于改进广义bdf算法的mebdf；以经典bdf算法为基础设计并被广泛使用的difsub和lsode；python的开源计算包scipy.integrate.ode(其底层仍为调用lsode和odepack)；r语言中基于显式runge-kutta法的desolve工程包(根据官网描述，其底层仍为调用fortran语言的开源odepack计算包、基于显式euler/runge-kutta法的程序包以及基于隐式runge-kutta法的程序包radau)；用于微分方程高性能求解器的julia套件“differentialequations.jl”等等。
7.除了开源求解器，还有很多非常优秀的商业仿真求解工具也在上世纪七十年代后逐渐出现在大众眼前。在所有商业系统仿真工具中，普遍被认为综合性能最好的就是matlab/simulink，其ode求解器包含了几乎所有能见到的常规算法，如bdf法、改进bdf法、
adams法、显式或隐式runge-kutta法等；此外，dymola、simulationx、motion、adams和amesim也是很好的动力学求解工具，由于这些工具处理的模型主要是多学科或多体动力学等问题，这一类模型普遍具有较大刚性，所以这些软件的求解器基本都采用隐式多步法，最常见的就是bdf及其改进方法。
8.我国的mworks也是一款用于modelica模型动力学求解的仿真工具，但其求解器仍是通过调用一些国外的算法包来实现的。与国外现有技术积累相比，我国在以下两方面存在明显不足：一方面，国内没有国产自主的ode-ivp求解器，大都调用国外程序包(含商业和开源)，如cvode、dassl、radau5、mebdf、lsode等，国内在这一领域完全空白；另一方面，国外的开源程序包形成年代较早，大都没有包含近年来新的方法。此外，国外的商业工具进入中国市场普遍在2000年前后，求解器内核并没有大的改动，同开源程序包一样，没有将近些年新的算法纳入其中。


技术实现要素：

9.本发明解决的技术问题在于提供一种仿真中刚性常微分方程初值问题的数值积分求解方法，为系统级仿真(尤其是大规模、多学科系统仿真)求解提供底层支持。
10.本发明是通过以下技术方案来实现：
11.一种仿真中刚性常微分方程初值问题的数值积分求解方法，包括以下操作：
12.1)针对待求解的问题，选取包括积分求解算法、初始积分步长、收敛残差在内的必要参数；
13.其中积分求解算法从bdf算法、ibdf算法、sdisrk34算法和sdisrk44算法中选择；
14.2)根据所选取的求解算法进行判断：如采用多步法，则再判断当前求解器采用的是bdf或ibdf法，判断完成后继续执行3)；
15.如采用单步法，则再判断当前求解器采用的是sdisrk34算法还是sdisrk44算法，判断完成后跳转执行5)；
16.3)根据选取的求解算法，构建目标非线性方程组；
17.4)利用多种预测步算法分别执行得到多个真实解的预测值，将其分别代入目标非线性方程组，从中选取误差最小的一组解向量作为预测步的解向量；
18.5)将预测步得到的解向量作为迭代初值，求解目标非线性方程组的最小二乘解；
19.在求解时先构建迭代公式，并引入迭代步长的阻尼系数；然后更新迭代步长，判断是否满足收敛要求；如满足则输出结果作为真实解向量，否则继续重复迭代过程不断得到新的解向量，直至满足收敛误差要求；
20.6)执行变步长积分策略：分别执行提供积分步长的步长寻优机制、判断改变积分步长的切换机制、保证计算结果正确的变步长数值回插机制；
21.变步长积分策略执行完成后，判断误差是否满足要求；如满足则执行7)；
22.如不满足则返回1)，并按照变步长策略重新设置积分步长；
23.7)时间轴向前推进一个时间步，判断是否需要结束仿真；
24.如满足则执行8)；如不满足则返回1)，按照变步长策略重新设置积分步长；
25.8)结束仿真计算，存储并显示计算结果。
26.与现有技术相比，本发明具有以下有益的技术效果：
27.本发明提供的仿真中刚性常微分方程初值问题的数值积分求解方法，对标国外主流商业和开源系统仿真求解器，提供了适用于求解刚性和非刚性常微分方程初值问题(ode-ivp)的求解器，为系统级仿真(尤其是大规模、多学科系统仿真)求解提供底层支持。
28.本发明在经典levenberg-marquardt优化算法的基础上融合阻尼牛顿法，提出了一种damped-levenberg-marquardt优化方法(简称d-lm方法)；这种算法能够兼顾运行效率与稳定性，通过与国外商业工具matlab中的同类计算程序相比，本方法在效率上与国外商业工具相当。
29.本发明在预测-校正法中通过多种相对简单的显式算法进行解的预测，目的是给出一个更加靠近真实解的点，保证下一步的隐式求解更容易收敛，这一步是出于保护求解稳定性考虑的，但会略微牺牲一点点效率；这样的处理方法尽管会在预测步环节损失一些计算效率，但能够为校正步提供一个更好的迭代初值，进而减少迭代次数m，通过综合权衡与大量测试后，发现这样的处理对全局效率是有明显提升的。
30.本发明还采用了可大幅提升积分求解效率的变步长新策略，在填补国内空白的基础上，在运算效率上与国外同类产品处于同一水平；进一步的在保证求解精度基础上提出了d-lm方法，提高了运算稳定性。与国外商用系统仿真软件工具相比，本发明的数值方法更新颖，对1990年代之后的新的数值求解方法进行了实现与应用，数值求解过程更加稳定。
附图说明
31.图1为本发明的整体流程示意图，
32.图2-图5分别为图1中第一部分-第四部分的放大图；
33.图6为3级4阶单对角隐式辛rk方法的butcher矩阵示意图；
34.图7为4级4阶单对角隐式辛rk方法的butcher矩阵示意图；
35.图8为变步长积分的工作流程示意图；
36.图9为求解器核心部分代码的组织架构示意图；
37.图10为求解器中d-lm算法的计算结果示意；
38.图11为利用matlab中fsolve命令的计算结果；
39.图12为d-lm方法和fsolve命令的计算耗时比较；
40.图13为利用simulink商业软件求解刚性方程的结果；
41.图14为本求解器和simulink商业软件计算结果的比较；
42.图15为本求解器计算效率与国外商业工具的比较；
43.图16为y1和y2的数值解及其和解析解的误差；
44.图17为y3和y4的数值解及其和解析解的误差；
45.图18为lorenz系统三维演化轨迹示意图；
46.图19为管路模型示意图；
47.图20为管路模型的管路左右两端的压力示意图；
48.图21为管路模型和压力边界条件下本发明与simulink计算比较示意图。
具体实施方式
49.下面结合实施例对本发明做进一步详细描述，所述是对本发明的解释而不是限
定。
50.要实现仿真中高效、稳定的求解器，其主要关注三个方面：
51.首先，要选择合适的计算方法，既要具有较好的稳定性，否则不易收敛；同时还要兼顾数值求解精度，选取的数值算法阶数不能过低；
52.其次，还要充分考虑选取的算法是否适用于典型的刚性问题，包括但不限于：刚性振荡问题、非线性延迟微分问题、带有离散事件的刚性问题等，这些典型问题在动力系统仿真中广泛存在；
53.最后，实现算法的过程中，要尽可能提高代码的运行效率，否则没有工程意义。
54.为此，本发明提出以下方法和策略来进行求解。
55.1、基于线性多步法和runge-kutta法开发ode求解器
56.在线性多步法和runge-kutta法的理论基础上，本求解器选取了bdf、改进bdf以及两种4阶的单对角隐式辛rk这四种稳定性较好且精度也能满足要求的数值算法的作为底层核心处理方法。所提出的适用于刚性常微分方程初值问题的数值积分求解器，可以高效、稳定地求解ode-ivp问题，尤其是针对刚性问题具有较好的效果。本求解器对上述四种求解方法均提供变步长机制，可充分发挥其大稳定域的优势。
57.2、求解非线性方程组最小二乘问题的新方法
58.在应用隐格式算法求解常微分方程初值问题时，需要计算非线性方程组最小二乘问题。本求解器在经典levenberg-marquardt优化算法的基础上融合阻尼牛顿法，提出了一种damped-levenberg-marquardt优化方法，简称d-lm法。这种算法能够兼顾运行效率与稳定性，通过与国外商业工具matlab中的同类计算程序相比，本方法在效率上与国外商业工具相当。
59.3、预测步算法的优选策略
60.线性多步法是一种隐格式求解方法，求解隐式方程最常用的方法是“预测-校正”法，这一类方法也被称作pc方法类。“预测-校正”法的一般步骤是p(ec)m，其中p表示预测(predictor)，e表示预测步结果求得的导数值，c表示校正(corrector)，m表示多次迭代。
61.在本求解器中，为了保证m次隐式迭代的收敛性和求解效率，需要在预测步中尽可能给出一个靠近真实解的预测结果，本求解器提供了5种计算方法，分别应用这5种计算方法得到5组预测结果，在其中优选最靠近真实解的一组作为校正步的输入。这样的处理方法尽管会在预测步环节损失一些计算效率，但能够为校正步提供一个更好的迭代初值，进而减少迭代次数m，通过综合权衡与大量测试后，发现这样的处理对全局效率是有明显提升的。
62.4、可大幅提升积分求解效率的变步长新策略
63.传统的变步长求解策略一般先尝试以h为积分步长进行求解得到下一步的数值y
(h)
，然后按照给定比例缩小步长(如h/2或h/3)再次求解得到下一步的数值y
(h/2)
或y
(h/3)
。随后，通过判断||y
(h)-y
(h/2)
||或||y
(h)-y
(h/3)
||的某种范数是否小于预先给定的误差ε，如小于则认为当前步长合适，时间轴向前推进；如不小于则认为当前时间步不合适，需要按照某种步长缩小策略，重新计算新的积分步长h
new
。
64.本发明在上述策略的基础上额外增加一个判断策略，在变步长过程中引入“步长调节”和“步长稳定”概念，当求解过程中连续若干个积分步长h和误差ε都保持在一个提前
约定的范围内，那么就认为当前被求解模型进入“步长稳定”状态，在该状态下无需进行上述积分步长的求解与调整，直接使用前一步积分步长即可；反之，则认为当前被求解模型进入“步长调节”状态，则按照传统的变步长求解策略进行求解。这种处理方法的最大好处就是在保证计算结果准确性的前提下，大幅减少了计算次数，进而提升求解效率。
65.下面对本发明的仿真求解方法具体说明。
66.本发明提供的仿真中刚性常微分方程初值问题的数值积分求解方法，目的可以等效为求解一个实数域的n维常微分方程的初值问题，其表达式为：
[0067][0068]
其中y为微分变量，y’表示y的导数，t为时间，t0和t
end
表示起始积分时刻和结束积分时刻，n为方程组维度(即y的维度)，y0表示t0时刻y的值。
[0069]
本发明在线性多步法和runge-kutta法的理论基础上，选取了四种稳定性较好且精度也能满足要求的数值算法的作为底层核心处理方法：
[0070]
1.经典bdf(3～5阶)方法，又称作gear方法；
[0071]
2.改进bdf法，即ibdf(3～5阶)算法；
[0072]
3.3级4阶单对角隐式辛runge-kutta方法，即sdisrk34算法；
[0073]
4.4级4阶单对角隐式辛runge-kutta方法，即sdisrk44算法。
[0074]
在上述方法的基础上，构建一套适用于刚性常微分方程初值问题的数值积分求解器，可以高效、稳定地求解ode-ivp问题，尤其是针对刚性问题具有较好的效果；并对上述四种求解方法均提供变步长机制，可充分发挥其大稳定域的优势。
[0075]
参见图1-图5，仿真中刚性常微分方程初值问题的数值积分求解方法，包括以下操作：
[0076]
第1步：选取积分求解算法、初始积分步长、收敛残差等必要参数，结束后继续执行第2步；
[0077]
其中积分求解算法从bdf算法、ibdf算法、sdisrk34算法和sdisrk44算法中选择；
[0078]
第2步：根据所选取的求解算法进行判断：
[0079]
如采用多步法，则判断当前求解器采用的是bdf或ibdf法(一种改进bdf法)，判断完成后继续执行第3步；
[0080]
如采用单步法，则判断当前求解器采用的是sdisrk34或sdisrk44法，判断完成后跳转执行第5步；
[0081]
第3步：根据上一步的求解算法的选择，构建目标非线性方程组，结束后继续执行第4步；
[0082]
第4步：执行预测步结果的优选策略，利用多种预测步算法分别执行得到多个真实解的预测值，将其分别代入目标非线性方程组，从中选取误差最小的一组解向量作为预测步的解向量；结束后继续执行第5步；
[0083]
第5步：将预测步得到的解向量作为迭代初值，应用d-lm方法求解目标非线性方程组的最小二乘解；
[0084]
在求解时先构建迭代公式，并引入迭代步长的阻尼系数；然后更新迭代步长，判断
是否满足收敛要求；如满足则输出结果作为真实解向量，否则继续重复迭代过程不断得到新的解向量，直至满足收敛误差要求；
[0085]
第6步：执行变步长积分策略：分别执行提供积分步长的步长寻优机制、判断改变积分步长的切换机制、保证计算结果正确的变步长数值回插机制；
[0086]
判断误差是否满足要求；如满足，执行步骤7；
[0087]
如不满足，则返回步骤1，按照变步长策略重新设置积分步长。
[0088]
第7步：时间轴向前推进一个时间步，判断是否需要结束仿真；
[0089]
如满足，执行步骤8；
[0090]
如不满足，则返回第1步，按照变步长策略重新设置积分步长。
[0091]
第8步：结束仿真计算，存储并显示计算结果。
[0092]
下面给出更具体的步骤描述。
[0093]
a、在上述求解的第2步中，本求解器共提供了四种求解算法，即两种线性多步方法和两种单步高精度方法，分别将其命名为“方法a”～“方法d”。
[0094]
方法a(3～5阶传统bdf)：
[0095][0096]
方法b(3～5阶改进bdf)：
[0097][0098]
方法c(3级4阶单对角隐式辛runge-kutta)：
[0099][0100]
3级4阶单对角隐式辛rk方法的各项系数一般以butcher矩阵的形式给出，如图6所示。
[0101]
方法d(4级4阶单对角隐式辛runge-kutta)：
[0102][0103]
4级4阶单对角隐式辛rk方法表达式与方法c，区别仅在与butcher矩阵的的各项系数。4级4阶单对角隐式辛rk方法的butcher矩阵如图7所示。
[0104]
对于一般方程用“方法a”或“方法b”即可，这两种方法的求解过程计算量较小，求解效率较高。
[0105]
对于一些刚性振荡、随机、离散-连续混合或者复杂非线性等问题可以尝试使用“方法c”或“方法d”，这两种方法的求解过程计算量较大，求解效率相对较低，但好处是其数值稳定性好(“方法c”和“方法d”是a-稳定，“方法a”或“方法b”是a(
ɑ
)-稳定)，同时，这两种
方法都属于“辛格式(symplectic schemes)”，较之“方法a”和“方法b”，其具有更小的数值耗散，因而具备更好的数值保真度。
[0106]
b、目标非线性方程组的构建
[0107]
本发明提供的四种求解算法均为隐式求解，所以都需要构建目标非线性方程组。
[0108]
方法a的目标非线性方程组表达式为：
[0109][0110]
方法b的目标非线性方程组表达式为：
[0111][0112]
在上述非线性方程组中，y
n+k
是待求解的未知数，其余均为已知数据，f(y
n+k
)为目标非线性方程组。
[0113]
方法c和方法d的目标非线性方程组求解的是隐式rk中的各项系数，其表达式为：
[0114][0115]
在上述非线性方程组中，ki是待求解的未知数，其余均为已知数据，f(ki)为目标非线性方程组。
[0116]
c、预测步结果的优选策略
[0117]
在方法a和方法b中会用到预测步计算。预测步中，求解器提供了5种算法用来进行真实解的预估：
[0118]
1.显式euler法作为预测步，得到真实解的预测值x1；
[0119]
2.显示2阶runge-kutta法作为预测步，得到真实解的预测值x2；
[0120]
3.显示3阶runge-kutta法作为预测步，得到真实解的预测值x3；
[0121]
4.显示4阶runge-kutta法作为预测步，得到真实解的预测值x4；
[0122]
5.用上一步结果作为预测步，得到真实解的预测值x5。
[0123]
将上面5种解向量分别代入“目标非线性方程组”，从中选取误差最小的一组解向量x
pre
作为给“校正步”的输入。
[0124][0125]
d、应用d-lm法求解目标非线性方程组
[0126]
将预测步得到的解向量x
pre
作为迭代初值，应用d-lm法求解“目标非线性方程组”的最小二乘解，最终结果即为真实解向量。
[0127]
d-lm法的求解过程为：
[0128]
1)引入目标非线性方程组，基于levenberg-marquardt方法构建迭代公式
[0129]
(j
t
j+ηi)h
lm
＝-j
t
f,g＝j
t
f,η≥0
ꢀꢀꢀ
(10)
[0130][0131]
[0132][0133]
其中f为待求解目标非线性方程组，j为目标非线性方程组f的jacobian矩阵，j
t
为j的转置矩阵，η为梯度下降法(一阶优化算法)和牛顿迭代法(二阶优化算法)之间的调节系数，i为单位矩阵，ρ为信赖域半径，h
lm
为每一个迭代步的收敛步长，v是一个用于计算信赖域半径的过程变量。
[0134]
2)计算阻尼系数
[0135]
牛顿法最突出的优点是收敛速度快，具有局部二阶收敛性，但是，基本牛顿法初始点需要足够“靠近”极小点，否则有可能导致算法不收敛。为了确保计算过程的稳定性，借鉴阻尼牛顿法引入迭代步长阻尼系数。
[0136][0137]
其中λ为待求解的迭代步长阻尼系数，h
lm
为每一个迭代步的收敛步长，x是当前时刻的解向量。
[0138]
3)更新迭代步长，判断是否满足收敛要求，如满足则输出结果，否则继续重复迭代过程不断得到新的解向量x
new
，直至满足收敛误差要求。
[0139]hlm
＝(j
t
j+ηi)-1
(-j
t
f)
ꢀꢀꢀ
(15)
[0140]
x
new
＝x+λh
lm
ꢀꢀꢀ
(16)
[0141]
e、变步长积分策略(参见图8)
[0142]
本发明给出的a-d四种求解方法均允许使用变步长积分策略，该策略能够显著提升积分求解效率；变步长积分策略主要有以下三个机制：
[0143]
第一个是积分步长的寻优机制；该机制中，步长递减为指数递减关系：
[0144][0145][0146]
其中dt
new
为更新后的积分步长，dt为当前积分步长，error为当前计算误差，ε为提前给定的允许误差，y
p
为p阶p步法求解得到的解向量，z
p+1
为p+1阶p+1步法求解得到的解向量；
[0147]
分别计算得到y
p
和z
p+1
，当error≤1.0则表明计算结果符合给定的允许误差ε要求，当前计算结果可以被接受，当前积分步结束；
[0148]
第二个是“步长调节”和“步长稳定”两种状态之间的判断与切换机制
[0149]
每一个积分步长理论上都应给定一个初始积分步长，然后通过公式(17)和(18)迭代得到一个恰当的积分步长。这样的补偿调节机制就是“步长调节”状态。如果每一个积分步长都按照这种机制进行计算，其计算结果一定是正确的，然而其求解效率一定非常慢，为了进一步提升求解效率这里又提出一种“步长稳定”状态。
[0150]
当连续若干个(本求解器认为连续5个)积分步长都按照“步长调节”状态进行计
算，那么下一个时间步长可以无需迭代而是直接采用前面若干个积分步长的平均值进行计算，但是此时需要利用公式(18)进行监测，如此时满足0.75≤error≤1的条件，测直接认为当前积分步完成计算，时间轴可以直接向前推进。
[0151]
第三个是变步长数值回插机制，此机制仅在多步法中被激活生效，单步法中无需执行此机制。采用变步长求解时，由于每个积分步长都不相同，本求解器提供的两种线性多步法(方法a和b)要求前几步的积分步长均相同，否则会引起数值发散，所以此处引入数值回插的方法，其中interp_1d为一维插值函数：
[0152][0153]
上述三种机制对求解器的有益效果分别是：
[0154]
步长寻优机制为求解器提供了恰当的积分步长，既不会过大导致结果失真，也不会太小导致计算次数增加进而拖慢整个求解效率；
[0155]
步长调节和步长稳定两种状态的判断与切换机制，能够合理判断是否需要改变积分步长，如无必要就尽可能不调整积分步长，节省了步长调整带来的额外计算量，进一步提升了求解效率；
[0156]
变步长数值回插机制保证了线性多步法中计算结果的正确性，避免了步长变化带来的结果失真。
[0157]
以上三种机制共同支撑起了变步长积分策略，缺一不可。
[0158]
具体的，本发明的求解器采用c++作为开发语言，依赖开源线性代数计算包eigen(v3.3.9)，编译器采用intel数学核心函数库(mkl,math kernel library)优化加速，其基本的代码组织架构如图9所示。
[0159]
下面给出具体的实施例。
[0160]
实施例1：d-lm方法和fsolve命令的效率比较
[0161]
本求解器中的bdf和ibdf是非常高效的两种求解方法，求解非线性方程组是其中不可缺少的一个重要环节，本求解器创新地提出了d-lm方法；为了验证其求解效率，此处选取matlab作为对照验证基准，比较本求解器给出的基于c++开发的d-lm方法和matlab中求解非线性方程组的fsolve命令之间的差别。
[0162]
本算例选取如下的非线性方程组作为被测试对象。
[0163][0164]
上述方程的解向量初值为[0,0,0]，分别采用本求解器的d-lm和matlab的fsolve进行求解，计算结果分别如图10和图11所示。
[0165]
根据上面的计算结果可以看出，针对同样的数学模型得到相同的解，即x＝[10.0449,4.9796,-0.0653]，正确性得到了验证。
[0166]
此外，可以看出matlab中的fsolve需要迭代9步可以使方程组达到1e-12求解精
度，本求解器提供的d-lm方法只需要7步就可以达到1e-12求解精度，更少的迭代次数意味着更高效率的可能性。
[0167]
图12所示为重复20次求解计算时，本发明的d-lm方法和fsolve的求解耗时统计。由计算结果可以看出，d-lm方法和fsolve在同一水平，但fsolve命令在最开始计算时效率较低，重复多次后效率会有所提升。
[0168]
实施例2：典型刚性问题实例
[0169]
在化学反应中经常会遇到下面类型的反应速度方程，其中y1，y2，y3表示反应物的浓度，该模型是典型的刚性问题，常被用于各种数值求解方法的仿真验证。
[0170][0171]
利用matlab/simulink商业工具求解上面的常微分方程，可以得到图13所示的结果。以此为比较基准，利用本求解器中的4步4阶改进bdf算法求解上面的常微分方程，将得到的结果与simulink结果做差得到如图14所示的误差曲线，可以看出两者一致性非常好。
[0172]
当前算例中的变步长收敛误差为1e-5，所以得出的解的误差也基本围绕在这一量级。如果希望得到更准确的解，就需要将变步长的收敛误差进一步收紧，但这肯定会影响计算的效率，计算效率和准确度在求解过程中是一对矛盾，需要综合考虑。
[0173]
除了比较准确度之外，本发明还比较了不同求解方法的效率差别。图15所示为国外主流商业工具求解效率和本求解器计算耗时的比较。
[0174]
可以看出，simulink和simulationx的计算耗时分布在0.4～2.0s，本求解器的bdf方法计算耗时主要分布在0.7～1.0s。总的来说，本求解器和国外主流商业工具在计算效率上基本相当。其中sdisrk方法本质上是一种隐式runge-kutta法，由于其本身算法的计算量较大，这也就导致了计算耗时较多。
[0175]
实施例3：典型hamilton系统振荡问题
[0176]
hamilton系统是典型的动力学系统，在天体运行和机构运动等领域被广泛应用于描述系统动力学过程。
[0177]
kepler方程是典型的hamilton系统，其表达式为：
[0178][0179]
其hamilton函数为：
[0180][0181]
当其积分初值为y1(0)＝0,y2(0)＝1,y3(0)＝1,y4(0)＝0时，其解析解的数学表达式为：
[0182]
y1＝-sin(t)，y2＝cos(t),y3＝cos(t),y4＝sin(t)
[0183]
利用本求解器对公式(22)、(23)所描述的震荡方程进行数值积分求解，得到的结果如图16和图17所示。
[0184]
可以看出数值解和理论解析解有很好的一致性，两者的误差为1e-5量级，这与求解器的预先设置有关，但更准确的数值解需要更多的计算耗时。
[0185]
实施例4：lorenz混沌系统
[0186]
lorenz系统是典型的混沌系统，其动力学过程比较复杂，是检验微分方程求解算法正确性的典型手段。其微分方程的表达式为：
[0187][0188]
利用本求解器对上述lorenz混沌系统进行数值求解，并将结果绘制成三维演化轨迹与matlab/simulink计算结果进行比对参见图18，可以看出两者计算结果所体现的混沌吸引子完全一致，可证明本求解器的正确性。
[0189]
实施例5：r-c-r一维流体管路系统
[0190]
在多学科系统仿真中，管路模型中主要考虑流阻(r)和流容(c)，流阻主要体现管路的阻性，流容主要体现管路的压缩性。当管路较长时，可将其分段处理，对应的模型可用r-c-r-···-c-r表示。
[0191]
对于一个简单的等截面圆形管路，可简单将其等效为r-c-r模型(如图19所示)，其动力学模型为：
[0192][0193]
其中p
in,out
为管路进出口压力，为管路质量流量，r为管路流阻，ρ为流体介质密度，η为流量系数，a为管路内径，v为管路内容积，a为流体介质声速。计算时，取介质密度ρ＝1000kg/m3，管路长度12.73m，管路内径10mm，介质声速a＝1400m/s。
[0194]
计算得到的管路左右两端的压力如图20所示，在相同的管路模型和压力边界条件下，应用本求解器和simulink计算并比较两者的计算结果，如图21所示。可以看出，本求解器的计算结果与simulink的计算结果完全一致。
[0195]
以上给出的实施例是实现本发明较优的例子，本发明不限于上述实施例。本领域的技术人员根据本发明技术方案的技术特征所做出的任何非本质的添加、替换，均属于本发明的保护范围。

技术特征：
1.一种仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，包括以下操作：1)针对待求解的问题，选取包括积分求解算法、初始积分步长、收敛残差在内的必要参数；其中积分求解算法从bdf算法、ibdf算法、sdisrk34算法和sdisrk44算法中选择；2)根据所选取的求解算法进行判断：如采用多步法，则再判断当前求解器采用的是bdf或ibdf法，判断完成后继续执行3)；如采用单步法，则再判断当前求解器采用的是sdisrk34算法还是sdisrk44算法，判断完成后跳转执行5)；3)根据选取的求解算法，构建目标非线性方程组；4)利用多种预测步算法分别执行得到多个真实解的预测值，将其分别代入目标非线性方程组，从中选取误差最小的一组解向量作为预测步的解向量；5)将预测步得到的解向量作为迭代初值，求解目标非线性方程组的最小二乘解；在求解时先构建迭代公式，并引入迭代步长的阻尼系数；然后更新迭代步长，判断是否满足收敛要求；如满足则输出结果作为真实解向量，否则继续重复迭代过程不断得到新的解向量，直至满足收敛误差要求；6)执行变步长积分策略：分别执行提供积分步长的步长寻优机制、判断改变积分步长的切换机制、保证计算结果正确的变步长数值回插机制；变步长积分策略执行完成后，判断误差是否满足要求；如满足则执行7)；如不满足则返回1)，并按照变步长策略重新设置积分步长；7)时间轴向前推进一个时间步，判断是否需要结束仿真；如满足则执行8)；如不满足则返回1)，按照变步长策略重新设置积分步长；8)结束仿真计算，存储并显示计算结果。2.如权利要求1所述的仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，在积分求解算法选取时，对于一般方程采用bdf算法或ibdf算法；对于刚性振荡、随机、离散-连续混合或者复杂非线性的问题则使用sdisrk34算法或sdisrk44算法。3.如权利要求1所述的仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，所述目标非线性方程组的构建为：bdf算法的目标非线性方程组表达式为：ibdf算法的目标非线性方程组表达式为：在上述非线性方程组中，y
n+k
是待求解的未知数，其余均为已知数据，f(y
n+k
)为目标非线性方程组；sdisrk34算法和sdisrk44算法的目标非线性方程组求解的是隐式rk中的各项系数，其
表达式为：其中，k
i
是待求解的未知数，其余均为已知数据，f(k
i
)为目标非线性方程组。4.如权利要求1所述的仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，所述预测步算法通过以下算法用来进行真实解的预估：401)显式euler法作为预测步，得到真实解的预测值x1；402)显示2阶runge-kutta法作为预测步，得到真实解的预测值x2；403)显示3阶runge-kutta法作为预测步，得到真实解的预测值x3；404)显示4阶runge-kutta法作为预测步，得到真实解的预测值x4；405)用上一步结果作为预测步，得到真实解的预测值x5；将上面5种解向量分别代入目标非线性方程组，从中选取误差最小的一组解向量x
pre
作为给校正步的输入：5.如权利要求1所述的仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，采用以下操作求解目标非线性方程组的最小二乘解：501)引入目标非线性方程组，基于levenberg-marquardt方法构建迭代公式：(j
t
j+ηi)h
lm
＝-j
t
f,g＝j
t
f,η≥0f,η≥0if ρ＞0η＝η*max{1/3,1-(2ρ-1)3}；ν＝2elseη＝η*ν；ν＝2*ν其中f为待求解目标非线性方程组，j为目标非线性方程组f的jacobian矩阵，j
t
为j的转置矩阵，η为梯度下降法(一阶优化算法)和牛顿迭代法(二阶优化算法)之间的调节系数，i为单位矩阵，ρ为信赖域半径，h
lm
为每一个迭代步的收敛步长，v是一个用于计算信赖域半径的过程变量；502)引入迭代步长阻尼系数：其中，λ为待求解的迭代步长阻尼系数，h
lm
为每一个迭代步的收敛步长，x是当前时刻的解向量；503)更新迭代步长，重复迭代过程不断得到新的解向量x
new
，直至满足收敛误差要求，并输出结果：h
lm
＝(j
t
j+ηi)-1
(-j
t
f)x
new
＝x+λh
lm
。
6.如权利要求1所述的仿真中刚性常微分方程初值问题的数值积分求解方法，其特征在于，所述变步长积分策略中各机制的执行为：601)寻优机制，步长递减为指数递减关系：01)寻优机制，步长递减为指数递减关系：其中dt
new
为更新后的积分步长，dt为当前积分步长，error为当前计算误差，ε为提前给定的允许误差，y
p
为p阶p步法求解得到的解向量，z
p+1
为p+1阶p+1步法求解得到的解向量；若error≤1.0则表明计算结果符合给定的允许误差ε要求，当前计算结果可以被接受，当前积分步结束；若error>1.0则表明当前计算结果不满足给定的允许误差要求，应重新计算新的dt
new
；602)判断改变积分步长的切换机制，是在步长调节和步长稳定两种状态之间的判断与切换的机制：所述步长调节状态为：每一个积分步长上都给定一个初始积分步长，然后通过(17)式和(18)式迭代得到一个积分步长；所述的步长稳定状态为：当连续若干个积分步长都按照步长调节状态进行计算，那么下一个时间步长无需迭代而是直接采用前面若干个积分步长的平均值进行计算，并利用(18)式进行监测；监测时若满足0.75≤error≤1的条件，则直接认为当前积分步完成计算，时间轴直接向前推进；603)变步长数值回插机制：该机制仅在多步法中被激活生效，引入数值回插的方法，其中interp_1d为一维插值函数：

技术总结
本发明公开了一种适用于刚性常微分方程初值问题的数值积分求解方法，在保证求解精度基础上提出了d-LM方法，进一步提高了隐式运算稳定性。在填补国内空白的基础上，其运算效率接近国外同类产品；与国外商用系统仿真软件工具相比，本发明的数值方法更新颖，对1990年代之后的新的数值求解方法进行了实现与应用，数值求解过程更加稳定。值求解过程更加稳定。


技术研发人员：王珺 郝康康 付翔 尚永权 吕文军
受保护的技术使用者：西安中锐创联科技有限公司
技术研发日：2021.12.06
技术公布日：2022/3/8