/* euler.c This code performs a euler integration of the ordinary diff. equ.: d2x/dt2 + x = 0 It is mainly useful for seeing how unstable this method is. J. Watlington, 10/16/95 Usage: euler [ ] */ #include #include #include #define DEFAULT_FINAL_TIME ( 100.0 * 3.1415926535 ) #define NUM_STATE_VARIABLES 2 #define DEFAULT_FINAL_ERROR 0.01 #define DEFAULT_STEP_SIZE 0.001 #define STEP_SIZE_SCALING_FACTOR 0.40 #define MAX_STEP_SIZE 10.0 /* just some sanity checks */ #define MIN_STEP_SIZE 0.000001 #define MONITOR_START_TIME 0.0 #define MONITOR_END_TIME 10.0 float state[ NUM_STATE_VARIABLES ]; float initial_state[ NUM_STATE_VARIABLES ] = { 1.0, 0.0 }; double run_simulation( double step, double final_time ); void main( int argc, char **argv ) { double step_size = DEFAULT_STEP_SIZE; double final_time = DEFAULT_FINAL_TIME; float final_res = DEFAULT_FINAL_ERROR; float simulation_error = 1000.0; double actual_time; int i; if( argc > 2 ) { printf( "usage: %s []\n", argv[0] ); exit( -1 ); } else if( argc > 1 ) { sscanf( argv[1], "%lf", &step_size ); if( (step_size < MIN_STEP_SIZE) || (step_size > MAX_STEP_SIZE) ) { printf("%s: step_size (%f) out of range !\n", argv[0], step_size ); exit( -1 ); } } printf( "running till final time %f\n", final_time ); printf( "running with tolerance %f\n", final_res ); while( final_res < simulation_error ) { printf( "using step size %f\n", step_size ); printf( "time > x dx/dt\n" ); /* problem dependent */ actual_time = run_simulation( step_size, final_time ); /* Print final result */ printf( "%3.5f> ", actual_time ); for( i = 0; i < NUM_STATE_VARIABLES; i++) printf( "%2.5f ", state[ i ] ); printf( "\n\n correct output: %f\n", cos( actual_time ) ); simulation_error = fabs( state[ 0 ] - cos( actual_time )); step_size *= STEP_SIZE_SCALING_FACTOR; } } double run_simulation( double step, double final_time ) { int i; register double step_size = step; double t = 0.0; float new_state[ NUM_STATE_VARIABLES ]; for( i = 0; i < NUM_STATE_VARIABLES; i++ ) state[ i ] = initial_state[ i ]; while( t <= final_time ) { /* The following, along with the initial conditions specified above, determine the diff. equation set being solved... */ new_state[ 0 ] = state[ 0 ] + (step_size * state[ 1 ]); new_state[ 1 ] = state[ 1 ] - (step_size * state[ 0 ]); /* Now update the "real" state variables with our new ones */ for( i = 0; i < NUM_STATE_VARIABLES; i++ ) state[ i ] = new_state[ i ]; #ifdef MONITOR_SIMULATION /* Print out information about the simulation */ if( (t >= MONITOR_START_TIME) && (t <= MONITOR_END_TIME) ) { printf( "%3.5f> ", t ); for( i = 0; i < NUM_STATE_VARIABLES; i++) printf( "%2.4f ", state[ i ] ); printf( "\n" ); } #endif MONITOR_SIMULATION /* Increment the time index by the step size */ t += step_size; } return( t - step_size ); }