/***************************************************************************************************************** This program was written by Kostas Chlouverakis ----------------------------------------------- It calculates the Lyapunov exponents spectrum for a DDE. It is done using the Euler integration, following Farmer's work and ideas: Farmer, J. D. "Chaotic Attractors of an Infinite-Dimensional Dynamical System." Physica D 4 (1982): 366-393. Variable N should be as high as possible. Try 100-400. If N>100 then it will need many-many hours... The integration step should be equal to: h=tau/N. Time-step h should have usually a maximum value at 0.1-0.2. So, if tau is high then N should be high as well in order the calculation to be accurate. Discussions with Prof. J.C. Sprott, Dr. Colet and Dr. P. Grassberger are greatly acknowledged!!!!!!!! *******************************************************************************************************************/ #include #include #include #include #include #include #include #include #define pi 3.1415926535 #define value_zero 0.005 /* If an LE is really zero, numerically it will be: LE=|value_zero| and never LE=0.00...0 It is up to you to set it as low as you want to. */ #define N 100 #define Na (N+1)*(N+1)+N+2 /* Don't mess with this!! */ #define NN (N+1)*(N+1)+N /* Don't mess with this too!! */ double tau = 5. ; double b = 5.; double fi= -pi/4.; double imax=pow(10,5); int num=0; double X[Na], Xnew[Na] ; double CUM[N+1], Znorm[N+1], GSC[N+1], lyp_exp[N+1] ; void main(void) { double h=tau/N, t=0; int i, k, j, l; ofstream fout; ofstream fout1; ofstream fout2; ofstream fout3; fout.open("xOF.txt"); fout2.open("yOF.txt"); fout1.open("tOF.txt"); fout3.open("LEs.txt"); /* Init. Condts */ for (int mal=1; mal<=N; mal++) X[mal]=1. ; /* Init. Condts for linear system (orthonormal frame) */ for (i = N+1; i <= NN; i++) X[i] = 0.0; for (i = 1; i <= N; i++) { X[(N+1)*i] = 1.0; CUM[i] = 0.0; } /* Start iterating... */ for (int d=1; d=0 ) { kmax=j1; lsum=sum; } } Dky = kmax - lsum/lyp_exp[kmax+1]; /* The value of the Jacobian trace, that is equal to the sum of all LEs */ for (int j2=1; j2<=N; j2++) LEsum = LEsum + lyp_exp[j2] ; /* Check positive, zero and negative LEs */ int pos=0 ,zero=0, negative=0, jl; for (jl=1; jl<=N; jl++) { if(lyp_exp[jl]>value_zero) pos=pos+1; if(fabs(lyp_exp[jl])