/*------------------------------------------------------------------------
中村・上田写像のリアプノフスペクトラムの計算
calculating the Lyapunov spectrum of the Nakamura-Ueta map
x[t] = A0 + A1 * y[t-2] + A2 * x[t-1] x[t-3]
x[t+1] = A0 + A1 * y[t] + A2 * x[t] z[t]
y[t+1] = x[t]
z[t+1] = y[t] = x[t-1]
nakamura-ueta_lyap.c
------------------------------------------------------------------------*/
#include "nakamura-ueta_lyap.h"
int main(void)
{
int i;
int j;
int t;
double x;
double y;
double z;
double lyap_dim;
double sum_lyap;
double R_diag;
double lyap[DIM]; // リアプノフ指数, Lyapunov exponent
double dfQ[DIM][DIM]; // ヤコビ行列dfと行列Qの積
double Q[DIM][DIM]; // QR分解
double R[DIM][DIM]; // QR分解
/*----------------------------------------------------------------------------*/
// 初期値の設定
// set the initial values
x = X0;
y = Y0;
z = Z0;
/*----------------------------------------------------------------------------*/
// 直交化した長さ１のベクトル（単位行列）の作成
// create orthogonalized vectors with length 1 (unit matrix)
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
Q[i][j] = ( i == j ? 1.0 : 0.0 );
}
lyap[i] = 0.0;
}
/*----------------------------------------------------------------------------*/
// 中村・上田写像の値とそのヤコビ行列の計算
// calculate the Nakamura-Ueta mapl and the Jacobian matrix
for( t = 0; t < SKIP; t++ ){
nakamura_ueta_eqs( &x, &y, &z, dfQ, Q );
// printf("%.10lf %.10lf %.10lf\n", x, y, z);
}
/*----------------------------------------------------------------------------*/
for( t = 0; t < ITERATION; t++ ){
// 中村・上田写像の値とそのヤコビ行列の計算
// calculate the Nakamura-Ueta mapl and the Jacobian matrix
nakamura_ueta_eqs( &x, &y, &z, dfQ, Q );
// printf("%.10lf %.10lf %.10lf\n", x, y, z); /* output the generated data */
// dfQをQR分解
// QR decomposition of dfQ
householder_QR_decomp( dfQ, Q, R );
// 行列Rの対角成分R_diagを使ってリアプノフスペクトラムを計算
// calculate the Lyapunov exponents using the diagonal component of the matrix R
for ( i = 0; i < DIM; i++ ) {
R_diag = R[i][i]; // 行列Rの対角成分を取得, get the diagonal components of the matrix R
lyap[i] += log( fabs( R_diag ) );
}
}
/*----------------------------------------------------------------------------*/
sum_lyap = 0.0;
for ( i = 0; i < DIM; i++ ) {
lyap[i] = lyap[i]/(double)ITERATION;
sum_lyap += lyap[i];
}
// リアプノフ次元の計算
// calculate Lyapunov dimension
lyap_dim = calc_lyap_dim( lyap );
printf("#----------------------------------------------\n" );
printf("# ITERATION: %d\n", (int)ITERATION);
printf("# SKIP: %d\n", (int)SKIP);
printf("# initial condition\n");
printf("# x0: %g\n", X0);
printf("# y0: %g\n", Y0);
printf("# z0: %g\n", Z0);
printf("#----------------------------------------------\n" );
printf("# A0 = %g, A1 = %g, A2 = %g\n", (double)A0, (double)A1, (double)A2 );
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("# lyap[3]: %lf\n", lyap[2] );
printf("# Lyapunov Dimension: %lf\n", lyap_dim );
printf("#----------------------------------------------\n" );
printf("# sum of lyap: %lf\n", sum_lyap);
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;
}