/* * lorenz.c -- Calculation of Lorenz attractor using GSL * Copyright (C) 2009 Telford Tendys * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * * -------------------------------------------------------------------------------- * * Calculation of the Lorenz attractor. * Equations reference from Wikipedia: * http://en.wikipedia.org/w/index.php?title=Lorenz_attractor&oldid=255197908 * Solving the ODE using the GNU scientific library gsl_odeiv_ system of functions. * * Compile with: gcc lorenz.c -lgsl -lgslcblas * */ #include #include #include #include #include struct lorenz_params_S { double sigma; /* Prandtl number */ double rho; /* Rayleigh number */ double beta; }; /* * This one is the basic definition of the lorenz attractor. * No way to check the dimension is correct here so just trust that is must be 3. */ int lorenz_definition( double t, const double y[], double f[], void *params ) { struct lorenz_params_S *P; P = params; f[ 0 ] = P->sigma * ( y[ 1 ] - y[ 0 ]); f[ 1 ] = y[ 0 ] * ( P->rho - y[ 2 ]) - y[ 1 ]; f[ 2 ] = y[ 0 ] * y[ 1 ] - P->beta * y[ 2 ]; return( GSL_SUCCESS ); } /* * This is some derivatives from the above definition which assist the solver * (but in theory they could be found by numerical means, which would be less * accurate and take longer). */ int lorenz_jacobian( double t, const double y[], double *dfdy, double dfdt[], void *params) { struct lorenz_params_S *P; gsl_matrix_view V; gsl_matrix *M; P = params; /* * The system of equations does not change over time. */ dfdt[ 0 ] = 0; dfdt[ 1 ] = 0; dfdt[ 2 ] = 0; V = gsl_matrix_view_array( dfdy, 3, 3 ); M = &V.matrix; /* * First derivatives of all the equations from lorenz_definition() * against all the input variables (9 results in total) */ gsl_matrix_set( M, 0, 0, -P->sigma ); /* df[ 0 ] / dy[ 0 ] */ gsl_matrix_set( M, 0, 1, P->sigma ); /* df[ 0 ] / dy[ 1 ] */ gsl_matrix_set( M, 0, 2, 0 ); /* df[ 0 ] / dy[ 2 ] */ gsl_matrix_set( M, 1, 0, P->rho - y[ 2 ]); /* df[ 1 ] / dy[ 0 ] */ gsl_matrix_set( M, 1, 1, -1 ); /* df[ 1 ] / dy[ 1 ] */ gsl_matrix_set( M, 1, 2, -y[ 0 ]); /* df[ 1 ] / dy[ 2 ] */ gsl_matrix_set( M, 2, 0, y[ 1 ]); /* df[ 2 ] / dy[ 0 ] */ gsl_matrix_set( M, 2, 1, y[ 0 ]); /* df[ 2 ] / dy[ 1 ] */ gsl_matrix_set( M, 2, 2, -P->beta ); /* df[ 2 ] / dy[ 2 ] */ return( GSL_SUCCESS ); } #define NUMBER_DIMS 3 /* -------------------------------------------------------------------------------- */ /* * Just some high-level wrapper to make it a bit easier */ struct odeiv_wrap_S { gsl_odeiv_step *stepper; gsl_odeiv_control *controller; gsl_odeiv_evolve *evolver; gsl_odeiv_system sys; double t; double t_stop; double h; }; int odeiv_wrap_apply( struct odeiv_wrap_S *W, double y[]) { int e; double h; double t_stop; h = W->h; t_stop = W->t + h; if( t_stop > W->t_stop ) t_stop = W->t_stop; for(;;) { e = gsl_odeiv_evolve_apply( W->evolver, W->controller, W->stepper, &W->sys, &W->t, t_stop, &h, y ); if( e != GSL_SUCCESS ) return( e ); if( W->t == t_stop ) return( GSL_SUCCESS ); } return( e ); } int odeiv_wrap_free( struct odeiv_wrap_S *W ) { gsl_odeiv_evolve_free( W->evolver ); gsl_odeiv_control_free( W->controller ); gsl_odeiv_step_free( W->stepper ); W->evolver = 0; W->controller = 0; W->stepper = 0; return( GSL_SUCCESS ); } /* -------------------------------------------------------------------------------- */ /* * Collect statistics, we don't actually look at the time, just collect equal steps * in the time axis to ensure a standard set of sample points. The statistics have no * understanding of chaos, nor of correlation over a local time period. */ struct stats_S { int N; double Ex; double Ex2; }; int stats_clear( struct stats_S *S ) { S->N = 0; S->Ex = 0; S->Ex2 = 0; return( GSL_SUCCESS ); } int stats_update( struct stats_S *S, double x ) { ++S->N; S->Ex += x; S->Ex2 += ( x * x ); return( GSL_SUCCESS ); } double stats_mean( struct stats_S *S ) { return( S->Ex / S->N ); } /* * This gets the same results as the R function "sd" to * six decimal places so that's probably close enough to correct */ double stats_sd( struct stats_S *S ) { register double tmp1; tmp1 = S->Ex * S->Ex; tmp1 /= S->N; tmp1 = S->Ex2 - tmp1; tmp1 /= ( S->N - 1 ); return( sqrt( tmp1 )); } /* -------------------------------------------------------------------------------- */ struct stats_S stats[ NUMBER_DIMS ]; struct stats_S Bstats[ NUMBER_DIMS ]; /* * We handle some basic command line parameters. * Such as setting for the param variables and initial conditions as well * as time limit and statistical measurements. */ int main( int argc, char *argv[]) { double y[ NUMBER_DIMS ] = { 0, 1, 1.05 }; struct lorenz_params_S par = { 10, 28, 8/3 }; struct odeiv_wrap_S oediv = {0}; int i; int average_count = 1; char summary = 0; oediv.t_stop = 100; oediv.h = 0.005; for( i = 1; i < argc; i++ ) { char *p = argv[ i ]; if( !strcmp( p, "--summary" )) { summary = 1; continue; } if( !strncmp( p, "--sigma=", 8 )) { par.sigma = strtod( p + 8, 0 ); continue; } if( !strncmp( p, "--rho=", 6 )) { par.rho = strtod( p + 6, 0 ); continue; } if( !strncmp( p, "--beta=", 7 )) { par.beta = strtod( p + 7, 0 ); continue; } if( !strncmp( p, "--y0=", 5 )) { y[ 0 ] = strtod( p + 5, 0 ); continue; } if( !strncmp( p, "--y1=", 5 )) { y[ 1 ] = strtod( p + 5, 0 ); continue; } if( !strncmp( p, "--y2=", 5 )) { y[ 2 ] = strtod( p + 5, 0 ); continue; } if( !strncmp( p, "--stop=", 7 )) { oediv.t_stop = strtod( p + 7, 0 ); continue; } if( !strncmp( p, "--ave=", 6 )) { average_count = (int)( strtod( p + 6, 0 ) + 0.5 ); continue; } fprintf( stderr, "UNKNOWN: %s\n", p ); exit( -1 ); } /* Decide what type of solver we want to use (the "step" method) */ oediv.stepper = gsl_odeiv_step_alloc( gsl_odeiv_step_rk8pd, NUMBER_DIMS ); /* Standard stepsize controller */ oediv.controller = gsl_odeiv_control_standard_new( 1e-7, 1e-10, 1, 1 ); /* Thing that can step forward and back until a workable step size is found */ oediv.evolver = gsl_odeiv_evolve_alloc( NUMBER_DIMS ); oediv.sys.function = &lorenz_definition; oediv.sys.jacobian = &lorenz_jacobian; oediv.sys.dimension = NUMBER_DIMS; oediv.sys.params = ∥ if( average_count < 1 ) average_count = 1; for( oediv.t = 0; oediv.t < oediv.t_stop; ) { int e; e = odeiv_wrap_apply( &oediv, y ); if( e != GSL_SUCCESS ) break; for( i = 0; i < NUMBER_DIMS; i++ ) { stats_update( stats + i, y[ i ]); } if( stats[ 0 ].N >= average_count ) { if( summary ) { for( i = 0; i < NUMBER_DIMS; i++ ) { stats_update( Bstats + i, stats_mean( stats + i )); stats_clear( stats + i ); } } else { /* * Of course you must always output to 6 decimal places because we know * that if you take an old printout with only 3 decimal place accuracy * the system will produce errors. */ printf("%.06f\t%.06f\t%.06f\t%.06f\n", oediv.t, stats_mean( stats + 0 ), stats_mean( stats + 1 ), stats_mean( stats + 2 )); for( i = 0; i < NUMBER_DIMS; i++ ) { stats_clear( stats + i ); } } } } /* * When in summary mode, just output the standard deviation of each variable */ if( summary ) { printf("%d\t%.06f\t%.06f\t%.06f\n", average_count, stats_sd( Bstats + 0 ), stats_sd( Bstats + 1 ), stats_sd( Bstats + 2 )); } odeiv_wrap_free( &oediv ); return 0; }