- #1
syberraith
- 40
- 0
I'm working on a little project that involve some mechanical modeling. So, I wrote a little c program to do some numerical differentiation and integration. I'm using a slight twist on the central difference algorithm for the differentiation and a standard trapezoid algorithm for the integration.
I have a rather comical result that baffles me as to from where it originates. Basically if I increase the number of data points by a factor of ten, the variable d_cnt in my program, the overall result decreases by a factor of ten, and vice-versa. Here's the code:
I have a rather comical result that baffles me as to from where it originates. Basically if I increase the number of data points by a factor of ten, the variable d_cnt in my program, the overall result decreases by a factor of ten, and vice-versa. 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 = 100000; /* 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 = 0; /* tangential result */
double h_sum = 0; /* horizontal result */
double v_sum = 0; /* 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 slope values stored in array "m" with simple rectangular algorithm, summing the total into t_sum */
for ( i = 0 ; i < d_cnt - 1 ; i++ )
t_sum += h * ( m[i] + m[i+1] ) / 2;
}
/* 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 h_vel() */
for ( i = 0 ; i < d_cnt ; i++ ) {
m[i] = ( h_vel(t + g) - h_vel(t - g) ) / h;
t = t + h;
}
/* Integrate */
for ( i = 0 ; i < d_cnt - 1 ; i++ )
h_sum += h * ( m[i] + m[i+1] ) / 2;
}
/* 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 v_vel() */
for ( i = 0 ; i < d_cnt ; i++ ) {
m[i] = ( v_vel(t + g) - v_vel(t - g) ) / h;
t = t + h;
}
/* Integrate */
for ( i = 0 ; i < d_cnt - 1 ; i++ )
v_sum += h * ( m[i] + m[i+1] ) / 2;
}
/* 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 );
}