henon-heiles_lyap.c

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

エノン・ハイレス方程式のリアプノフスペクトラムの計算

calculating the Lyapunov spectrum of the Henon-Heiles equations

dx/dt =  z

dy/dt =  w

dz/dt = -x - 2 * x * y

dw/dt = -y - x^2 + y^2

henon-heiles_lyap.c

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

#include "henon-heiles_lyap.h"

int main(void)

{

int i;

int j;

int t;

double total_T;

double d2;

double lyap_dim;

double sum_lyap;

double v[DIM]; /* 変数 variables */

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

double u[DIM][DIM];

double e[DIM][DIM];

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

// 初期値の設定

// set the initial values

// v[0]: x    v[1]: y    v[2]: z    v[3]: w

v[0] = X0;

v[1] = Y0;

v[2] = Z0;

v[3] = W0;

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

// 直交化した長さ１のベクトル（単位行列）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;

}

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

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

// エノン・ハイレス方程式の値とヤコビ行列の計算

// calculate the Hyperchaos equations and the Jacobian matrix

solve_henon_heiles_eqs_using_RK4( v, u );

// output the data, v[0]: x, v[1]: y, v[2]: z, v[3]: w

// printf("%.10lf   %.10lf   %.10lf   %.10lf\n", v[0], v[1], v[2], v[3] );

}

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

// 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 Hyperchaos equations and the Jacobian matrix

solve_henon_heiles_eqs_using_RK4( v, u );

// output the data, v[0]: x, v[1]: y, v[2]: z, v[3]: w

// printf("%.10lf   %.10lf   %.10lf   %.10lf\n", v[0], v[1], v[2], v[3] );

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

// Gram-Schmidt orthogonalisation

gram_schmidt_orth( u, e );

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");

}

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

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

// calculate the Lyapunov exponents

sum_lyap = 0.0;

total_T = ITERATION * DELTA_T;

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

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

sum_lyap += lyap[i];

}

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

// 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("#    DELTA_T: %g\n", (double)DELTA_T);

printf("# TOTAL_TIME: %g\n", total_T);

printf("# initial condition\n");

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

printf("#         y0: %#g\n", (double)Y0);

printf("#         z0: %#g\n", (double)Z0);

printf("#         w0: %#g\n", (double)W0);

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

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

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

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

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

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

printf("#    lyap[4]: %lf\n", lyap[3] );

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

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

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

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

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;

}