When I implement the following code, I get wrong answers for maybe the first dozen or two runs, and then the system hits a maximum and seems to be accurate. This is a simple dampened pendulum. The code in question follows:
struct params{
double g;
double length;
double damping;
double startPosition;
int startVelocity;
};
int func(double t, const double y[], double f[], void *params){
(void) (t);
struct params *init = (struct params*)params;
// Starting velocity
f[0] = y[1];
// Free Fall
f[1] = -(init->damping)*y[1] - (init->g / init->length)*sin(y[0]);
return GSL_SUCCESS;
}
int main(){
// Paramter initialization
struct params *inputs = malloc(sizeof(struct params));
if (inputs != NULL){
inputs->g=9.8;
inputs->length=2;
inputs->damping=0.1;
inputs->startPosition = 0.2;
inputs->startVelocity = 0;
}
int i;
double t = 0.0, t1 = 100;
double dt = 1;
const int step_total = 101;
double* timeX = malloc(step_total*sizeof(double));
double* yOne = malloc(step_total*sizeof(double));
double* yTwo = malloc(step_total*sizeof(double));
clock_t start_time = clock(); //Record Start Time
double timeout = 5.0; // Timeout in seconds
int SCREEN_WIDTH = 1280;
int SCREEN_HEIGHT = 720;
// System initialization
gsl_odeiv2_system system = {func, NULL, 2, inputs};
// ODE solver initialization
const gsl_odeiv2_step_type *type = gsl_odeiv2_step_rk8pd;
gsl_odeiv2_step *step = gsl_odeiv2_step_alloc(type, 2);
gsl_odeiv2_control *control = gsl_odeiv2_control_y_new(1e-6,0.0);
gsl_odeiv2_evolve *evolve = gsl_odeiv2_evolve_alloc(2);
double y[2] = {inputs->startPosition, inputs->startVelocity};
//printf ( " Step method is '%s'\n", gsl_odeiv2_step_name(step));
i = 0;
while(t < t1){
int status = gsl_odeiv2_evolve_apply(evolve, control, step, &system, &t, t1, &dt, y);
if (status != GSL_SUCCESS){
printf("error, return value=%d\n", status);
break;
}
//if (y[0] <= 0) break;
//printf("At time t=%.3f s, position=%.3f m, velocity=%.3f m/s\n", t, y[0], y[1]);
timeX[i] = t;
yOne[i] = y[0];
yTwo[i] = y[1];
i++;
clock_t current_time = clock();
double elapsed_time = (double)(current_time - start_time) / CLOCKS_PER_SEC;
if (elapsed_time > timeout){
printf(" Infinite loop timeout reached. \n");
break;
}
}
int total = i-1;
i = 0;
gsl_odeiv2_evolve_reset(evolve);
gsl_odeiv2_control_free(control);
gsl_odeiv2_step_free(step);
gsl_odeiv2_evolve_free(evolve);
When I output the results in an array this is what I see: If you look at the first dozen or so output lines you can tell that my initial values are not correct for position and velocity. Where does the system find the initial position of 73.609? Starting velocity is close to zero so that may be accurate. Once the position hits 100. Then the results look accurate and can be plotted.
At time t=0.575 s, position=73.609 m, velocity=0.004 m/s
At time t=1.185 s, position=74.530 m, velocity=0.000 m/s
At time t=1.721 s, position=75.452 m, velocity=-0.004 m/s
At time t=2.257 s, position=76.375 m, velocity=0.003 m/s
At time t=2.792 s, position=77.298 m, velocity=0.001 m/s
At time t=3.326 s, position=78.221 m, velocity=-0.004 m/s
At time t=3.910 s, position=79.166 m, velocity=0.003 m/s
At time t=4.499 s, position=80.111 m, velocity=0.001 m/s
At time t=5.088 s, position=81.056 m, velocity=-0.003 m/s
At time t=5.637 s, position=82.017 m, velocity=0.002 m/s
At time t=6.185 s, position=82.978 m, velocity=0.001 m/s
At time t=6.772 s, position=83.939 m, velocity=-0.003 m/s
At time t=7.371 s, position=84.919 m, velocity=0.002 m/s
At time t=7.970 s, position=85.899 m, velocity=0.000 m/s
At time t=8.525 s, position=86.879 m, velocity=-0.002 m/s
At time t=9.080 s, position=87.884 m, velocity=0.002 m/s
At time t=9.686 s, position=88.889 m, velocity=-0.001 m/s
At time t=10.296 s, position=89.894 m, velocity=-0.001 m/s
At time t=10.905 s, position=90.899 m, velocity=0.002 m/s
At time t=11.532 s, position=91.924 m, velocity=-0.001 m/s
At time t=12.158 s, position=92.949 m, velocity=-0.000 m/s
At time t=12.803 s, position=93.974 m, velocity=0.002 m/s
At time t=13.449 s, position=95.020 m, velocity=-0.002 m/s
At time t=14.094 s, position=96.067 m, velocity=0.001 m/s
At time t=14.739 s, position=97.114 m, velocity=0.001 m/s
At time t=15.384 s, position=98.161 m, velocity=-0.001 m/s
At time t=16.029 s, position=99.233 m, velocity=0.001 m/s
At time t=16.675 s, position=100.000 m, velocity=0.000 m/s
At time t=17.320 s, position=0.071 m, velocity=-0.103 m/s
At time t=17.965 s, position=-0.033 m, velocity=-0.163 m/s
At time t=18.618 s, position=-0.075 m, velocity=0.054 m/s
At time t=19.271 s, position=0.013 m, velocity=0.166 m/s
At time t=19.937 s, position=0.074 m, velocity=-0.015 m/s
At time t=20.602 s, position=0.002 m, velocity=-0.158 m/s
At time t=21.268 s, position=-0.068 m, velocity=-0.016 m/s
At time t=21.934 s, position=-0.015 m, velocity=0.145 m/s
At time t=22.600 s, position=0.061 m, velocity=0.042 m/s
At time t=23.266 s, position=0.026 m, velocity=-0.127 m/s
At time t=23.931 s, position=-0.052 m, velocity=-0.064 m/s
At time t=24.603 s, position=-0.033 m, velocity=0.108 m/s
At time t=25.275 s, position=0.044 m, velocity=0.077 m/s
At time t=25.950 s, position=0.038 m, velocity=-0.089 m/s
At time t=26.625 s, position=-0.035 m, velocity=-0.085 m/s
At time t=27.305 s, position=-0.040 m, velocity=0.071 m/s
At time t=27.986 s, position=0.028 m, velocity=0.089 m/s
At time t=28.673 s, position=0.041 m, velocity=-0.057 m/s
At time t=29.361 s, position=-0.022 m, velocity=-0.088 m/s
At time t=30.056 s, position=-0.040 m, velocity=0.046 m/s
At time t=30.751 s, position=0.018 m, velocity=0.085 m/s
At time t=31.455 s, position=0.038 m, velocity=-0.039 m/s
At time t=32.158 s, position=-0.016 m, velocity=-0.081 m/s
At time t=32.870 s, position=-0.035 m, velocity=0.036 m/s
At time t=33.581 s, position=0.015 m, velocity=0.075 m/s
At time t=34.299 s, position=0.033 m, velocity=-0.035 m/s
At time t=35.017 s, position=-0.015 m, velocity=-0.068 m/s
At time t=35.741 s, position=-0.030 m, velocity=0.036 m/s
At time t=36.465 s, position=0.016 m, velocity=0.062 m/s
At time t=37.193 s, position=0.027 m, velocity=-0.037 m/s
At time t=37.922 s, position=-0.017 m, velocity=-0.054 m/s
At time t=38.655 s, position=-0.023 m, velocity=0.040 m/s
At time t=39.388 s, position=0.018 m, velocity=0.046 m/s
At time t=40.126 s, position=0.020 m, velocity=-0.042 m/s
At time t=40.864 s, position=-0.019 m, velocity=-0.038 m/s
At time t=41.607 s, position=-0.016 m, velocity=0.044 m/s
At time t=42.349 s, position=0.020 m, velocity=0.029 m/s
At time t=43.099 s, position=0.011 m, velocity=-0.045 m/s
At time t=43.848 s, position=-0.020 m, velocity=-0.020 m/s
At time t=44.606 s, position=-0.007 m, velocity=0.045 m/s
At time t=45.364 s, position=0.020 m, velocity=0.009 m/s
At time t=46.135 s, position=0.002 m, velocity=-0.044 m/s
At time t=46.906 s, position=-0.019 m, velocity=0.003 m/s
At time t=47.676 s, position=0.003 m, velocity=0.040 m/s
At time t=48.447 s, position=0.017 m, velocity=-0.013 m/s
At time t=49.217 s, position=-0.008 m, velocity=-0.034 m/s
At time t=49.993 s, position=-0.014 m, velocity=0.021 m/s
At time t=50.768 s, position=0.011 m, velocity=0.025 m/s
At time t=51.561 s, position=0.009 m, velocity=-0.027 m/s
At time t=52.353 s, position=-0.013 m, velocity=-0.014 m/s
At time t=53.151 s, position=-0.004 m, velocity=0.030 m/s
At time t=53.949 s, position=0.013 m, velocity=0.001 m/s
At time t=54.760 s, position=-0.002 m, velocity=-0.028 m/s
At time t=55.570 s, position=-0.012 m, velocity=0.011 m/s
At time t=56.381 s, position=0.007 m, velocity=0.021 m/s
At time t=57.194 s, position=0.008 m, velocity=-0.019 m/s
At time t=58.008 s, position=-0.010 m, velocity=-0.011 m/s
At time t=58.834 s, position=-0.003 m, velocity=0.023 m/s
At time t=59.661 s, position=0.010 m, velocity=-0.001 m/s
At time t=60.504 s, position=-0.003 m, velocity=-0.020 m/s
At time t=61.347 s, position=-0.008 m, velocity=0.012 m/s
At time t=62.190 s, position=0.007 m, velocity=0.012 m/s
At time t=63.045 s, position=0.003 m, velocity=-0.018 m/s
At time t=63.901 s, position=-0.008 m, velocity=-0.000 m/s
At time t=64.763 s, position=0.002 m, velocity=0.016 m/s
At time t=65.626 s, position=0.006 m, velocity=-0.010 m/s
At time t=66.488 s, position=-0.006 m, velocity=-0.009 m/s
At time t=67.366 s, position=-0.002 m, velocity=0.015 m/s
At time t=68.244 s, position=0.007 m, velocity=-0.003 m/s
At time t=69.134 s, position=-0.003 m, velocity=-0.012 m/s