最大李雅谱诺夫指数的数值计算
最大李雅谱诺夫指数的数值计算J.C.Sprott
Department of Physics, University of Wisconsin, Madison, WI 53706, USA
October 15, 1997
(Revised August 31, 2004)
袁来翻译——上海大学——2005/03/22;；Email: liguoy@citiz.net
原文出处：http://sprott.physics.wisc.edu/chaos/lyapexp.htm
通
常通过计算最大的李雅谱诺夫指数来检验系统是否混沌，一个正的最大李雅谱诺夫指数标致了该系统是否混沌，当你已经获得已知产生混沌的方程，这是相对容易做
的。当你仅仅只有试验数据记录时，这种计算的困难程度是几乎不可能的，这里我们不考虑这种情形。计算的一般思路是跟踪两个较近的轨道，计算它们分离的平均
对数率，无论它们分离得多远，其中之一的轨道将最终回到另一条轨道的附近。在每次迭代过程中，一程序将实现这个工作，完整的程序过程如下：
1： 在吸引域内，开始于任何初值条件。
最好开始于吸引子上的一个已知点，这种情形，第2步可以省略掉。
2： 迭代直至轨道落到吸引子上。
这需要自己的判断或通过研究得到的有关该系统的一些知识，对绝大多数系统，只要迭代几百次，我们就认为足够了。通常情况下，除非你选的初值点刚好靠近分岔点，要不在几百次后它不会远离吸引子的。
3：选择（几乎任意的）靠近的点（设两点的距离为d0）
适合的d0的选择是至少要比计算机中的浮点数的1000倍大，如：在（8字节）的双精度的（对这些计算的最低推荐），变量有52－比特的尾数，故，精度是：2-52=2.22x10-16，因此，我们只要有d0=10-12 就足够了。
4： 两轨道向前迭代一次，并计算分离量d1
分离量的计算：两轨道的每个分量上的差的平方再和再开平方，如：对二维的系统有变量x和y，那么分离量是: d =[(xa-xb)2+(ya-yb)2]1/2, 这里下标（a和b）分别标出了两轨道。
5： 以任何方便的底数估计log|d1/d0|
通常是采用自然对数（e为底数），但对映射，在每次少许的迭代中，李雅谱诺夫指数被引用，在这种情形下，我们需要以2做为底数（注：log2x=1.4427logex）。
如果d1近似等于0，那么在估计对数时，将会得到错误，在这种情形下，尝试用更大的d0,如果这还不够，那么你不得不忽略这里发生的值。但是这样做的话，你的李雅谱诺夫指数计算可能有错误。
6： 调整一条轨道，使得在方向d1上的分离量为d0
这也许是最难的一步，且易出错，如（在二维的情形），假如轨道b是要调整的，在一次调整后它的值是(xb1,yb1),那么它将重新初始化为：xb0=xa1+d0(xb1-xa1)/d1 和yb0=ya1+d0(yb1-ya1)/d1.
7： 重复4－6步，计算第5步的平均值
你
或许舍掉些前面的值去保证轨道已经进入了最大的扩展方向，这也是一好的想法去计算动态的平均值作为值是否已经趋于一固定的唯一值，得到可靠的计算，有时
候，结果收敛得很慢，但是，在几千次迭代后，通常能得到两个数的精确估计，这是不错的想法去证明你的结果是独立于初值条件的，d0的值，迭代的次数包含在
平均中，你也许会想检验无界的轨道，由于你将得到数值误差，在这种情形下，李雅谱诺夫指数是没有意义的。
在一般的2维的二次映射中，当寻找混沌解时，其Basic源程序和可执行编码实现(source and executable).
算Hénon映射的李雅谱诺夫指数（－0.5到0.5）Xn+1 = 1 - CXn2 + BXn-1 这里B=0.3,0<C<2,( BASIC source and executable),程序的结果输出如下：
如果系统包含常微分方程（流），而不是差分方程（映射），除了所得指数被迭代步长相除外，程序是一样的。下面给出Lorenz attractor的一个计算lyapunov exponent的算法。（available）。
[ 本帖最后由 yejet 于 2006-11-3 21:29 编辑 ] 翻译得很烂 Numerical Calculation of Largest Lyapunov Exponent
J. C. Sprott
Department of Physics, University of Wisconsin, Madison, WI 53706, USA
October 15, 1997
(Revised August 31, 2004)
The usual test for chaos is calculation of the largest Lyapunov exponent. A positive largest Lyapunov exponent indicates chaos. When one has access to the equations generating the chaos, this is relatively easy to do. When one only has access to an experimental data record, such a calculation is difficult to impossible, and that case will not be considered here. The general idea is to follow two nearby orbits and to calculate their average logarithmic rate of separation. Whenever they get too far apart, one of the orbits has to be moved back to the vicinity of the other along the line of separation. A conservative procedure is to do this at each iteration. The complete procedure is as follows:
1. Start with any initial condition in the basin of attraction.
Even better would be to start with a point known to be on the attractor, in which case step 2 can be omitted.
2. Iterate until the orbit is on the attractor.
This requires some judgement or prior knowledge of the system under study. For most systems, it is safe just to iterate a few hundred times and assume that is sufficient. It usually will be, and in any case, the error incurred by being slightly off the attractor is usually not large unless you happen to be very close to a bifurcation point.
3. Select (almost any) nearby point (separated by d0).
An appropriate choice of d0 is one that is about 1000 times larger than the precision of the floating point numbers that are being used. For example, in (8-byte) double-precision (minimum recommended for such calculations), variables have a 52-bit mantissa, and the precision is thus 2-52 = 2.22 x 10-16. Therefore a value of d0 = 10-12 will usually suffice.
4. Advance both orbits one iteration and calculate the new separation d1.
The separation is calculated from the sum of the squares of the differences in each variable. So for a 2-dimensional system with variables x and y, the separation would be d = [(xa - xb)2 + (ya - yb)2]1/2, where the subscripts (a and b) denote the two orbits respectively.
5. Evaluate log |d1/d0| in any convenient base.
By convention, the natural logarithm (base-e) is usually used, but for maps, the Lyapunov exponent is often quoted in bits per iteration, in which case you would need to use base-2. (Note that log2x = 1.4427 loge x). You may get run-time errors when evaluating the logarithm if d1 becomes so small as to be indistinguishable from zero. In such a case, try using a larger value of d0. If this doesn't suffice, you may have to ignore values where this happens, but in doing so, your calculation of the Lyapunov exponent will be somewhat in error.
6. Readjust one orbit so its separation is d0 in the same direction as d1.
This is probably the most difficult and error-prone step. As an example (in 2-dimensions), suppose orbit b is the one to be adjusted and its value after one iteration is (xb1, yb1). It would then be reinitialized to xb0 = xa1 + d0(xb1 - xa1) / d1 and yb0 = ya1 + d0(yb1 - ya1) / d1.
7. Repeat steps 4-6 many times and calculate the average of step 5.
You might wish to discard the first few values you obtain to be sure the orbits have oriented themselves along the direction of maximum expansion. It is also a good idea to calculate a running average as an indication of whether the values have settled down to a unique number and to get an indication of the reliability of the calculation. Sometimes, the result converges rather slowly, but a few thousand iterates usually suffices to obtain an estimate accurate to about two significant digits. It is a good idea to verify that your result is independent of initial conditions, the value of d0, and the number of iterations included in the average. You may also want to test for unbounded orbits, since you will probably get numerical errors and the Lyapunov exponent will not be meaningful in such a case.
Sample software that implements this procedure while searching for chaotic solutions in general 2-D quadratic maps is available in (DOS) BASIC source and executable code. Sample software that calculates the Lyapunov exponent (-0.5 to 0.5) for the Hénon Map Xn+1 = 1 - CXn2 + BXn-1 for B = 0.3 and 0 < C < 2 is also available in (DOS) BASIC source and executable code. Output from that program is shown below:
If the system consists of ordinary differential equations (a flow) instead of difference equations (a map), the procedure is the same except that the resulting exponent is divided by the iteration step size so that it has units of inverse seconds instead of inverse iterations. An example for the Lorenz attractor is available.
Ref: J. C. Sprott, Chaos and Time-Series Analysis (Oxford University Press, 2003), pp.116-117. http://sprott.physics.wisc.edu/chaostsa/
哈，这本书可以下载下来。谁来下载？ 有一点不明白，为什么计算出的是最大李雅普诺夫指数？N维系统有N个李雅普诺夫指数，这样每个时间步前后的偏差应该是所有李雅普诺夫指数一起作用的结果啊是不是? 我看了这本书了，很不错，可是怎么下载啊，我搞了半天也没有下载下来！！！！
心灯能不能打个包发布出来
心灯能不能把这本书打个包发布出来，这样方便大家下载啊 不是有中文版的吗？ 原帖由 ddb_2005 于 2006-8-21 23:42 发表不是有中文版的吗？
中文的是翻译过来的，翻译的不怎么样 谢谢各位的分享
页:
[1]