# Numerical Calculus in C

1. Dec 3, 2011

### syberraith

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:

Code (Text):
# 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_vel = 1256;       /* gryoscopic velocity */
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) {

}

/* 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 */

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 );

}

2. Dec 3, 2011

### Staff: Mentor

Looks like there's a bug in your program. Increasing the number of points should give more correct results, up to a point, but it doesn't seem reasonable that the result should change that drastically.

Do you have a debugger to work with? If so, do you know how to use it?

To isolate the source of the problem, start with a small number of subintervals, so that you can manually calculate all of the values. After doing that, run your program and compare the calculated results with those produced by your program. If both sets of results agree, double the number of subintervals. Hand-calculate the results again, and compare with those produced by your program. At some point, the results you calculate should differ from those produced by the program, and that should help pinpoint which function is producing incorrect results. With a debugger, you can put in a breakpoint, and single-step through your code to find the exact line that is causing trouble.

1) Instead of hardcoding a value for pi, you can get a value that is more precise using 4.0 * atan(1.0).
2) Some of your variable and function names are nearly inscrutable, most notably d_the. A much better name would be delta_theta. Instead of h_vel and v_vel, hor_vel and vert_vel are much more descriptive, and would give a naive reader (like myself) a better idea as to what they are supposed to return.

3. Dec 3, 2011

### syberraith

I could expand the function and variable names easy enough. I see you got the general idea.

I am using gcc to compile the program, although I'm writing it in a text editor. I might still have an eclipse programing environment on this machine, or I could get some other programming environment to do debugging.

I did imagine that it was simple enough to forgo all that... Anyway, I did some basic debugging, reduced the data points and putting printf statements in here and there to print out intermediate results as they were being calculated. From that quick look, all seems to be working properly right up to the summation total. As far as I can tell, that's the only value the fails to jive.

I could omit the delta theta correction and hand calculate for a small number of data points. That correction only seems to affect the result by a factor of 2.

I did find a 32bit value to use for pi, although perhaps calculating it would be best for portability.

Last edited: Dec 3, 2011