Tuesday, February 12, 2013

Numberic Integation methods for calculus - Trapazoidal approximation method.

You can integrate functions without using calculus, but it would be amazingly difficult to get good results by hand.  You would have to grind through functions putting in values for each interval you were trying to calculate and then total those values up and divide by the change in x over the interval by 2.  But what we do have now are computers that will happily run through a million calculations before you can take your finger off the enter key.  Maybe you want to double check an integration you are doing to see if the results are in the ball park, a sanity check, if you will.

These functions will require you to tweak them if you go changing them around, and no guarantee is made for their accuracy.  I have not systematically tested them for anything and do not warranty them for any use other than possibly to double check a homework assignment or two. 

The first program I tied I got amazingly close to the results in my calculus book.  


.#include <stdio.h>
#include <stdlib.h>

// return x squared
double f (double x){
        return (x * x );
}

int main (){

        long double b = 1;      // interval begin
        long double e = 2;      // interval end
        long double n = 1000000;        // number of steps., denomenator
        long double r = 0;      // result

        long double c=(e-b); // The change in x per unit
        long double ca=c/n/2; // The multiplier for final result
        long double ba = b * n + c; // the begining numerator
        long double ea = e * n;     // the ending numerator

        for (;ba<ea; ba+=c)
                r += f(ba/n);  // sum all the intervals
        r *= 2;  // double all the interval values.
        r += f(b) + f(e); //handle the first and last interval
        r *= ca; // apply the final multiplier
        printf("Answer is %.20Lf \n", r);
}
Save the above as program.c 
Compile it with gcc program.c -o name 
Then run it with   ./name
 I got this as my answer:
Answer is 2.33333333333350000744


The following code is similar to the above code, but I changed the square function to return the result of a polynomial. I don't have a good way of running this as

#include <stdio.h>
#include <stdlib.h>

// return a polynomial
double f (double x){
        return (2 * x * x +4*x -5);
}

int main (){

        long double b = 1;      // interval begin
        long double e = 2;      // interval end
        long double n = 1000000;        // number of steps., denomenator
        long double r = 0;      // result

        long double c=(e-b); // The change in x per unit
        long double ca=c/n/2; // The multiplier for final result
        long double ba = b * n + c; // the begining numerator
        long double ea = e * n;     // the ending numerator

        for (;ba<ea; ba+=c)
                r += f(ba/n);  // sum all the intervals
        r *= 2;  // double all the interval values.
        r += f(b) + f(e); //handle the first and last interval 
        r *= ca; // apply the final multiplier
        printf("Answer is %.20Lf \n", r);
}

This gave me a result of   5.66666666666699991731  but I am not sure how to check if that is right.  I will integrate the function, and then calculate for the interval using those numbers to see if it gives similar results.

I  am not sure if this program will work with fractional intervals, will require additional testing to see if that works.


A more accurate method is to use Simpson's rule, Approximating using Parabolas, which I will attempt to perform in the next few days.

I updated the interval calculation to work right when the range is not an interval of 1.

No comments:

Post a Comment