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 lyap_dim;

double sum_lyap;

double R_diag; // 行列Rの対角成分, the diagonal component of the matrix R

double v[DIM]; // 変数 variables

double lyap[DIM]; // リアプノフ指数, Lyapunov exponent 

double dfQ[DIM][DIM]; // ヤコビ行列dfと行列Qの積, Product of Jacobi matrix df and matrix Q

double Q[DIM][DIM]; // QR分解のQ, Q of QR decomposition

double R[DIM][DIM]; // QR分解のR, R of QR decomposition

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

// 初期値の設定

// set the initial values

// v[0]: x    v[1]: y    v[2]: z

v[0] = X0;

v[1] = Y0;

v[2] = 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;

}

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

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

// ローレンツ方程式の値とヤコビ行列の計算

// calculate the Lorenz equations and the Jacobian matrix

solve_lorenz_eqs_using_RK4( v, dfQ, Q );

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

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

}

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

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

solve_lorenz_eqs_using_RK4( v, dfQ, Q );

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

}

}

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

// リアプノフ指数の計算

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

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

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

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

printf("# SIGMA = %g, R = %g, B = %g\n", (double)SIGMA, (double)RR, (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;

}