高维非线性映射系统的不稳定流形计算方法研究

贾 蒙

(新乡学院 机电工程学院,河南 新乡 453003)

高维非线性映射系统的不稳定流形计算方法研究

贾 蒙

(新乡学院 机电工程学院,河南 新乡 453003)

提出了一种基于曲率约束条件的计算离散动力系统鞍型不动点一维不稳定流形的新算法，并以Hénon映射为例进行了计算。新算法以增长流形为基本思想，通过曲率约束和距离控制来确定离散点间的距离；提出流形的偏转角度可以通过流形上的已知点来预测，解决了流形上新点的原像位置快速确定的困难。仿真发现：Hénon映射的一维不稳定流形在标准参数下与Hénon映射产生的散点图分布一致，在其它几组参数下，一维不稳定流形的两个分支之间保持着某种程度的对称性，该研究对Hénon映射的进一步研究打下基础。

离散动力系统;双曲不动点;不稳定流形;Hénon映射;混沌

动力系统按照其状态变量随时间的变化是否连续可分为两类：连续动力系统和离散映射动力系统。其中离散动力系统可看作是对连续动力系统采样后的结果，而离散动力系统经常以迭代函数的形式出现。科学研究及生活实际中的很多过程都可以用迭代函数的形式来描述，Hénon映射[1]就是一个比较著名的例子，它具有一定的混沌行为，在图像加密[2]、分形研究[3]等领域都有广泛地应用。其表达式为

(1)

将式(1)写成映射函数的形式

(2)

无论是连续时间系统，还是离散时间系统，双曲平衡点或双曲周期轨道的稳定、不稳定流形的计算对动力学问题的全局性理解起到了非常重要的作用。 对于一个给定的非线性动力学问题，进行严格的理论分析是非常复杂和艰难的，通常采用数值方法计算其稳定和不稳定流形。

Krauskopf等[4]提出利用测地线距离作为参数，不需要对向量场进行调整，而是通过求解一系列边值问题以计算流形。 它能够很好的控制流形计算中的增长问题，但是边值问题的求解计算量很大，特别当与向量场近似相切的时候所需要的步骤会更多，耗时严重． 另外，该方法可以以径向的轨线构建流形，但不能绘制出动力学系统的流． 其他方法，如 PDE 方法[5]、轨道参数化方法[6]、细分法[7]等，在计算流形过程中没有给出有效控制流形各个方向增长的方法。

本文提出了一种基于曲率约束条件的计算离散动力系统鞍型不动点一维不稳定流形的新算法，解决了流形上新点的原像位置快速确定的困难。

1 流形计算的基本方法

首先我们将研究范围扩展到一般映射。设F:Rn→Rn为一个保向的微分同胚映射函数，

xm+1=F(xm)

(3)

取x0为映射的初始点，则系统满足

Fk+l(x0)=Fk(Fl(x0))

(4)

式中:k,l为任意正整数;当F可逆时;k,l可以取任意整数。

对于映射式(3)，若存在x0∈Rn，满足F(x0)=x0，则称x0为系统的不动点。在x0处对系统进行线性化，将式(3)转化为线性映射x→Ax，其中A为F在x0处的雅各比矩阵A=DF(x0)=[∂fi/∂xj](x0)。若矩阵A的特征值的模都不等于1，那么x0就是一个双曲不动点。其中模小于1的特征值叫做稳定特征值，其对应的特征向量{v1,v2,…,vl}张成稳定特征空间Es；模大于1的特征值叫做不稳定特征值，它们对应的特征向量{vl+1,vl+2,…,vn}张成不稳定特征空间Eu。全空间

Rn=Es⊕Eu

(5)

设x0是微分同胚映射函数F的一个双曲不动点，则在x0的邻域U内存在局部稳定和不稳定流形：

(6)

(7)

定理的证明见文献[4]。全局稳定和不稳定流形定义为

稳定流形和不稳定流形在分析系统的动力学特性中起着非常重要的作用，它们充当不同吸引子吸引域的边界，将全空间划分为多个具有不同动力特性的不变子空间，而当稳定与不稳定流形相交时，就会引起同宿、异宿以及混沌[3-6]等复杂动力学行为的出现。不变流形计算对于上述问题的研究都有着积极地促进作用。

离散映射动力系统的流形计算比较困难，对于一维流形的计算，比较著名的算法见文献[7-8]。Krauskopf等[9-10]通过每步向流形中加入一个离散点来增长流形，离散点间的距离受流形局部曲率的控制，该算法能够有效展示流形的细节。该算法的不足之处在于，它用二分法搜索已有流形区间段来确定新离散点的前像位置，比较繁琐，这也是影响该算法计算速度的一个关键因素。

本文提出的算法从增长流形[11-12]这个思想出发的，提出了一种预测流形偏转角度的方法，新方法可以快速确定流形上新离散点的原像位置，同时采用了中曲率控制[13-14]和距离控制条件，算法简洁明了，计算速度快，且易于编程实现。

2 一维流形计算方法的误差分析

贾蒙介绍的各种一维流形计算方法，都没有对算法的计算精度进行讨论。因为计算映射的一维流形时，经常会遇到混沌映射，而混沌的一个显著特点就是对于初始误差的敏感性，意味着计算误差将随计算的进行而指数放大，对于研究流形计算方法的人来说，这不啻为一场灾难。

既然误差会指数放大，所以采用精度有限、具有截断误差的计算机来，计算流形似乎是不可能的。但其实没有这么绝望。混沌系统虽然对于初值具有敏感依赖性，长期行为近似随机，但短期行为还是可以预测的，也就是说，对长期行为进行计算的话误差较大，但计算短期行为时误差仍然是保持有界的。而我们在计算映射系统的一维流形时，计算的弧长都非常有限，因此迭代的次数也很少，所以其误差能够保持在有界范围内。

并且，一切形式为

L={(x,y):|x|<δ,y=p(x)}

的流形L的迭代都将在C1拓扑下收敛于不稳定流形，即

(10)

设M是Rn空间中的一个m维子流形，在M的每个点x∈M处定义维数为n-m的超平面N(x)，该超平面在x处与M横截相交。不失一般性，假设N(x)是光滑的，而且N(x)的选择并不唯一。集合E={N(x),x∈M}在M上具有向量纤维丛(vector fiber bundle)结构，因此可以把M看作零截面(null section){(x,0)}。定义(x⊕y)=x+y，其中x∈M、y∈N(x)。对于任意固定的小量ε>0，称集合Uε={(x⊕y):x∈M,y∈N(x),|y|<ε}为子流形M的一个管状邻域。

设Nε是流形M的一个管状邻域，流形M1⊂Nε的形式为

M1={x⊕h(x):h(x)∈N(x)}

其中x⊕h(x):M→Rn是一个光滑映射。满足该条件的M1被称作是处于M的C1邻域内。C1距离定义为映射的C1范数(C1norm)

(11)

令原点0是微分同胚F:Rn→Rn的一个双曲不动点，则任意形式为

L={(x,y):|x|<δ,y=p(x)}

的流形L将C1拓扑收敛于不稳定流形Wu的紧致部分，即当m→∞时，Wm=Fk(L)∩Nε→M。迭代将逐点收敛于全局不稳定流形Wu。

当用逆函数F-1替换函数F时，以上理论也可适用于稳定流形。

3 基于曲率约束条件的流形计算方法

考虑轨道上相邻三个离散点之间的角度，如图2所示，由图易得

(12)

其中

(13)

(14)

检查α是否满足下列约束条件

αmin<α<αmax

(Δα)min<Δkα<(Δα)max

(15)

将pk+1加入到序列M中，并判断M′的长度arcl是否已经达到所需弧长ARC，满足则结束，不满足则继续计算。

上述步骤都结束后，再沿A的不稳定特征向量的反方向计算Wu(x0)的另一分支，两个分支之和就是Wu(x0)。

4 高维不稳定流形计算与仿真

4.1三维Hénon映射的一维不稳定流形

三维Hénon映射的表达式为

(16)

当M1=1.4,M2=0.2,B=0.1时，具有与二维Hénon映射类似的混沌吸引子，如图2(a)所示。易得该映射具有不动点

x0=(x*,y*,z*)，其中

该不动点是二次映射F2的一个双曲不动点，且具有一维不稳定流形，计算结果如图2(b)所示。

(a) 混沌吸引子

(b) 一维不稳定流形，弧长为30，计算耗时0.6 s

观察比较图2中的两幅图，发现二者几乎一模一

样，不同之处在于：(a)图中的点是混沌的，近似随机排列的，如果按照迭代的顺序将这些点顺序连接起来的话，整个图将是杂乱无章的，而(b)图中的点是有序的。我们发现这点同样适用于二维Hénon映射。

当M1=0.19、M2=0.999 1、B=0时,x0也是一个双曲不动点，一维流形计算结果如图3所示。

图3 三维Hénon映射的一维不稳定流形

Fig.3 One dimensional unstable manifold of three dimensional Hénon system

图3中间似乎“杂乱”却又左右对称的部分为三维Hénon映射的吸引子，可见该吸引子完全被一维不稳定流形所“包围”。

4.2四维映射的一维流形计算

为了说明新算法能够用于高维映射的流形计算，特将三维Hénon映射增加了一维，变为四维Hénon映射，表达式为

(17)

(a) w=0

(b) z=0

(c) y=0

(d) x=0

图4 四维Hénon map的一维不稳定流形

Fig.4 One dimensional unstable manifold of four dimensional Hénon system

5 总 结

本文提出了一种计算离散动力系统不动点一维不稳定流形的新算法，该算法通过每步加入一个离散点来增长流形，相邻离散点间的距离通过曲率约束和距离控制条件来确定。提出了一种不稳定流形上角度的预测方法，可以用来快速地确定当前所要加入点的原像的位置，这是相对于文献[9]中提出的算法的优越之处。此外，本算法可同时用于一维稳定和不稳定流形的计算，并且导数传递的推导在高维空间下依旧成立，算例5.1和5.2说明了本章算法可以用于高维空间下的一维流形计算。在计算一维稳定流形时，本章算法不需要映射的逆函数F-1的显式表达式，所以通用性也较传统算法有了提高。

本算法需要计算映射函数的Jacobian矩阵，所以得考虑该矩阵的计算难度。对于显式定义的同胚映射函数，其Jacobian矩阵可以写成显式表达式，对于某个特定点，只要将坐标带入即可求得，计算比较简单，耗时与函数进行一次迭代相当；但对于向量场的Poincaré映射，其Jacobian矩阵需要用数值方法进行求解，比较复杂，在这种情况下，本算法可能是不具优势的。对于非同胚映射，或者函数不可逆以及有多个逆，本算法都不能够使用。

另外，当公式(1)参数a和b取其它更多变化时，Hénon映射的不稳定流形可能会表现出更加丰富的性质以及更加漂亮的图形，这有待我们进一步的探索验证。我们希望本文的计算结果能对Hénon映射的进一步研究产生一些启发。

[1] 贾蒙. 映射动力系统一维流形并行计算方法[J]. 振动与冲击, 2014,33(9):40-47.

JIA Meng. A parallel algorithm for approximating 1-D manifold of map[J].Journal of Vibration and Shock, 2014,33(9):40-47.

[2] GAO H, ZHANG Y, LIANG S, et al. New chaotic algorithm for image encryption[J].Chaos Solutions and Fractals,2006(29): 393-399.

[3] GRASSBERGER P. On the fractal dimension of the Henon attractor[J]. Physics Letters A, 1983, 97(6):224-226.

[4] PALIS J, MELO W D. Geometric theory of dynamical Systems[M].Springer, 1982.

[5] YOU Z, KOSTELICH E J, YORKE J A. Calculating stable and unstable manifolds[J]. Int. J. Bifurc. Chaos Appl.Sci. Eng.,1991, 1(3): 605-623.

[7] PARKER T S, CHUA L O. Practical numerical algorithms for Chaotic Systems[M]. Springer, Berlin, 1989.

[8] HOBSON D. An efficient method for computing invariant manifolds[J]. J. Comput. Phys. 1991,104: 14-22.

[9] KRAUSKOPF B, OSINGA H M. Growing unstable manifolds of planar maps[J]. 1517, 1997, http://www.ima.umn.edu/preprints/OCT97/1517.ps.gz.

[10] KRAUSKOPF B, OSINGA H M. Growing 1D and quasi-2D unstable manifolds of maps[J]. J. Comput. Phys.1998,146: 406-419.

[11] ENGLAND J P, KRAUSKOPF B, OSINGA H M. Computing One-Dimensional Stable Manifolds and Stable Sets of Planar Maps without the Inverse[J]. SIAM J. Appl. Dyn. Syst. 2004,3(2):161-190.

[12] DELLNITZ M, HOHMANN A. A subdivision algorithm for the computation of unstable manifolds and global attractors[J]. Numer. Math, 1997,75: 293-317.

[13] FUNDINGER D. Toward the calculation of higher-dimensional stable manifolds and stable sets for noninvertible and piecewise-smooth maps[J]. J Nonlinear Sci, 2008, 18: 391-413.

[14] KRAUSKOPF B, OSINGA H M. Globalizing two-dimensional unstable manifolds of maps[J]. Int. J. Bifurc.Chaos,1998, 8(3): 483-503.

Computationmethodforunstablemanifoldofhénonmap

JIA Meng

(College of Mechanical & Electrical Engineering， Xinxiang College， Xinxiang 453003, China)

A new algorithm was presented for computing one dimensional unstable manifold for saddle fixed point of a discrete dynamic system based on the constraint condition of curvature. Hénon map was taken as an example to check the performance of the algorithm. The manifold growing was taken as the basic idea of the new algorithm. Curvature constraint and distance control were used to determine the distance between discrete points. The unstable manifold grew with new point added at each step and the distance between consecutive points was adjusted according to the local curvature. It was proved that the gradient of the manifold at the new point can be predicted with the known points on the manifold and in this way the preimage of the new point can be located immediately. The simulation showed that the one dimensional unstable manifold of Hénon map coincides with the scatter point diagram distribution produced by itself under standard parameters; under the other several groups of parameters, two branches of the unstable manifold are nearly symmetric, and they serve as the borderline of Hénon map iteration sequence. The results laid a foundation for further studying Hénon map.

discrete dynamic system; hyperbolic fixed point; unstable manifold; Hénon map; chaos

国家自然科学基金(61501391);河南省高等学校重点科技项目(15A510035)

2016-03-14 修改稿收到日期：2016-07-05

TP301.6

： A

10.13465/j.cnki.jvs.2017.17.038

作 者 贾蒙 男，博士，副教授，1981年12月生