I'm working on a little program to do some mechanical modeling involving numeric differentiation and integration. I implemented my own version of a central difference algorithm and a simple rectangular integration algorithm.
The rather comical problem I ran into is that when I increase the number of data points, the variable d_cnt in my code, by a factor of ten, the final result of the program decreases by a factor of ten.
Here's the code:
Code:# include <stdio.h> # include <math.h> void process_tangential(void); void process_horizontal(void); void process_vertical(void); void process_results(void); double d_the(double); double f_vel(double); double t_vel(double); double h_vel(double); double v_vel(double); main () { process_tangential(); process_horizontal(); process_vertical(); process_results(); } const long d_cnt = 10000; /* division_count */ const double g_rad = 2.50e-2; /* gryoscopic radius */ const double g_vel = 1256; /* gryoscopic velocity */ const double x_rad = 1.50e-1; /* cross-arm radius */ const double x_vel = 31.4; /* cross-arm velocity */ const double g_den = 8.33e3; /* gyroscopic density */ const double g_are = 8.06e-5; /* gyroscopic area */ const double pi = 3.1415926; double t_sum; /* tangential result */ double h_sum; /* horizontal result */ double v_sum; /* vertical result */ /* process tangential data */ void process_tangential() { long i; /* Iteration counter */ double t = 0; /* Interval position */ /* Width for numeric differentiation and integration */ double h = ( (2 * pi) / g_vel ) / (double) d_cnt; /* Half-width for numeric differentiation and integration */ double g = h/2; double m[ d_cnt ]; /* array for derivative values */ /* First, numerically differentiate t_vel with central difference algorithm and put derivative values into array "m" */ for ( i = 0 ; i < d_cnt ; i++ ) { m[i] = ( t_vel(t + g) - t_vel(t - g) ) / h; t = t + h; } /* Second, numerically integrate the derivative values stored in the array with a simple rectangular algorithm, summing the total into t_sum */ t_sum = m[0] * g; for ( i = 1 ; i < d_cnt - 1 ; i++ ) { t_sum = t_sum + ( m[i] * h ); } t_sum = t_sum + ( m[i] * g ); } /* process horizontal data */ void process_horizontal() { /* Similar to process_tangential() */ long i; /* Iteration counter */ double t = 0; /* Interval position */ /* Width for numeric differentiation and integration */ double h = ( (2 * pi) / g_vel ) / (double) d_cnt; /* Half-width for numeric differentiation and integration */ double g = h/2; double m[ d_cnt ]; /* array for derivative values */ /* Differentiate */ for ( i = 0 ; i < d_cnt ; i++ ) { m[i] = ( h_vel(t + g) - h_vel(t - g) ) / h; t = t + h; } /* Integrate */ h_sum = m[0] * g; for ( i = 1 ; i < d_cnt - 1 ; i++ ) { h_sum = h_sum + ( m[i] * h ); } h_sum = h_sum + ( m[i] * g); } /* process vertical data */ void process_vertical() { /* Similar to process_tangential() */ long i; /* Iteration counter */ double t = 0; /* Interval position */ /* Width for numeric differentiation and integration */ double h = ( (2 * pi) / g_vel ) / (double) d_cnt; /* Half-width for numeric differentiation and integration */ double g = h/2; double m[ d_cnt ]; /* array for derivative values */ /* Differentiate */ for ( i = 0 ; i < d_cnt ; i++ ) { m[i] = ( v_vel(t + g) - v_vel(t - g) ) / h; t = t + h; } /* Integrate */ v_sum = m[0] * g; for ( i = 1 ; i < d_cnt - 1 ; i++ ) { v_sum = v_sum + ( m[i] * h ); } v_sum = v_sum + ( m[i] * g); } /* Calculate reference frame velocity */ double f_vel(double x) { return sqrt( ( pow((-g_rad*g_vel*sin(g_vel*x) + x_rad*x_vel),2) + pow((g_rad*g_vel*cos(g_vel*x)),2) ) ); } /* Calculate detla theta angle */ double d_the(double x) { double f_ang; /* frame tangent angle */ double r_ang; /* rotor tangent angle */ double y; /* return value */ f_ang = atan2( (-g_rad*g_vel*sin(g_vel*x) + x_rad*x_vel), (g_rad*g_vel*cos(g_vel*x)) ); r_ang = atan2( sin(g_vel*x), cos(g_vel*x) ); y = f_ang + r_ang; if ( y > pi ) y = y - ( 2 * pi ); return y; } /* Calculate tangential velocity */ double t_vel(double x) { return cos(d_the(x)) * f_vel(x); } /* Calculate horizontal velocity */ double h_vel(double x) { return -sin(g_vel*x) * cos(d_the(x)) * f_vel(x); } /* Calculate vertical velocity */ double v_vel(double x) { return cos(g_vel*x) * cos(d_the(x)) * f_vel(x); } /* print out the resutlts */ void process_results() { printf( "The tangential sum is: %1.6e\n", t_sum ); printf( "The horizontal sum is: %1.6e\n", h_sum ); printf( "The vertical sum is: %1.6e\n", v_sum ); }
