tent_lyap.c

/*-----------------------------------------------------------------

テント写像のリアプノフ指数の計算

calculating the Lyapunov exponent of the Tent map

x[t+1] = 1.0 - | 1.0 - A * x[t] |

x[t+1] = A * x[t] (0 < x < 0.5)

x[t+1] = A * ( 1.0 - x[t] ) (0.5 <= x < 1)

tent_lyap.c

-----------------------------------------------------------------*/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

/*---------------------------------------------------------------------*/

#define ITERATION 100000

#define SKIP 10000

#define A 2.0

/*---------------------------------------------------------------------*/

int main(void)

{

int t;

double k0 = 43597;

double m = 9999991;

double k = k0;

double x;

double x0;

double pre_x; // xの１つ前の値, previous value of x

double lyap;

double temp;

/*------------------------------------------------------------------*/

// A=2のとき、数式をそのまま計算すると、コンピュータで使える浮動小数点が不十分のため、

// 50回ほどの計算で値がゼロになる。それを防ぐため、式を分母と分子に分けて計算を行う。

// When A=2, if we calculate the formula as it is, the value will be zero
// after about 50 calculations, because of the insufficient number of
// floating points available in computers.

// To prevent this, the calculations are performed by dividing

// the formula into the denominator and the numerator.

// x[t+1] = 1.0 - fabs( 1.0 - A * x[t] )

// k/m = m/m - fabs( m/m - A * k/m )

/*------------------------------------------------------------------*/

// 0 < x0 < 1.0

x0 = k/(double)m; // 初期値, initial condition

x = x0;

/*------------------------------------------------------------------*/

for( t = 0; t < SKIP; t++ ) {

k = m - fabs( m - A * k );

x = k/(double)m;

}

/*------------------------------------------------------------------*/

lyap = 0.0;

for( t = 0; t < ITERATION; t++ ) {

pre_x = x;

k = m - fabs( m - A * k );

x = k/(double)m;

// printf("%.10lf %.10lf\n", pre_x, x); // show tent map

// 線形化したテント写像を利用してリアプノフ指数の計算

// テント写像を微分した式 (temp = ・・・) から、テント写像のリアプノフ指数がlog(A)となることは明らかだが、

// この式を備忘録として書いておくことにした。

// calculate the Lyapunov exponents using the linearized Tent map

// From the formula for the derivative of the tent map, it is obvious that the Lyapunov exponent

// of the tent map is log(A), but I write down this formula as a memorandum.

temp = A * ( 1.0 - A * x )/(double)( fabs(1.0 - A * x) );

lyap += log( fabs( temp ) );

// lyap += log2( fabs( temp ) );

}

/*------------------------------------------------------------------*/

lyap /= (double)ITERATION;