/*-----------------------------------------------------------------
マッキー・グラス遅延微分方程式
Mackey-Glass delay differential equation
dx/dt = A * x(t-tau)/( 1 + x(t-tau)^C ) - B * x(t)
mackey-glass_DDE_lyap.c
-----------------------------------------------------------------*/
#include "mackey-glass_DDE_lyap.h"
double x[DIM];
double u[DIM][DIM];
double e[DIM][DIM];
double lyap[DIM];
int main(void)
{
int i;
int j;
int t;
int kmax;
double d2;
double dxdt;
double total_T;
double lyap_dim;
double lyap_sum;
/*---------------------------------------------------------------*/
// 初期値の設定
// set the initial values
for ( i = 0; i < DIM; i++ ) {
x[i] = 0.9;
lyap[i] = 0.0;
/* 0 < x < 1 */
// x[i] = rand()/(double)RAND_MAX;
}
/*---------------------------------------------------------------*/
// 直交化した長さ１のベクトル（単位行列）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 );
}
}
/*---------------------------------------------------------------*/
for ( t = 0; t < SKIP; t++ ) {
// マッキー・グラス遅延微分方程式の値とヤコビ行列の計算
// calculate the Mackey-Glass delay differential equation and the Jacobian matrix
solve_mackey_glass_DDE_using_RK4( x, &dxdt, u );
// printf("%lf %lf %lf\n", x[0], dxdt, x[X_TAU] ); /* 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 Mackey-Glass delay differential equation and the Jacobian matrix
solve_mackey_glass_DDE_using_RK4( x, &dxdt, u );
// printf("%lf %lf %lf\n", x[0], dxdt, x[X_TAU] ); /* output the data */
// グラム・シュミットの正規直交化法
// 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
total_T = ITERATION * DELTA_T;
for ( i = 0; i < DIM; i++ ) {
lyap[i] = lyap[i]/(double)total_T;
}
// リアプノフ指数を大きい順に並べ替える
// sort Lyapunov exponents in ascending order
quick_sort( (int)DIM, lyap );
// リアプノフ次元の計算
// calculate the Lyapunov dimension
calc_lyap_dim( lyap, &lyap_dim, &kmax );
lyap_sum = 0.0;
for ( i = 0; i < DIM; i++ ){
lyap_sum += lyap[i];
}
printf("#---------------------------------------------------\n" );
printf("# ITERATION: %d\n", (int)ITERATION);
printf("# SKIP: %d\n", (int)SKIP);
printf("# DELTA_T: %g\n", (double)DELTA_T);
printf("# DIM: %d\n", (int)DIM);
printf("# TOTAL_TIME: %g\n", (double)total_T);
printf("#---------------------------------------------------\n" );
printf("# TAU = %g, A = %g, B = %g, C = %g\n", (double)TAU, (double)A, (double)B, (double)C);
printf("#---------------------------------------------------\n" );
printf("# sum of lyap: %lf\n", lyap_sum);
printf("# Lyapunov Dimension: %lf\n", lyap_dim);
printf("#----------------------------------\n");
printf("# Lyapunov Spectrum (log with base-e)\n" );
for ( i = 0; i < 5; i++ ){ /* the first five exponents */
// for ( i = 0; i <= kmax; i++ ){ /* necessary exponents to calculate Lyapunov dimension */
// for ( i = 0; i < DIM; i++ ){ /* all exponents */
printf("%lf\n", lyap[i]);
}
return 0;
}
/*-----------------------------------------------------------------------------------------*/
double log2( double x )
{
return log(x)/(double)log(2.0);
}
/*-----------------------------------------------------------------------------------------*/
// リアプノフ次元の計算
// calculate Lyapunov dimension
void calc_lyap_dim( double lyap[DIM], double (*lyap_dim), int (*iii) )
{
int i;
double sum;
(*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] );
(*iii) = i;
break;
}
}
}
/*-----------------------------------------------------------------------------------------*/
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 */
}