[Calcul scientifique / technique par Python] Intégration numérique, loi trapézoïdale / Simpson, calcul numérique, scipy

L'intégration numérique de données discrètes est réalisée en utilisant la méthode cumtrapz (loi trapézoïdale) et la méthode simps (loi de Simpson) de scipy.integrate. À titre d'exemple, considérons $ \ int_0 ^ 1 \ frac {4} {1 + x ^ 2} dx = \ pi $.

Contenu

(1.A) Code utilisant scipy. ** C'est quand vous êtes pressé. ** ** (2) Addendum sur la précision de calcul de la loi trapézoïdale et de la loi de Simpson. (3) Implémentation simple de règles trapézoïdales et Simpson avec du code python. Utile lorsque vous souhaitez connaître la procédure de calcul.

(1) Code utilisant scipy

(1) Code d'utilisation Scipy. Concis.

from scipy import integrate
import numpy as np

x=np.linspace(0,1,5)  # [0,1]Stockez la grille divisée en 5 parties égales en x
y=4/(1+x**2)          #Fonction intégrale de l'intégration numérique

y_integrate_trape = integrate.cumtrapz(y, x)  #Calcul numérique intégral par loi trapézoïdale
y_integrate_simps = integrate.simps(y, x) #Calcul de l'intégrale numérique par la loi de Simpson

print(y_integrate_trape[-1]) #Voir les résultats
print(y_integrate_simps)     #Voir les résultats

Résultat (1): Utilisez scipy

3.13117647059   #Règle trapézoïdale
3.14156862745   #Loi de Simpson
3.14159265358979...Valeur exacte(π)

(2) [Addendum] Comparaison de la précision et de la vitesse de convergence de la loi trapézoïdale et de la loi de Simpson

Comme cela est bien connu, pour un pas h suffisamment petit, L'erreur d'intégration selon la loi trapézoïdale est $ O (h ^ 2) $, et celle de la loi de Simpson est $ O (h ^ 3) $. En supposant que le nombre de grilles est N, ce sont respectivement $ O (N ^ {-2}) $ et $ O (N ^ {-3}) $. Il est éducatif de vérifier cela directement à travers cet exemple.

Voici le code pour le confirmer. Les deux tracés logarithmiques sont réalisés avec l'erreur relative due à l'intégration numérique sur l'axe y et le nombre de grilles N sur l'axe horizontal. Dans la zone où la relation ci-dessus se vérifie, $ log y \ propto n logN $, $ n = -2 $ correspond à la loi trapézoïdale, et $ n = -3 $ correspond à la loi de Simpson.

(2)Comparaison de la précision des calculs

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] #Stocker le nombre de grilles de 4 à 2048 dans la liste nn
for j in nn:
    x=np.linspace(0,1, j)
    y=4/(1+x**2)
    y_integrate_trape = integrate.cumtrapz(y, x) #Calcul numérique intégral par loi trapézoïdale
    y_integrate_simps = integrate.simps(y, x)     #Calcul de l'intégrale numérique par la loi de Simpson
    err_trape.append(abs(np.pi-y_integrate_trape[-1])/np.pi)  #Évaluation de l'erreur relative de l'intégration numérique par loi trapézoïdale
    err_simps.append(abs(np.pi-y_integrate_simps)/np.pi)     #Évaluation de l'erreur relative de l'intégration numérique par la loi de Simpson

# for plot
ax = plt.gca()
ax.set_yscale('log')  #Dessinez l'axe y sur l'échelle logarithmique
ax.set_xscale('log')  #Dessinez l'axe des x sur l'échelle logarithmique
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") #Affichage de grille."both"Dessine une grille sur les deux axes xy.

plt.legend(loc='upper right')


plt.show()

La pente de la droite sur la figure correspond à $ n $ d'ordre de convergence. Comme prévu, la loi trapézoïdale (ligne bleue) est $ n \ sim-2 $ et la loi de Simpson (ligne rouge) est $ n \ sim-3 $.

(3) Implémentation simple de règles trapézoïdales et Simpson avec du code python.

Utilisation: préparer un fichier de données XY ('xy.dat')

python3 fuga.py xy.dat

Ensuite, le résultat du calcul numérique par la règle trapézoïdale et la règle Simpson s'affiche.

fuga.py

import os, sys, math


def integrate_trape(f_inp): #Fonction d'intégration numérique par loi trapézoïdale
    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): #Fonction d'intégration numérique par la loi de Simpson
    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()