0

Im currently having a problem with the program i developed to calculate the value of my integral from 0 to infinity for the function 1/((1+x^2)(x^(2/3))). my code is as follows:

#include <stdio.h>
#include <math.h>
#include <float.h>

double integral_Function(double x){ //Define a function returning a double, this funtion calculates f(x) for a given x value
    double value;
    value = 1/((1+ pow(x, 2))*(pow(x, (2/3)))); //pow(x, 2) = x^2
    return value;
}

double trapezium_Rule(double lower, double upper, int n){ //Define a function containing an algorithm for the trapezium rule

    double h;
    double total;
    int i;

    h = (upper - lower)/n; //h = strip width

    total = (h/2) * (integral_Function(upper) + integral_Function(lower)); // Add the two end points

    for(i = 1; i < n; i++){
    total += (h * integral_Function(lower + (i * h))); //add the volume of each of the middle points to the total
    }

    return total;

}

int main() {

    double sum = 0;
    double upper = DBL_EPSILON * 1.001; //DBL_EPSILON is the smallest non-zero value of a double
    double lower = DBL_EPSILON;
    int accuracy;

    printf("Enter the degree of accuracy required: \n");
    scanf("%d", &accuracy);

    while (upper <= DBL_MAX){
        sum +=  trapezium_Rule(lower, upper, 2000);
        lower = upper;
        upper = upper * 1.02;
    }

    printf("The Integral = %.*f to %d decimal places.", accuracy,sum, accuracy);

    return 0;
}

In my program the integral is split up into many small chunks using the while loop in int main(), I have specified the number of strips for each integral.

Currently my code produces an answer that is exactly half what it should be, however I cannot see where this problem occurs. It also produces an incorrect value for other functions.

Magisch
  • 7,312
  • 9
  • 36
  • 52
Sam Handley
  • 1
  • 1
  • What exactly does `h = (upper - lower)/n;` do? How do you get the height with this statement? – Marievi Mar 17 '16 at 14:01
  • h is the strip width i.e. the width of each trapezium. And n is the number of strips – Sam Handley Mar 17 '16 at 14:02
  • `pow(x, (2/3))` doesn't do what you expected – phuclv Mar 17 '16 at 15:04
  • Possible duplicate of [Division result is always zero](http://stackoverflow.com/questions/2345902/division-result-is-always-zero) – phuclv Mar 17 '16 at 15:04
  • and `x*x` is more efficient than `pow(x, 2)` – phuclv Mar 17 '16 at 15:04
  • Quick note, replace `pow(x, (2/3)` with `pow(x, (2/3.0)` In C, dividing two integers results in integer division, i.e., the fraction is chopped off. 2/3 will be just 0. If at least one of the operands is has a decimal point, the whole division becomes a real division. – Elan Mar 17 '16 at 15:06
  • Thanks very much for the help, thats solved the issue of the result being half the correct value. However this has dramatically reduced the accuracy of my result, any idea why? – Sam Handley Mar 17 '16 at 15:25
  • 1
    Why don't you try debugging by redefining integral_Function to a simpler one, like y = x, reduce the number of integration points to 5 or so and see if you can't do exact integration? – Elan Mar 17 '16 at 15:40
  • what does it mean by `dramatically reduced the accuracy of my result`? How can it reduce accuracy. What's your input output and expected value? – phuclv Mar 18 '16 at 05:20

0 Answers0