Numerical Integration in Python
Numerical integration aims to find the area under a curve without using analytical methods.
Rectangular
The rectangle rule states that:
In other words, the area under a curve \(f(x)\) between a point \(a\) and a point \(b\) is roughly equal to the area of a rectangle with width and height \(\frac{f(a)+f(b)}{2}\).
This approach is easily implemented in Python:
def rectint(f,a,b,rectangles):
cumulative_area=0
a=float(a)
b=float(b)
rectangles=float(rectangles)
i=(b-a)/rectangles
trailing_x=a
leading_x=a+i
while (a<=leading_x<=b) or (a>=leading_x>=b):
area=f((trailing_x+leading_x)/2)*i
cumulative_area+=area
leading_x+=i
trailing_x+=i
return cumulative_area
Trapezoidal
The trapezoid rule is a small departure from the rectangle rule and tends to produce more accurate approximations:
The Python code is itself similar to that of the rectangle rule:
def trapint(f,a,b,trapezoids):
cumulative_area=0
a=float(a)
b=float(b)
trapezoids=float(trapezoids)
i=(b-a)/trapezoids
trailing_x=a
leading_x=a+i
while (a<=leading_x<=b) or (a>=leading_x>=b):
area=(f(trailing_x)+f(leading_x)))*i/2
cumulative_area+=area
leading_x+=i
trailing_x+=i
return cumulative_area
In each code example, rules for negativity are preserved in line with the FTC. Predictably, the more rectangles or trapezoids requested, the lower the error. Indeed, the limit definition of an integral is based around infinite subintervals.
Comparison
Let’s test the functions with an integral which is easily evaluated analytically:
With rectint
:
>>> def func(x):
... return x**2-2*x+3
>>> rectint(func,-5,10,15):
343.75
>>> rectint(func,-5,10,100):
332.72353125
>>> rectint(func,-5,10,1000):
344.99971875
>>> rectint(func,-5,10,10000):
344.875517437
>>> rectint(func,-5,10,100000):
344.999999971
With trapint
:
>>> def func(x):
... return x**2-2*x+3
>>> trapint(func,-5,10,15):
347.5
>>> trapint(func,-5,10,100):
332.8070625
>>> trapint(func,-5,10,1000):
345.0005625
>>> trapint(func,-5,10,10000):
344.875525873
>>> trapint(func,-5,10,100000):
345.000000056
QED.