Lorenz振子的快速求解方法
技术领域
本发明涉及数字信号处理领域，具体涉及一种Lorenz振子的快速求解方法。
背景技术
Lorenz振子是一种重要的混沌振子，该振子广泛应用于工程领域，比如文献 “基于Lorenz系统Lyapunov指数的管道超声导波检测”、“基于Lorenz混沌同步 系统的未知频率微弱信号检测”、“基于Lorenz系统的微弱谐和信号检测”等。
在应用Lorenz振子时，必须进行有效的数值求解。由于Lorenz方程是非刚 性方程，所以通常利用显式的龙格库塔法进行求解，尤其是定步长四阶龙格库塔 法兼具高精度和小运算量的优点，所以是首选数值求解方法。但是定步长四阶龙 格库塔法在求解一个数据点时，需要将微分方程求解4次，也就是说N个数据 点的求解需要将微分方程求解4N次，可见该方法的运算速度是有提升空间的。 因此，研究数值求解方法具有实际的意义。
发明内容
针对上述问题，本发明提供了一种基于半隐式方法的数值计算方法，可以将Lorenz振子求解时应用的定步长四阶龙格库塔法的运算时长减少四分之一。
本发明的技术方案是：Lorenz吸引子动力学系统的快速数值求解方法，具体 步骤包括如下：
步骤一、将Lorenz吸引子动力学系统方程变量的前一递推值，执行两个并 行计算步骤1a、1b，分别得到方程变量的下一递推值后，执行一个平均计算步 骤1c，得到方程变量的下一递推值1；
步骤二、将Lorenz吸引子动力学系统方程变量的前一递推值，执行两个并 行计算步骤2a、2b，分别得到方程变量的递推中间值后，执行一个平均计算步 骤2c，得到方程变量的递推中间值；然后将步骤2c得到的方程变量的递推中间 值，执行两个并行计算步骤2d、2e，分别得到方程变量的下一递推值后，执行 一个平均计算步骤2f，得到方程变量的下一递推值2；
步骤三、将步骤一求出的方程变量下一递推值1和步骤二求出的相应方程变 量下一递推值2带入一个公式进行计算，可得方程变量下一递推值。
进一步的，所述的Lorenz吸引子动力学系统方程具体如下式所示：
式(1)中，x、y、z表示该方程变量，分别表示该方程变量的一 阶导数，σ、r、b表示常数的系统参数。
进一步的，在步骤一中，将步骤一划分为1a、1b、1c三个子步骤，其中，
在步骤1a中，对于Lorenz吸引子动力学系统方程，将方程变量的前一递推 值带入公式(2)中求解得到下一递推值所述公 式(2)如下式所示：
在步骤1b中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(3)中求解得到方程向量下一递推值所述公式 (3)如下式所示：
在式(2)和(3)中，变量的上标n表示变量的第n个数据，变量的下标1 表示递推过程的第一步；变量的上标n+1表示变量的第n+1个数据，变量的下 标1表示递推过程的第一步；h表示递推步长H；另外，步骤1a和1b是单独并 行计算的；
在步骤1c中，将步骤1a、步骤1b求解出的变量下一递推值求平均，得到 方程变量下一递推值1。
进一步的，在步骤二中，将步骤二划分为2a、2b、2c、2d、2e及2f六个子 步骤，其中，
在步骤2a中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(4)中求解得到方程变量递推中间值
所述公式(4)如下式所示：
在步骤2b中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(5)求解得到方程变量递推中间值所述公式 (5)如下式所示：
在公式(4)和(5)中，变量的上标n表示变量的第n个数据，变量的下标 2表示递推过程的第二步；变量的上标n+1/2表示变量的第n+1/2个中间变量数 据，变量的下标2表示递推过程的第二步；h表示递推步长H/2；
另外，步骤2a、2b步骤是单独并行计算的；
在步骤2c中，将步骤2a、2b求解出的方程变量递推中间值求平均，得到方 程变量递推中间值；
在步骤2d中，将步骤2c中求得的中间变量数据值带 入公式(6)求解得到方程变量下一递推值所述公式(6)如 下式所示：
在步骤2e中，将步骤2c中求得的中间变量数据值带 入公式(7)求解得到方程变量下一递推值所述公式(7)如 下式所示：
在公式(6)和(7)中，变量的上标n+1表示变量的第n+1个数据，变量 的下标2表示递推过程的第二步；变量的上标n+1/2表示变量的第n+1/2个中间 变量数据，变量的下标2表示递推过程的第二步；h表示递推步长H/2；
另外，步骤2d、2e是单独并行计算的；
在步骤2f中，将子步骤2d、2e中求出的对应变量值求平均，得到方程变量 的下一递推值2。
进一步的，在步骤三中，将步骤一求出的方程变量下一递推值1和步骤二求 出的相应方程变量下一递推值2带入公式(8)进行计算，可得方程变量下一递 推值；
所述公式(8)如下式所示，
其中，步骤三的计算结果就是方程变量的下一递推值。
通过以上三个步骤，即可通过变量的第n个数据值递推出第n+1个数据值。 且由于本发明方法的递推算法也是4阶精度的，和4阶龙格库塔法精度一样，其 对Lorenz方程的求解结果精度也一样的。
从以上的半隐式递推步骤可以看出，每个变量的递推需要计算六次微分方程 才能得到，但是由于采用了并行处理方法，其计算时间缩减了一半，所以每个变 量的计算时间为计算三次微分方程的时间。而定步长4阶龙格库塔法，每个变量 的计算时间为计算四次微分方程的时间，所以本发明方法比定步长4阶龙格库塔 法的计算时间缩短了四分之一。
本发明的有益效果是：利用本发明提出的递推算法，提升了Lorenz方程数 字求解的速度，本发明方法相对于定步长4阶龙格库塔法的计算时间缩短了四分 之一。
附图说明
图1是本发明的结构流程图；
图2是通过定步长4阶龙格库塔法求解出的Lorenz方程的三维相图；
图3是本发明中求解出的Lorenz方程的三维相图；
图4是是求解出的Lorenz方程的x变量随时间变化的波形示意图；
图5是本发明中求解出的Lorenz方程的x变量随时间变化的波形图；
图6是本发明实施例中利用定步长4阶龙格库塔法求解出的Lorenz方程的y 变量随时间变化的波形图；
图7是本发明实施例中利用本发明方法求解出的Lorenz方程的y变量随时 间变化的波形图。
具体实施方式
为了更清楚地说明本发明的技术方案，下面结合附图对本发明的技术方案做 进一步的详细说明：
如图1所述；本发明采用半隐式方法构造Lorenz吸引子动力学系统的快速 数值求解方法，具体步骤包括如下：
步骤一、将Lorenz吸引子动力学系统方程变量的前一递推值，执行两个并 行计算步骤1a、1b，分别得到方程变量的下一递推值后，执行一个平均计算步 骤1c，得到方程变量的下一递推值1；
步骤二、将Lorenz吸引子动力学系统方程变量的前一递推值，执行两个并 行计算步骤2a、2b，分别得到方程变量的递推中间值后，执行一个平均计算步 骤2c，得到方程变量的递推中间值；然后将步骤2c得到的方程变量的递推中间 值，执行两个并行计算步骤2d、2e，分别得到方程变量的下一递推值后，执行 一个平均计算步骤2f，得到方程变量的下一递推值2；
步骤三、将步骤一求出的方程变量下一递推值1和步骤二求出的相应方程变 量下一递推值2带入一个公式进行计算，可得方程变量下一递推值。
进一步的，所述的Lorenz吸引子动力学系统方程具体如下式所示：
式(1)中，x、y、z表示该方程变量，分别表示该方程变量的一 阶导数，σ、r、b表示常数的系统参数；取不同的系统参数，振子会表现出丰 富的非线性特征。
对于Lorenz吸引子的动力学方程可以采用定步长4阶龙格库塔法求解，其 递推公式如下：
式中：xn、yn、zn是方程变量x、y、z的龙格库塔法的第n个递推值，xn+1、 yn+1、zn+1是方程变量x、y、z的龙格库塔法的第n+1个递推值，h是递推步长。
从定步长4阶龙格库塔法的递推公式可以看出，每个变量需四次计算微分方 程才能得到。
进一步的，在步骤一中，将步骤一划分为1a、1b、1c三个子步骤，其中，
在步骤1a中，对于Lorenz吸引子动力学系统方程，将方程变量的前一递推 值带入公式(2)中求解得到下一递推值所述公 式(2)如下式所示：
在步骤1b中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(3)中求解得到方程向量下一递推值所述公式 (3)如下式所示：
在式(2)和(3)中，变量的上标n表示变量的第n个数据，变量的下标1 表示递推过程的第一步，如表示变量x在递推过程第一步中的第n个数据值， 该值对应于龙格库塔法中的变量的第n个数据xn；变量的上标n+1表示变量的 第n+1个数据，变量的下标1表示递推过程的第一步，如表示变量x在递推 过程第一步中的第n+1个数据值；h表示递推步长H；其中，步骤1a和1b是单 独并行计算的；
在步骤1c中，将步骤1a、步骤1b求解出的变量下一递推值求平均，得到 方程变量下一递推值1；例如，将1a步骤求出的和1b步骤求出的相加 除2，即可得到变量x的下一递推值
进一步的，在步骤二中，将步骤二划分为2a、2b、2c、2d、2e及2f六个子 步骤，其中，
在步骤2a中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(4)中求解得到方程变量递推中间值
所述公式(4)如下式所示：
在步骤2b中，将Lorenz吸引子动力学系统方程变量的前一递推值 带入公式(5)求解得到方程变量递推中间值所述公式 (5)如下式所示：
在公式(4)和(5)中，变量的上标n表示变量的第n个数据，变量的下标 2表示递推过程的第二步，如表示变量x在递推过程第二步中的第n个数据值； 变量的上标n+1/2表示变量的第n+1/2个中间变量数据，变量的下标2表示递推 过程的第二步，如表示变量x在递推过程第二步中的第n+1/2个中间变量 数据值；h表示递推步长H/2；其中，步骤2a、2b步骤是单独并行计算的；
在步骤2c中，将步骤2a、2b求解出的方程变量递推中间值求平均，得到方 程变量递推中间值；例如，将2a步骤求出的和2b步骤求出的相加除 2，即可得到变量x的变量递推中间值
在步骤2d中，将步骤2c中求得的中间变量数据值带 入公式(6)求解得到方程变量下一递推值所述公式(6)如 下式所示：
在步骤2e中，将步骤2c中求得的中间变量数据值带 入公式(7)求解得到方程变量下一递推值所述公式(7)如 下式所示：
在公式(6)和(7)中，变量的上标n+1表示变量的第n+1个数据，变量 的下标2表示递推过程的第二步，如表示变量x在递推过程第二步中的第n+1 个数据值；变量的上标n+1/2表示变量的第n+1/2个中间变量数据，变量的下标 2表示递推过程的第二步，如表示变量x在递推过程第二步中的第n+1/2 个中间变量数据值；h表示递推步长H/2；其中，步骤2d、2e是单独并行计算的；
在步骤2f中，将子步骤2d、2e中求出的对应变量值求平均，得到方程变量 的下一递推值2；例如，将子步骤2d求出的和步骤2e求出的相加除2， 即得到变量x的下一递推值
进一步的，在步骤三中，将步骤一求出的方程变量下一递推值1和步骤二求 出的相应方程变量下一递推值2带入公式(8)进行计算，可得方程变量下一递 推值；
所述公式(8)如下式所示，
其中，步骤三的计算结果就是方程变量的下一递推值；
例如：变量x的第n个数据值通过步骤一求出和步骤二求出后， 通过公式(8)可以计算得到变量x的第n+1个数据值这样就实现了变量x的递推。
通过以上三个步骤，即可通过变量的第n个数据值递推出第n+1个数据值。
如附图2所示，为定步长4阶龙格库塔法求解出的Lorenz方程的三维相图。 求解的参数为(公式(1)中)σ＝10，r＝28，b＝8/3，此时振子处于混沌态。
如附图3所示，为本发明方法求解出的Lorenz方程的三维相图。求解参数 与附图2的求解参数一致。对比可以发现，附图2和附图3的求解结果完全一致。
如附图4所示，为定步长4阶龙格库塔法求解出的Lorenz方程的x变量随 时间变化的波形。
如附图5所示，为本发明方法求解出的Lorenz方程的x变量随时间变化的 波形；对比可以发现，附图4和附图5的求解结果完全一致。
实施例：
依据“基于Lorenz系统的微弱谐和信号检测”构建的信号检测模型，其 Lorenz振子信号检测系统微分方程为：
式中，x、y、z表示该方程变量，分别表示该方程变量的一阶导数， t表示时间，Fbp表示增益为5dB的带通滤波器，s(t)表示待检测的信号，n(t)表示 高斯白噪声；
该检测模型中，先把混合有噪声的待检测信号进行带通滤波，然后输入 Lorenz振子进行检测；在实验中，将n(t)噪声的功率设置为0.04W，信号设置为 s(t)＝0.002cos(70t)；
将信号和噪声输入公式10中，利用定步长4阶龙格库塔法求解，得到如图 6所示的系统输出y的波形图，从中可以明显看到系统输出波形在1秒左右迅速 衰减为很小的值，说明存在有70弧度/秒的角频率待检测信号；该信号检测结果 与文献“基于Lorenz系统的微弱谐和信号检测”中的检测结果是一致的。
同样，利用本发明的快速数值求解方法，对该信号检测问题进行求解，其求 解结果如图7所示。对比图6和图7，可以看出利用定步长4阶龙格库塔法求解 的结果和本发明方法求解结果几乎完全一样。说明本发明方法可以正确应用于信 号检测等实际工程应用中。
最后，应当理解的是，本发明中所述实施例仅用以说明本发明实施例的原则； 其他的变形也可能属于本发明的范围；因此，作为示例而非限制，本发明实施例 的替代配置可视为与本发明的教导一致；相应地，本发明的实施例不限于本发明 明确介绍和描述的实施例。
