from __future__ import division
from scipy import *
from scipy import integrate
import pylab

def f(x):
    return 4*x**3 - 6*x**2 + 1

def intf(x):
    return x**4 - 2*x**3 + x

# Plot the function
x = linspace(-1, 1.5, 200)
pylab.plot(x, f(x))
pylab.plot([-1, 1.5], [0, 0], 'k')
pylab.plot([0, 0], [-1, 1.5], 'k')
pylab.ylim(-1, 1.5)
pylab.show()

# One-dimensional integration
a = 0.5
b = 1
print 'True value:      ', intf(b) - intf(a)
print 'Numerical method:', integrate.quad(f, a, b)

# Two-dimensional integration
def circle(y, x):
    return 1

def negative_limit(x):
    return -sqrt(1-x**2)

def positive_limit(x):
    return sqrt(1-x**2)

print integrate.dblquad(circle, -1, 1, negative_limit, positive_limit)