Formule de Simpson

On rappelle la formule de Simpson pour le calcul approché d’une intégrale :

\[\begin{equation*} \begin{split} \int_a^bf(x)\mathrm dx \approx \Delta_x\sum_{i=0}^{N-1}\frac{f(x_i) + 4f(x_i+\frac{\Delta_x}2) + f(x_i + \Delta_x)}6 \equiv I(f;a,b,N) \end{split} \end{equation*}\]

où \(\Delta_x =\frac{b-a}N\) et \(x_i = a + i\Delta_x\).

1. Ecrivez une fonction simpson_slow qui calcule l’intégrale d’une function \(f\) entre \(a\) et \(b\) avec \(N\) pas avec la méthode de Simpson en utilisant une boucle (for loop).

2. Ecrivez une autre fonction simpson_fast qui fait la même chose sans utiliser de boucle (on suposera que la fonction \(f\) est vectorisée).

3. Calculez l’intégrale de \(f(x) = \frac{1}{1+x^2}\) entre 0 et 1 pour \(N=1000\) et comparez le temps entre les deux fonctions en ecrirant %timeit avant la commande.

4. Calculez la valeur théorique \(I^\star\) de l’intégrale et tracez en échelle logarithmique le residu \(|I^\star - I(f,0,1,N)|\) par rapport à \(N\). Qu’elle est la vitesse de convergence de cette intégrale ?

---

Solution

# 1.

def simpson_slow(f, a, b, N):
    N = int(N)
    dx = (b-a)/N
    x = a + np.arange(N)*dx
    integral = 0
    for i in range(N):
        integral += (f(x[i]) + 4*f(x[i]+dx/2) + f(x[i]+dx))*dx/6
    return integral

# 2.

def simpson_fast(f, a, b, N):
    N = int(N)
    dx = (b-a)/N
    x = a + np.arange(N)*dx
    return np.sum(f(x) + 4*f(x+dx/2) + f(x+dx))*dx/6

# 3.

def f(x):
    return 1/(1+x**2)

print("Temps d'exécution pour simpson_slow:")
%timeit simpson_slow(f, 0, 1, 1000)
print('\n')
print("Temps d'exécution pour simpson_fast:")
%timeit simpson_fast(f, 0, 1, 1000)

Temps d'exécution pour simpson_slow:
5.59 ms ± 250 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)


Temps d'exécution pour simpson_fast:
60.6 µs ± 708 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

La solution analytique de l’intégrale est : \(\int_0^1 \frac{dx}{1+x^2} = \arctan(x)|_0^1 = \pi/4\)

N_array = np.arange(6,30)
residuum = []
for N in N_array:
    residuum.append(np.pi/4-simpson_fast(f, a=0, b=1, N=N))
# Remarque : utiliser une boucle for ici ne ralentit pas le code car
# 1/ le calcul de simpson prend relativement beaucoup de temps
# 2/ il y a moins de 30 itérations

# 4.

plt.loglog(N_array, residuum, 'bo', label='Residuum')
plt.title('Convergence de la méthode de Simpson')
plt.grid(which = 'both')
plt.ylabel('Residus')
plt.xlabel('Nombre de pas N')

TN = np.logspace(0.7, 1.5)
plt.loglog(TN, residuum[0]*(N_array[0]/TN)**6, label='$1/N^6$')

plt.legend()
#gradient = np.polyfit(np.log(N_array), np.log(residuum), 1)[0]
#print(f"L'intégrale converge en res(N)~ N^({gradient:0.1f})")

<matplotlib.legend.Legend at 0x7fca3c2a9520>