Adaptive Trapezoidal integration in Python

I’ve decided to start blogging about my experiences in graduate school and some of the more interesting projects I’ve been working on. As my first post I’ll cover the idea of trapezoidal integration and the Python code I wrote to do it.

The fundamental idea of trapezoidal integration idea is to use a sequence of trapezoids to perform integration numerically. Numerical integration is also called numerical quadrature.


from math import *

f1 = lambda x: x**3-x**2+3*x-cos(x)*x+atan(sin(x)+1)

First we import the math library, nothing special about that.

The next line may be unfamiliar for people who have not explored Python. Python allows so called “Lambda Functions”. It’s a really cool feature that is great for declaring mathematical functions (as opposed to functions that perform tasks that are not mathematical)

The above code is equivalent to saying f(x) = x^3-x^2+3x-cos(x)*x+arctan(sin(x)+1)


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

def diff2(f,x,h=1E-6):
    r=(f(x-h)-2*f(x) + f(x+h))/float(h*h)
    return r
    
def trapezint(f,a,b,n):
    h=(b-a)/n
    sum=0    
    part1=(0.5)*h*(f(a)+f(b))
    for i in range(1,n):
        xi=a+i*h
        sum=sum+f(xi)        
    return part1+h*sum

def adaptivetrap(f,a,b,ep):
    
    max=0    
    step=float(abs(a-b)/1000)
    i=0
    while (i<1000):
        i=i+1        
        adj=a;
        adj=a+step*i;
        dval=diff2(f,adj)
        if(abs(dval)>max):
            max=abs(dval)
            
    h=sqrt(12*ep)*((b-a)*max)**.5
    n=(b-a)/h        
    return trapezint(f,a,b,int(ceil(n)))
    
print adaptivetrap(f1,0.0,10.0,1E-5)

The function, adaptivetrap, above takes the f1 we defined above and integrates it using the adaptive trapezoid rule from 0 to 10 with an error tolerance of 1E-5.