Numeric Differentiation and Integration

This is a discussion on Numeric Differentiation and Integration within the C Programming forums, part of the General Programming Boards category; I'm working on a little program to do some mechanical modeling involving numeric differentiation and integration. I implemented my own ...

  1. #1
    Registered User
    Join Date
    Dec 2011
    Posts
    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_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;           /* 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) { 
    
       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 );
            
    }
    Last edited by syberraith; 12-03-2011 at 02:25 AM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. numerical differentiation of data
    By lordbubonicus in forum C Programming
    Replies: 3
    Last Post: 10-25-2007, 04:01 PM
  2. perform numeric integration
    By Frank_Rye in forum C Programming
    Replies: 5
    Last Post: 10-18-2005, 02:17 PM
  3. Integration
    By westclok in forum C Programming
    Replies: 8
    Last Post: 02-01-2005, 05:06 AM
  4. Expression Manipulator v0.3 (differentiation)
    By ygfperson in forum A Brief History of Cprogramming.com
    Replies: 5
    Last Post: 06-01-2003, 10:15 PM
  5. differentiation?
    By HAcKADEmY in forum C Programming
    Replies: 2
    Last Post: 04-15-2002, 02:11 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21