Numerical integration with Python

In this short article I am going to post a simple Python script for numerical integration.

Numerical integration is a part of a family of algorithms for calculating the numerical value of a definite integral. The two simplest method for performing numerical integration are rectangle and trapezoidal rule. If you’d like to know more about these two methods you can check the wikipedia page which has been pretty helpful to me.

Now, the code below works fine if the functions you are using is not decreasing and positive, however:
-If your function is decreasing, the plot need to be adjusted to show the correct graph
-If the function is negative at some point, the code might not work, since it has not been designed to handle it. Perhaps in the future I’ll fix that.

Therefore, bottom line, square root of x is fine, logarithms (after 1),x^2,x^3(from 0) etc…

Here is an approximation with 50 points and the rectangle rule:

and the code output on the console:

	# Output
	# >>> ================================ RESTART ================================
	# >>>
	# Actual area: 5208.333333333333
	#
	# Approximating using rectangular rule...
	# Percentage error: -0.0104123281966 %
	# Approximation: (5207.7910245730936, -0.010412328196596356)
	#
	# Approximating using trapezoidal rule...
	# Percentage error: 0.0208246563932 %
	# Approximation: (5209.4179508538109, 0.020824656393175246)
	# #############################################################
	# Approximating using rectangular rule with: 15 points, percentage error: -1.5625 %
	# Approximating using rectangular rule with: 25 points, percentage error: -0.127551020408 %
	# Approximating using rectangular rule with: 35 points, percentage error: -0.0434027777778 %
	# 5206.07277199
	# Approximating using trapezoidal rule with: 15 points, percentage error: 3.125 %
	# Approximating using trapezoidal rule with: 25 points, percentage error: 0.255102040816 %
	# Approximating using trapezoidal rule with: 35 points, percentage error: 0.0868055555556 %
	# Approximating using trapezoidal rule with: 45 points, percentage error: 0.0432525951557 %
	# 5210.58607266
	# >>>

view raw numerical_integration_out.py hosted with ❤ by GitHub

Finally, below is the code I used to generate the output above:

	import numpy as np
	import matplotlib.pyplot as plt
	plt.style.use('ggplot')
	# Initial function given
	def f(x):
	return x**2
	# Definite integral of the function from a to b
	def F(a,b):
	if a > b:
	raise ValueError('b must be greater than a')
	elif a == b:
	return 0
	else:
	return (b**3-a**3)/3
	# Approximating function using numerical methods:
	# rectangular
	# trapezoidal
	def approximateNumerical(a,b,points=10,error=False,mod='rectangular',plt_data=False):
	if points < 2:
	raise ValueError('Number of points must be greater than 2')
	if a == b:
	return 0
	n = np.linspace(a,b,points)
	partialSum = 0
	if mod == 'rectangular':
	def miniArea(c,d):
	return (d-c)*f((c+d)/2)
	elif mod == 'trapezoidal':
	def miniArea(c,d):
	return (d-c)*(f(c)+f(d))/2
	else:
	raise ValueError('Method '+mod+' unknown')
	for i in range(1,len(n)):
	partialSum += miniArea(n[i-1],n[i])
	e = (partialSum-F(a,b))/F(a,b) *100
	if error:
	print('\nApproximating using '+ mod+' rule...')
	print('Percentage error: ',e,'%')
	if plt_data:
	plot_dat(a,b,points,mod=mod)
	return partialSum,e
	# Plotting function for a visual representation
	def plot_dat(a,b,points,mod='rectangular'):
	n = np.linspace(a,b,points)
	plt.plot(n,f(n),color='red')
	plt.xlabel('x')
	plt.ylabel('y')
	plt.title('Numerical approximation: '+mod)
	if mod == 'rectangular':
	for i in range(1,len(n)):
	c = n[i-1]
	d = n[i]
	plt.plot([c,c],[0,f((c+d)/2)],color='blue')
	plt.plot([d,d],[0,f(d)],color='blue')
	plt.plot([c,d],[f((c+d)/2),f((c+d)/2)],color='blue')
	plt.show()
	if mod == 'trapezoidal':
	for i in range(1,len(n)):
	c = n[i-1]
	d = n[i]
	plt.plot([d,d],[0,f(d)],color='blue')
	plt.plot([c,c],[0,f(c)],color='blue')
	plt.plot([c,d],[f(c),f(d)],color='blue')
	plt.show()
	return 0
	# Approximate area with a given precision
	def approxGivenPrecision(a,b,error=0.5,mod='rectangular',printf=False):
	e = 100
	p = 5
	itMax = 500
	it = 0
	while abs(e) > error:
	area,e = approximateNumerical(a,b,mod=mod,points=p)
	p += 10
	it += 1
	if printf:
	print('Approximating using '+mod+' rule with:',p,'points, percentage error:',e,'%')
	if it > itMax:
	print('Number of iterations exceeded: '+str(itMax))
	break
	return area
	#-----------------------------Run the program----------------------------------
	# Initial parameters:
	# xmin = m
	# xmax = M
	# points used = p
	m = 0
	M = 25
	p = 50
	print('Actual area:',F(m,M))
	print('Approximation:',approximateNumerical(m,M,p,error=True,plt_data=True))
	print('Approximation:',approximateNumerical(m,M,p,mod='trapezoidal',error=True,plt_data=True))
	print('#############################################################')
	print(approxGivenPrecision(m,M,error=0.05,printf=True))
	print(approxGivenPrecision(m,M,error=0.05,mod='trapezoidal',printf=True))

view raw numerical_integration.py hosted with ❤ by GitHub

Hope this is useful.

Tag di Technorati: maths,python,calculus,algebra,programming,integral,integrals