logistic_lyap.c

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

ロジスティック写像のリアプノフ指数の計算

calculating the Lyapunov exponent of the Logistic map

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

logistic_lyap.c

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

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#define ITERATION 100000000

#define SKIP      10000

#define A  4.0

#define X0 0.121 /* 初期値, initial value */

double log2( double x );

int main(void)

{

int t;

double x;

double lyap;

x = X0;

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

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

x = A * x * ( 1.0 - x );

}

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

lyap = 0.0;

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

// ロジスティック写像の値の計算

// calculate the Logistic map

x = A * x * ( 1.0 - x );

// printf("%.10lf\n", x); /* output the data */

// 線形化したロジスティック写像を利用してリアプノフ指数の計算

// calculate the Lyapunov exponents using the linearized Logistic map

lyap += log( fabs( A - 2.0 * A * x ) );

// lyap += log2( fabs( A - 2.0 * A * x ) );

// printf("%d   %.10lf\n", t, lyap/(double)(t + 1));

}

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

// リアプノフ指数の計算

// calculate the Lyapunov exponents

lyap = lyap/(double)ITERATION;

printf("#----------------------------------------------\n" );

printf("# ITERATION: %d\n", (int)ITERATION);

printf("#      SKIP: %d\n", (int)SKIP);

printf("#----------------------------------------------\n" );

printf("# initial condition\n");

printf("#        x0: %g\n", X0);

printf("#----------------------------------------------\n" );

printf("# A = %g\n", (double)A );

printf("#----------------------------------------------\n" );

printf("# Lyapunov exponent (log with base-e)\n" );

printf("#          %.10lf\n", lyap );

printf("#----------------------------------------------\n" );

printf("#   ln(2): %.10lf\n", log( 2.0 ) );

printf("#----------------------------------------------\n" );

return 0;

}

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

double log2( double x )

{

return log(x)/(double)log(2.0);

}