[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy

Numerical integration of discrete data is performed using the cumtrapz method (trapezoidal law) and simps method (Simpson's law) of scipy.integrate. As an example, consider $ \ int_0 ^ 1 \ frac {4} {1 + x ^ 2} dx = \ pi $.

Contents

(1.A) Code using scipy. ** This is when you are in a hurry. ** ** (2) Addendum on the calculation accuracy of trapezoidal rule and Simpson rule. (3) Simple implementation of trapezoidal rule and Simpson rule by python code. Useful when you want to know the calculation procedure.

(1) Code using scipy

(1) Scipy usage code. Concise.

from scipy import integrate
import numpy as np

x=np.linspace(0,1,5)  # [0,1]Store the grid divided into 5 equal parts in x
y=4/(1+x**2)          #Integrand of numerical integration

y_integrate_trape = integrate.cumtrapz(y, x)  #Numerical integration calculation by trapezoidal law
y_integrate_simps = integrate.simps(y, x) #Numerical integration calculation by Simpson's law

print(y_integrate_trape[-1]) #View results
print(y_integrate_simps)     #View results

Result (1): Use scipy

3.13117647059   #Trapezoidal rule
3.14156862745   #Simpson's law
3.14159265358979...Exact value(π)

(2) [Addendum] Comparison of accuracy and convergence speed of trapezoidal law and Simpson's law

As is well known, for a sufficiently small step size h, The integration error according to the trapezoidal rule is $ O (h ^ 2) $, and that of the Simpson rule is $ O (h ^ 3) $. Assuming that the number of grids is N, these are $ O (N ^ {-2}) $ and $ O (N ^ {-3}) $, respectively. It is educational to verify this directly through this example.

Below is the code to confirm it. A log-log plot is made with the relative error due to numerical integration on the y-axis and the grid number N on the horizontal axis. In the area where the above relationship holds, $ log y \ propto n logN $, $ n = -2 $ corresponds to the trapezoidal law, and $ n = -3 $ corresponds to the Simpson law.

(2)Comparison of calculation accuracy

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

err_trape=[]
err_simps=[]
nn=[4,8,16,32,64,128,256,512,1024,2048] #Store the number of grids from 4 to 2048 in list nn
for j in nn:
    x=np.linspace(0,1, j)
    y=4/(1+x**2)
    y_integrate_trape = integrate.cumtrapz(y, x) #Numerical integration calculation by trapezoidal law
    y_integrate_simps = integrate.simps(y, x)     #Numerical integration calculation by Simpson's law
    err_trape.append(abs(np.pi-y_integrate_trape[-1])/np.pi)  #Relative error evaluation of numerical integration by trapezoidal law
    err_simps.append(abs(np.pi-y_integrate_simps)/np.pi)     #Relative error evaluation of numerical integration by Simpson's law

# for plot
ax = plt.gca()
ax.set_yscale('log')  #Draw y-axis on log scale
ax.set_xscale('log')  #Draw x-axis on log scale
plt.plot(nn,err_trape,"-", color='blue', label='Trapezoid rule')
plt.plot(nn,err_simps,"-", color='red', label='Simpson rule')


plt.xlabel('Number of grids',fontsize=18)
plt.ylabel('Error (%)',fontsize=18)
plt.grid(which="both") #Grid display."both"Draws a grid on both xy axes.

plt.legend(loc='upper right')


plt.show()

The slope of the straight line in the figure corresponds to the convergence order $ n $. As expected, the trapezoidal rule (blue line) is $ n \ sim-2 $ and the Simpson rule (red line) is $ n \ sim-3 $.

(3) Simple implementation of trapezoidal rule and Simpson rule by python code.

Usage: Prepare an XY data file ('xy.dat')

python3 fuga.py xy.dat

Then, the numerical calculation result by the trapezoidal rule and Simpson's rule is displayed.

fuga.py

import os, sys, math


def integrate_trape(f_inp): #Function of numerical integration by trapezoidal law
    x_lis=[]; y_lis=[]
#   f_inp.readline()
    for line in f_inp:
        x_lis.append(float(line.split()[0]))
        y_lis.append(float(line.split()[1]))
    sum_part_ylis=sum(y_lis[1:-2]) 
    del_x=(x_lis[1]-x_lis[0])
    integrated_value = 0.5*del_x*(y_lis[0]+y_lis[-1]+2*sum_part_ylis) 
    print("solution_use_trapezoid_rule: ", integrated_value)

def integrate_simpson(f_inp): #Function of numerical integration by Simpson's law
    x_lis=[]; y_lis=[]
    n_counter = 0
    for line in f_inp:
        x_lis.append(float(line.split()[0]))
        if n_counter % 2 == 0:
            y_lis.append(2*float(line.split()[1]))
        else:
            y_lis.append(4*float(line.split()[1]))
        n_counter += 1
    sum_part_ylis=sum(y_lis[1:-2]) # sum from y_1 to y_n-1 
    del_x=(x_lis[1]-x_lis[0])
    integrated_value = del_x*(y_lis[0]+y_lis[-1]+sum_part_ylis)/3 
    print("solution_use_Simpson_formula: ", integrated_value)

##
def main():
    fname=sys.argv[1]

    f_inp=open(fname,"r")
    integrate_trape(f_inp) 
    f_inp.close()

    f_inp=open(fname,"r")
    integrate_simpson(f_inp) 
    f_inp.close()

if __name__ == '__main__':
    main()