Find the area under a function in the range from a to b

Question

For an assignment, I am trying to find the area function F(X) between the range from a to b, [a,b].

Using calculus, this wouldn't be so hard. I did base this on the theorems of calculus to find the area and worked my way around it to reach certain parts of code, like this:

Note: I am using f = x**2 for testing.

def integrate(a,b,tolerance_level):
    firsttrapezoid = simpleIntegrate(a,b)
    secondtrapezoid = simpleIntegrate(a,b/2) + simpleIntegrate(b/2,b)
    error_range = abs(firsttrapezoid - secondtrapezoid)
    if error_range < tolerance_level:
        return secondtrapezoid
    else:
        return integrate(a, b/2, tolerance_level/2) + integrate(b/2, b, tolerance_level/2)

def simpleIntegrate(a,b):
    return (b-a)*(f(a)+f(b))/2

def f(x):
    f = x**2
    return f

result = integrate(0,5,0.0001)

print(result)

The problem is, I should get a value around 41, but the value I get is around 44.

Check how you divide your [a,b] intervals into two smaller intervals. Check specifically how you compute your midpoint. If you just to have a more clear view: do not start with an example from [0, 5]. Start with an example [1, 5]. — joanolo, Dec 10 '16 at 23:57
@joanolo I was using [0,5] since it is one of the input examples I was given, but I'll give it a try. — tokyolerd, Dec 11 '16 at 00:11
=> and your midpoint is just 5/2. But if you work with [1, 5]... your midpoint is ... — joanolo, Dec 11 '16 at 00:17

score 2 · Accepted Answer · answered Dec 11 '16 at 01:00

change b/2 to the midpoint between a and b which is (a+b)/2

def integrate(a, b, tolerance_level):
    firsttrapezoid = simpleIntegrate(a, b)
    secondtrapezoid = simpleIntegrate(a, (a + b) / 2) + simpleIntegrate((a + b) / 2, b)
    error_range = abs(firsttrapezoid - secondtrapezoid)
    if error_range <= tolerance_level:
        return secondtrapezoid
    else:
        return integrate(a, (a + b) / 2, tolerance_level / 2) + integrate((a + b) / 2, b, tolerance_level / 2)


def simpleIntegrate(a, b):
    return (b - a) * (f(a) + f(b)) / 2


def f(x):
    f = x ** 2
    return f


def intf(x):
    int_f = (x ** 3) / 3
    return int_f


a = 0
b = 5
tolerance = 0.0001
result = integrate(a, b, tolerance)
exactly = intf(b) - intf(a)
error = abs(exactly-result)
print("aprox: {approx} exactly: {exactly} error:{error} max error:{max_error}"
      .format(approx=result, exactly=exactly, error=error, max_error=tolerance))

Output:

aprox: 41.66668653488159 exactly: 41.666666666666664 error:1.9868214927498684e-05 max error:0.0001

John Coleman · Answer 2 · 2016-12-11T02:05:15.470

0

@eyllanesc (as well as @joanolo) have pointed out the error in the midpoint calculation.

Two other remarks:

1) Hard-wiring in f as a global name of a function is poor design. What if a person wants to integrate two or more functions as part of a single calculation? Your approach would force them to repeatedly redefine f, which would be quite inconvenient. Instead, I would recommend changing

def integrate(a, b, tolerance_level):

to

def integrate(f, a, b, tolerance_level):

with similar changes to simpleIntegrate, and appropriate adjustments made to the lines in which either integrate or simpleIntegrate are called. The resulting function will be much more flexible. Among other things, it would allow the function to be integrated to be passed as an anonymous function.

2) The algorithm that you are implementing works for most integrals but spectacularly fails for some. After making the adjustments I recommend above,

>>> def f(x): return 150*x*(1-x)*(x+1)**2
>>> integrate(f,-1,1,0.001)
0.0

But the answer should be around 40. Just because a function takes on the same value at the endpoints and at the midpoint, it doesn't follow that the function is constant, but this algorithm in effect assumes that it is. On the other hand, your algorithm does work for most functions and most intervals, so if you were told to implement it that way I wouldn't worry about it too much.

edited Dec 11 '16 at 02:05

answered Dec 11 '16 at 01:33

John Coleman

51,337
7
54
119

Thank you so much for the feedback. I was told that my function "integrate()" is supposed to be recursive to receive full credit, but I believe I'm not supposed to be using the math module for this assignment – tokyolerd Dec 11 '16 at 01:39
I just confirmed that I am not allowed to import math, so I assume we are only trying to solve simple integrals – tokyolerd Dec 11 '16 at 01:46
The import of math is neither here nor there. I could find a polynomial that behaves in much the same way. Your `integrate` has the property that it returns `simpleIntegrate` whenever the points `(a,f(a))`, `(m,f(m))`, and `(b,f(b))` are collinear (where `m = (a+b)/2`). Your approach to integration seems to assume that any such function is linear. It need not be. – John Coleman Dec 11 '16 at 01:51
@Roy740 I edited so that it is a polynomial example. You can adjust the constant in the front so that the true value of the integral is any chosen real number, but your `integrate` will always return `0`. – John Coleman Dec 11 '16 at 02:06
So in short words, my `integrate` function would not work for, let's say, a quadratic function? – tokyolerd Dec 11 '16 at 02:11
@Roy740 For a quadratic it would work. It should also work with cubics. It will fail with *some* polynomials of degree 4 or higher on *some* intervals. My example has degree 4. – John Coleman Dec 11 '16 at 02:15
Makes sense. According to my assignment, I am only supposed to find the area under a trapezoid, which is w*(h1+h2)/2, which in other words is `(b-a) * (F(a) + F(b))/2` I guess he tried to point out that we must use the function like that so I don't think I am really supposed to get too deep about what kind of functions are inputted – tokyolerd Dec 11 '16 at 02:23
@Roy740 Your Python code looks pretty good (after than bug is fixed). Numerical integration is trickier than it seems at first, but such questions are more for a course in numerical analysis than an intro to programming course (or whatever it is that you are doing). – John Coleman Dec 11 '16 at 02:28
it is an intro to programming course. One last question, is there any way I can take the antiderivative of any function that the user inputs? – tokyolerd Dec 11 '16 at 02:46
@Roy740 No easy way. By the fundamental theorem of calculus, numerical integration with a variable upper limit of integration *is* an antiderivative. There are modules for symbolic integration: http://docs.sympy.org/dev/modules/integrals/integrals.html . For the special case of polynomial functions it is easy enough to write your own function which will return the antiderivative given the list of coefficients. – John Coleman Dec 11 '16 at 02:55
I see. I'll give it a try and make another post if nothing I do seems to work – tokyolerd Dec 11 '16 at 03:06

Find the area under a function in the range from a to b

2 Answers2