lorenz_lyap.c

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

ローレンツ方程式のリアプノフスペクトラムの計算

calculating the Lyapunov spectrum of the Lorenz equations

dx/dt = SIGMA * (y - x)

dy/dt = -y - x * z + R * x

dz/dt = x * y - B * z

lorenz_lyap.c

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

#include "lorenz_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[0] = X0;

v[1] = Y0;

v[2] = Z0;

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

solve_lorenz_eqs_using_RK4( v, u );

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

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

}

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

// 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++ ){

solve_lorenz_eqs_using_RK4( v, u );

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

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

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

// 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("#---------------------------------------------------\n" );

printf("# initial condition\n");

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

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

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

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

printf("# SIGMA = %g, R = %g, B = %g\n", (double)SIGMA, (double)R, (double)B);

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("# Lyapunov Dimension: %lf\n", lyap_dim );

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

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

printf("# diagonal elements of the jacobian: %lf\n", -SIGMA - 1.0 - B );

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;

}