double_pendulum_lyap.c

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

2重振り子の方程式のリアプノフスペクトラムの計算

calculating the Lyapunov spectrum of the double pendulum equations

double_pendulum_lyap.c

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

#include "double_pendulum_lyap.h"

int main(void)

{

int i;

int j;

int t;

double x1; // おもりM1のx座標, x-coordinate of the weight M1

double y1; // おもりM1のy座標, y-coordinate of the weight M1

double x2; // おもりM2のx座標, x-coordinate of the weight M2

double y2; // おもりM2のy座標, y-coordinate of the weight M2

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]: theta1    v[1]: theta2    v[2]: omega1, \dot{\theta1}    v[3]: omega2, \dot{\theta2}

v[0] = (THETA1 * M_PI)/(double)180.0; // 角度からラジアンへの変換, from degree to radian

v[1] = (THETA2 * M_PI)/(double)180.0; // 角度からラジアンへの変換, from degree to radian

v[2] = OMEGA1;

v[3] = OMEGA2;

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

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

}

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

// 初期状態

// initial condition

// x1 =  L1 * sin( v[0] );

// y1 = -L1 * cos( v[0] );

// printf("%.10lf   %.10lf   %.10lf   %.10lf\n", 0.0, 0.0, x1, y1 );

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

// 2重振り子の方程式の値とヤコビ行列の計算

// calculate the double pendulum equations and the Jacobian matrix

solve_double_pendulum_eqs_using_RK4( v, u );

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

/*

x1 =  L1 * sin( v[0] );

y1 = -L1 * cos( v[0] );

x2 =  x1 + L2 * sin( v[1] );

y2 =  y1 - L2 * cos( v[1] );

printf("%.10lf   %.10lf   %.10lf   %.10lf\n", x1, y1, x2, y2 );

*/

}

// return 0;

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

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

// 2重振り子の方程式の値とヤコビ行列の計算

// calculate the double pendulum equations and the Jacobian matrix

solve_double_pendulum_eqs_using_RK4( v, u );

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

/*

x1 =  L1 * sin( v[0] );

y1 = -L1 * cos( v[0] );

x2 =  x1 + L2 * sin( v[1] );

y2 =  y1 - L2 * cos( v[1] );

printf("%.10lf   %.10lf   %.10lf   %.10lf\n", x1, y1, x2, y2 );

*/

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

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

printf("#        theta2: %#g\n", (double)THETA2);

printf("#        omega1: %#g\n", (double)OMEGA1);

printf("#        omega2: %#g\n", (double)OMEGA2);

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

printf("# G = %g, L1 = %g, L2 = %g, M1 = %g, M2 = %g\n", (double)G, (double)L1, (double)L2, (double)M1, (double)M2);

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

printf("#  As there is no friction and damping in the model used, this system\n" );

printf("#  is conservative.\n" );

printf("#  Hence, the sum of all Lyapunov exponents is theoretically zero.\n" );

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 lyap_tmp[DIM+1];

double lyap_dim = 0.0;

lyap_tmp[0] = 0.0;

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

lyap_tmp[i+1] = lyap_tmp[i] + lyap[i];

if( lyap_tmp[i] > 0 && lyap_tmp[i+1] < 0 ){

lyap_dim = i + lyap_tmp[i]/(double)fabs( lyap[i] );

}

}

return lyap_dim;

}