Why is my numerical calculus program giving me incorrect results?

In summary: Half-width for numeric differentiation and integration... */ double g = h/2; /* Integrate */ for ( i = 0 ; i < d_cnt - 1 ; i++ ) h_sum += h * ( m[i] + m[i+1] ) / 2; }}
  • #1
syberraith
42
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:

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 );
        
}
 
Technology news on Phys.org
  • #2
syberraith said:
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.
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.

Some comments.
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
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:

1. What is numerical calculus?

Numerical calculus is a branch of mathematics that deals with using numerical methods to solve problems in calculus. It involves approximating solutions to calculus problems using algorithms and computer programming.

2. Why is numerical calculus important?

Numerical calculus is important because it allows us to solve complex calculus problems that may not have exact analytical solutions. It also has many practical applications in fields such as physics, engineering, and economics.

3. Can I use C for numerical calculus?

Yes, C is a commonly used programming language for numerical calculus. It has a wide range of built-in functions and libraries that make it suitable for performing complex mathematical calculations.

4. What are the advantages of using C for numerical calculus?

One advantage of using C for numerical calculus is its speed and efficiency. C is a low-level language, which means it can directly access and manipulate computer hardware, making it faster than higher-level languages. Additionally, C has a large community of developers and a vast library of mathematical functions that can be used for numerical calculus.

5. Are there any limitations to using C for numerical calculus?

One limitation of using C for numerical calculus is that it can be more challenging to learn and use compared to other programming languages. Additionally, C does not have built-in support for advanced data structures, which can make it more challenging to implement certain numerical algorithms.

Similar threads

  • Programming and Computer Science
Replies
3
Views
1K
  • Programming and Computer Science
2
Replies
36
Views
2K
  • Programming and Computer Science
Replies
1
Views
650
  • Programming and Computer Science
Replies
2
Views
1K
  • Programming and Computer Science
Replies
4
Views
1K
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
Replies
1
Views
2K
  • Programming and Computer Science
2
Replies
50
Views
4K
  • Programming and Computer Science
Replies
8
Views
1K
  • Programming and Computer Science
Replies
4
Views
1K
Back
Top