Comment réaliser une intégration numérique avec python ?

Exemple de comment calculer une intégrale à partir d'un tableau de données sous python (pour faire faire une intégration sur une fonction définie voir Calculer une intégrale simple avec python). Dans cet exemple on suppose que l'on a deux matrices x et y et on veut calculer l'intégrale y(x) (pour illustrer on a pris $y=x^2$). Pour calculer l'intégrale on peut alors utiliser la méthode des trapèzes numpy.trapz ou la méthode de Simpson scipy.integrate.simps:

from scipy.integrate import simps
from numpy import trapz

import numpy as np

def function(x):
    return x**2

x = np.arange(1,10,0.1)
y = function(x)

print x
print y

# primitive :

print "area: ", 1.0 / 3.0 * ( x[len(x)-1]**3 - x[0]**3 )

# using Trapezoidal rule:

area = trapz(y,x)
print 'area: ',area

# using Simpson's rule:

area = simps(y,x)
print 'area: ',area

donne:

area from the primitive:  323.099666667
area from Trapezoidal rule:  323.1145
area from Simpson's rule:  323.099833333

Références