# Thread: Numeric Differentiation and Integration

1. ## Numeric Differentiation and Integration

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_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;           /* 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) {

}

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

}```