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.(adsbygoogle = window.adsbygoogle || []).push({});

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

}

**Physics Forums | Science Articles, Homework Help, Discussion**

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Numerical Calculus in C

Loading...

Similar Threads for Numerical Calculus |
---|

C/++/# How to use C++ in studying calculus |

Python Numerical integration 'quad' error |

Numerical Integration with variable limits MATLAB |

Python Integrating to Infinity Numerically |

C/++/# How to integrate when one of the limits is a variable? |

**Physics Forums | Science Articles, Homework Help, Discussion**