4

I have a range of data that I have approximated using a polynomial of degree 2 in Python. I want to calculate the area underneath this polynomial between 0 and 1.

Is there a calculus, or similar package from numpy that I can use, or should I just make a simple function to integrate these functions?

I'm a little unclear what the best approach for defining mathematical functions is.

Thanks.

Mike Graham
  • 73,987
  • 14
  • 101
  • 130
djq
  • 14,810
  • 45
  • 122
  • 157
  • 2
    If it's a polynomial of degree 2, just integrate it by hand - there's no need to use code. –  Feb 28 '10 at 21:15
  • It's for a large set of polynomials, that I process in batches of 40. – djq Feb 28 '10 at 21:22
  • The area under the curve of a polynomial of degree 2 is a polynomial with a degree of 1. Just plug values into this equation. – S.Lott Feb 28 '10 at 21:52
  • @S.Lott It's degree 3, but yes. –  Feb 28 '10 at 22:08
  • @Paul Hankin: Ooops -- for some reason I was thinking of the derivative of a polynomial at each point. – S.Lott Mar 01 '10 at 00:41

5 Answers5

9

If you're integrating only polynomials, you don't need to represent a general mathematical function, use numpy.poly1d, which has an integ method for integration.

>>> import numpy
>>> p = numpy.poly1d([2, 4, 6])
>>> print p
   2
2 x + 4 x + 6
>>> i = p.integ()
>>> i
poly1d([ 0.66666667,  2.        ,  6.        ,  0.        ])
>>> integrand = i(1) - i(0) # Use call notation to evaluate a poly1d
>>> integrand
8.6666666666666661

For integrating arbitrary numerical functions, you would use scipy.integrate with normal Python functions for functions. For integrating functions analytically, you would use sympy. It doesn't sound like you want either in this case, especially not the latter.

Mike Graham
  • 73,987
  • 14
  • 101
  • 130
  • Great - thanks! That is very useful. So for calculating the area from 0 to 1, I can use: Area = i(1) - i(0) ? – djq Feb 28 '10 at 20:38
  • 1
    That will be the definite integral from 0 to 1. In this case, that is the same as the area, but in some situations (where part or all of the polynomial is negative), it is not. – Mike Graham Feb 28 '10 at 20:40
5

Look, Ma, no imports!

>>> coeffs = [2., 4., 6.]
>>> sum(coeff / (i+1) for i, coeff in enumerate(reversed(coeffs)))
8.6666666666666661
>>>

Our guarantee: Works for a polynomial of any positive degree or your money back!

Update from our research lab: Guarantee extended; s/positive/non-negative/ :-)

Update Here's the industrial-strength version that is robust in the face of stray ints in the coefficients without having a function call in the loop, and uses neither enumerate() nor reversed() in the setup:

>>> icoeffs = [2, 4, 6]
>>> tot = 0.0
>>> divisor = float(len(icoeffs))
>>> for coeff in icoeffs:
...     tot += coeff / divisor
...     divisor -= 1.0
...
>>> tot
8.6666666666666661
>>>
John Machin
  • 81,303
  • 11
  • 141
  • 189
  • +1, Good solution (practically the same as the one `integ` implements, I'm sure). I would make my denominator `float(i + 1)` if I wasn't in a file with `from __future__ import division`. – Mike Graham Feb 28 '10 at 23:22
3

It might be overkill to resort to general-purpose numeric integration algorithms for your special case...if you work out the algebra, there's a simple expression that gives you the area.

You have a polynomial of degree 2: f(x) = ax2 + bx + c

You want to find the area under the curve for x in the range [0,1].

The antiderivative F(x) = ax3/3 + bx2/2 + cx + C

The area under the curve from 0 to 1 is: F(1) - F(0) = a/3 + b/2 + c

So if you're only calculating the area for the interval [0,1], you might consider using this simple expression rather than resorting to the general-purpose methods.

Jim Lewis
  • 43,505
  • 7
  • 82
  • 96
1

'quad' in scipy.integrate is the general purpose method for integrating functions of a single variable over a definite interval. In a simple case (such as the one described in your question) you pass in your function and the lower and upper limits, respectively. 'quad' returns a tuple comprised of the integral result and an upper bound on the error term.

from scipy import integrate as TG

fnx = lambda x: 3*x**2 + 9*x    # some polynomial of degree two
aoc, err = TG.quad(fnx, 0, 1)

[Note: after i posted this i an answer posted before mine, and which represents polynomials using 'poly1d' in Numpy. My scriptlet just above can also accept a polynomial in this form:

import numpy as NP

px = NP.poly1d([2,4,6])
aoc, err = TG.quad(px, 0, 1)
# returns (8.6666666666666661, 9.6219328800846896e-14)
doug
  • 69,080
  • 24
  • 165
  • 199
  • Since polynomials can be analytically integrated trivially, it is best to use the more-efficient `integ` method rather than repeated Gauss quadrature for them. – Mike Graham Feb 28 '10 at 21:34
1

If one is integrating quadratic or cubic polynomials from the get-go, an alternative to deriving the explicit integral expressions is to use Simpson's rule; it is a deep fact that this method exactly integrates polynomials of degree 3 and lower.

To borrow Mike Graham's example (I haven't used Python in a while; apologies if the code looks wonky):

>>> import numpy
>>> p = numpy.poly1d([2, 4, 6])
>>> print p
   2
2 x + 4 x + 6
>>> integrand = (1 - 0)(p(0) + 4*p((0 + 1)/2) + p(1))/6

uses Simpson's rule to compute the value of integrand. You can verify for yourself that the method works as advertised.

Of course, I did not simplify the expression for integrand to indicate that the 0 and 1 can be replaced with arbitrary values u and v, and the code will still work for finding the integral of the function from u to v.