ikeda_lyap.c

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

池田写像のリアプノフスペクトラムの計算

calculating the Lyapunov spectrum of the Ikeda map

x = 1.0 + MYU * ( x * cos(theta) - y * sin(theta) )

y = MYU * ( x * sin(theta) + y * cos(theta) )

theta = A - B/(double)( 1.0 + x * x + y * y )

ikeda_lyap.c

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

#include "ikeda_lyap.h"

int main(void)

{

int i;

int j;

int t;

double x;

double y;

double d2;

double lyap_dim;

double lyap[DIM]; /* リアプノフ指数, Lyapunov exponent */

double u[DIM][DIM];

double e[DIM][DIM];

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

// 直交化した長さ１のベクトル（単位行列）uの作成

// create orthogonalized vectors u with length 1 (unit matrix)

for ( i = 0; i < DIM; i++ ) {

for ( j = 0; j < DIM; j++ ) {

u[i][j] = ( i == j ? 1.0 : 0.0 );

// printf("%lf   ", u[i][j]);

}

// printf("\n");

lyap[i] = 0.0;

}

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

// 初期値の設定

// set the initial values

x = X0;

y = Y0;

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

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

// 池田写像の値とヤコビ行列の計算

// calculate the Ikeda map and the Jacobian matrix

ikeda_eqs( &x, &y, u );

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

}

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

// uの再設定

// reset vector u

for ( i = 0; i < DIM; i++ ) {

for ( j = 0; j < DIM; j++ ) {

u[i][j] = ( i == j ? 1.0 : 0.0 );

}

}

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

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

// 池田写像の値とヤコビ行列の計算

// calculate the Ikeda map and the Jacobian matrix

ikeda_eqs( &x, &y, u );

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

// グラム・シュミットの正規直交化法

// Gram-Schmidt orthogonalisation

gram_schmidt_orth( u, e );

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

// calculate the Lyapunov exponents

for ( j = 0; j < DIM; j++ ) {

d2 = 0.0;

for ( i = 0; i < DIM; i++ ) {

d2 += e[i][j] * e[i][j];

}

lyap[j] += log( sqrt( d2 ) );

// lyap[j] += log2( sqrt( d2 ) );

// printf( "%.10lf ", lyap[j]/(double)(t+1) );

}

// printf("\n");

}

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

for ( i = 0; i < DIM; i++ ) {

lyap[i] = lyap[i]/(double)ITERATION;

}

// リアプノフ次元の計算

// calculate the Lyapunov dimension

lyap_dim = calc_lyap_dim( lyap );

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("#        y0: %g\n", Y0);

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

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

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

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

printf("#   lyap[1]:  %lf\n", lyap[0] );

printf("#   lyap[2]: %lf\n", lyap[1] );

printf("# Lyapunov Dimension: %lf\n", lyap_dim );

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

d2 = lyap[0] + lyap[1];

printf("#   lyap[1]+lyap[2]: %lf\n", d2 );

printf("#   2.0 * ln( MYU ): %lf\n", 2.0 * log( MYU ) );

// printf("# 2.0 * log2( MYU ): %lf\n", 2.0 * log2( MYU ) );

return 0;

}

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

double log2( double x )

{

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

}

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

// リアプノフ次元の計算

// calculate Lyapunov dimension

double calc_lyap_dim( double lyap[DIM] )

{

int i;

double sum;

double lyap_dim = (double)DIM;

sum = 0.0;

for ( i = 0; i < (int)DIM; i++ ) {

sum += lyap[i];

if ( sum < 0.0 ) {

lyap_dim = i + (sum - lyap[i])/(double)fabs( lyap[i] );

break;

}

}

return lyap_dim;

}