/*------------------------------------------------------------------------------------
ダフィング方程式のリアプノフスペクトラムの計算
calculating the Lyapunov spectrum of the Duffing equations
dx/dt = y
dy/dt = -x^3 - K * y + B * sin(z)
dz/dt = 1.0
duffing_lyap.c
------------------------------------------------------------------------------------*/
#include "duffing_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 < SKIPS; t++ ){
// ダフィング方程式の値とヤコビ行列の計算
// calculate the Duffing equations and the Jacobian matrix
solve_duffing_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++ ){
// ダフィング方程式の値とヤコビ行列の計算
// calculate the Duffing equations and the Jacobian matrix
solve_duffing_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];
}
// リアプノフ指数を大きい順に並べ替える
// sort Lyapunov exponents in ascending order
quick_sort( (int)DIM, lyap );
// リアプノフ次元の計算
// calculate the Lyapunov dimension
lyap_dim = calc_lyap_dim( lyap );
/*------------------------------------------------------------------*/
printf("#---------------------------------------------------\n" );
printf("# ITERATION: %d\n", (int)ITERATION);
printf("# SKIPS: %d\n", (int)SKIPS);
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("# B = %g, K = %g\n", (double)B, (double)K );
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", -K );
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;
}
/*----------------------------------------------------------------------------*/
void quick_sort( int n, double x[DIM] )
{
int i, j;
double pivot, t;
if (n < 2) {
return; /* no need for sorting */
}
pivot = x[0];
i = 0;
j = n-1; /* lower and upper bounds */
while (i <= j) {
while (x[i] > pivot){
i++; /* search for x[i] >= pivot */
}
while (x[j] < pivot){
j--; /* search for x[j] <= pivot */
}
if (i <= j) { /* swap x[i] and x[j] */
t = x[i];
x[i] = x[j];
x[j] = t;
i++;
j--;
}
}
quick_sort(j+1, x); /* recursive call for smaller half */
quick_sort(n-i, x+i); /* recursive call for larger half */
}