https://docs.scipy.org/doc/scipy/tutorial/integrate.html
https://en.wikipedia.org/wiki/Numerical_integration
https://en.wikipedia.org/wiki/Trapezoidal_rule
https://en.wikipedia.org/wiki/Simpson%27s_rule
The 'scipy.integrate' sub-package provides several integration techniques including an ordinary differential equation integrator.
import numpy as np import scipy.integrate as integrate # help(integrate) # The function quad() is provided to integrate a function of one variable # between two points. # \integ_x1^x2 dx f(x) - a definite integral # integrate.quad(f, x1, x2) # The first argument to quad is a callable Python object. # x1 and x2 are the limits of integration. # Infinite inputs are also allowed in quad by using 'np.inf'. f = lambda x: np.exp(-x) result = integrate.quad(f, 0, 1) # (0.6321205588285578, 7.017947987503856e-15) # 0.632_120_558_828_5578(70) # exact = (np.e - 1) / np.e = 0.632_120_558_828_5577(1) # np.e = 2.718_281_828_459_045 # 16 digits
from scipy.special import jv # Bessel function of the first kind f = lambda x: jv(0, x) result = integrate.quad(f, 0, 5) # (0.7153119177847678, 2.47260738289741e-14) # 0.715_311_917_784_768(25)
# Double integration with the dblquad() function. # \integ_y1^y2 dy \integ_x1(y)^x2(y) dx f(x,y) # integrate.dblquad(f, y1, y2, x1, x2) # Note that x1=g(y) an x2=h(y) can be functions of y. # y2 (0,1)+---+(1,1) # | | # y1 (0,0)+---+(1,0) # x1 x2 f = lambda x, y: np.exp(-x**2 -y**2) y1, y2 = 0, 1 x1 = lambda y: 0 x2 = lambda y: 1 area = integrate.dblquad(f, y1, y2, x1, x2) # (0.5577462853510337, 8.291374381535408e-15) # 0.557_746_285_351_0337(83)
# f(x,y) = x*y, the area is a triangle with boundaries: # x = 0, y = 0, x + 2*y = 1, # vertices are (0, 0), (1, 0), (0, 0.5). # (0,0.5)+ # | \ # | \ # (0,0)+---+(1,0) f = lambda x, y: x*y y1, y2 = 0, 0.5 x1 = lambda y: 0 x2 = lambda y: 1-2*y # x = 1 - 2*y area = integrate.dblquad(f, y1, y2, x1, x2) # (0.010416666666666668, 4.101620128472366e-16) # 0.010_416_666_666_666_67(41) # exact = 1. / 96. = 0.010_416_666_666_666_666(1)
# Triple integration with the tplquad() function.
# For n-fold integration, scipy provides the function nquad.
# integrate.nquad(f, iterable_with_bounds)
def bounds_y():
return [0, 0.5]
def bounds_x(y):
return [0, 1-2*y]
area3 = integrate.nquad(lambda x, y: x*y, [bounds_x, bounds_y])
The purpose of the numerical integration is to calculate the numerical value of a definite integral to a given degree of accuracy. The integrand is evaluated at a finite set of points called 'integration points' and a weighted sum of these values is used to approximate the integral.
I = \integ_x1^x2 dx f(x) The midpoint rule or rectangle rule (constant function approximation) I = (x2-x1) f((x1+x2) / 2) The trapezoidal rule (linear function approximation) I = (x2-x1) [f(x1) + f(x2)] / 2 Simpson's 1/3 rule (quadratic function approximation) I = [(x2-x1) / 6] [f(x1) + 4 f((x1+x2) / 2) + f(x2)]
# Integrating functions using fixed samples. # trapz - Newton-Coates formulas (straight lines between adjacent points) # simps - Simpson’s rule (parabolas between three adjacent points) # romb - Romberg integration (2**k+1 points are required) # \integ_x1^x2 dx f(x) # \integ_0^2 dx sqrt(x) f = lambda x: np.sqrt(x) x1, x2 = 0, 2 x = np.linspace(x1, x2, 2**5+1) y = f(x) res_trapz = integrate.trapz(y, x) res_simps = integrate.simps(y, x) # recommended res_romb = integrate.romb(y, dx=x[1]-x[0]) # 2**k+1 points are required res_exact = (2.0/3.0) * (x2-x1)**(3.0/2.0) # for 2**5+1=33 points # res_trapz 1.882_484_948_488_58 # res_simps 1.884_349_579_078_4088 # res_romb 1.884_546_200_832_0417 # res_exact 1.885_618_083_164_1267 # for 2**6+1=65 points # res_trapz 1.884_498_435_329_6215 # res_simps 1.885_169_597_609_9685 # res_romb 1.885_239_284_741_2136 # res_exact 1.885_618_083_164_1267