Intégration numériques

Ici \(h = 0{,}25\) et \(x_k = 0\,;\,0{,}25\,;\,0{,}5\,;\,0{,}75\,;\,1\).

1. Rectangles.

• A gauche : \(0{,}25\,(0 + 0{,}0625 + 0{,}25 + 0{,}5625) = 0{,}21875\). Erreur : \(0{,}1146\).
• A droite : \(0{,}25\,(0{,}0625 + 0{,}25 + 0{,}5625 + 1) = 0{,}46875\). Erreur : \(0{,}1354\).
• Point milieu (\(x = 0{,}125\,;\,0{,}375\,;\,0{,}625\,;\,0{,}875\)) : \(0{,}25\,(0{,}015625 + 0{,}140625 + 0{,}390625 + 0{,}765625) = 0{,}328125\). Erreur : \(0{,}0052\).

2. Trapezes. \(T_4 = 0{,}25\left(\frac{0}{2} + 0{,}0625 + 0{,}25 + 0{,}5625 + \frac{1}{2}\right) = 0{,}34375\). Erreur : \(0{,}0104\).

3. Simpson avec \(n = 2\), donc \(h = 0{,}5\) : \[ S_2 = \frac{0{,}5}{3}\bigl(f(0) + 4f(0{,}5) + f(1)\bigr) = \frac{0{,}5}{3}(0 + 1 + 1) = \frac13. \] Erreur : 0 (aux arrondis pres).

4. Bilan.

Trois observations :

• Le point milieu est exactement deux fois plus precis que les trapezes (\(5{,}2\cdot10^{-3}\) contre \(1{,}04\cdot10^{-2}\)), conformement aux constantes \(1/24\) et \(1/12\), et pour une evaluation de moins.
• Leurs erreurs sont de signes opposes : le point milieu sous-estime, les trapezes surestiment (car \(x^2\) est convexe). D'ou \(\frac{T + 2M}{3} = \frac13\) exactement : c'est Simpson.
• Simpson est exacte ici parce que \(x^2\) est de degre \(2 \le 3\), et avec seulement 3 evaluations.

import numpy as np
from math import exp, e

def rectangles_gauche(f, a, b, n):
    h = (b-a)/n
    return h * sum(f(a + k*h) for k in range(n))

def point_milieu(f, a, b, n):
    h = (b-a)/n
    return h * sum(f(a + (k+0.5)*h) for k in range(n))

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b)
                  + 4*sum(f(a+k*h) for k in range(1, n, 2))
                  + 2*sum(f(a+k*h) for k in range(2, n, 2)))

# 2. degres d'exactitude
print(f"{'degre':>6}{'gauche':>12}{'milieu':>12}{'trapezes':>12}{'Simpson':>12}")
for d in range(6):
    f = lambda x, d=d: x**d
    exact = 1/(d+1)
    ligne = [abs(m(f, 0, 1, 6) - exact) for m in
             (rectangles_gauche, point_milieu, trapezes, simpson)]
    print(f"{d:>6}" + "".join(f"{v:>12.2e}" for v in ligne))

 degre      gauche      milieu    trapezes     Simpson
     0    0.00e+00    0.00e+00    0.00e+00    0.00e+00
     1    8.33e-02    0.00e+00    0.00e+00    0.00e+00
     2    8.33e-02    2.31e-03    4.63e-03    0.00e+00
     3    1.24e-01    6.94e-03    1.39e-02    5.55e-17
     4    1.51e-01    1.30e-02    2.62e-02    5.14e-05
     5    1.70e-01    1.99e-02    4.05e-02    2.14e-04

2. Les degres d'exactitude sont confirmes : rectangles exacts jusqu'au degre 0, point milieu et trapezes jusqu'au degre 1, Simpson jusqu'au degre 3. Le point milieu, avec un seul point, atteint le meme degre que les trapezes qui en utilisent deux.

# 3-4. ordres experimentaux
I = e - 1
print(f"\n{'n':>6}{'gauche':>13}{'r':>6}{'milieu':>13}{'r':>6}"
      f"{'trapezes':>13}{'r':>6}{'Simpson':>13}{'r':>6}")
prec = None
for k in range(1, 11):
    n = 2**k
    E = [abs(m(exp, 0, 1, n) - I) for m in
         (rectangles_gauche, point_milieu, trapezes, simpson)]
    if prec is None:
        print(f"{n:>6}" + "".join(f"{v:>13.2e}{'':>6}" for v in E))
    else:
        s = f"{n:>6}"
        for v, p in zip(E, prec):
            s += f"{v:>13.2e}{(p/v if v > 0 else float('inf')):>6.1f}"
        print(s)
    prec = E

     n       gauche     r       milieu     r     trapezes     r      Simpson     r
     2     4.26e-01          2.15e-02          4.31e-02          5.79e-05
     4     2.28e-01   1.9   5.36e-03   4.0   1.08e-02   4.0   3.59e-06  16.1
     8     1.18e-01   1.9   1.34e-03   4.0   2.69e-03   4.0   2.23e-07  16.1
    16     6.01e-02   2.0   3.35e-04   4.0   6.72e-04   4.0   1.39e-08  16.0
    32     3.03e-02   2.0   8.38e-05   4.0   1.68e-04   4.0   8.70e-10  16.0
   256     3.80e-03   2.0   1.31e-06   4.0   2.62e-06   4.0   2.12e-13  16.0
  1024     9.51e-04   2.0   8.19e-08   4.0   1.64e-07   4.0   8.44e-15   9.3

4. Les rapports valent 2, 4, 4 et 16, soit des ordres \(\log_2 2 = 1\), \(\log_2 4 = 2\), 2 et \(\log_2 16 = 4\). Conforme a la theorie.

A la derniere ligne, le rapport de Simpson tombe a 9,3 : l'erreur \(8\cdot10^{-15}\) est proche de \(\varepsilon\), et l'erreur d'arrondi commence a dominer. C'est la limite de precision atteignable, conformement au chapitre 20.

1. Pour \(f(x) = e^x\) sur \([0,1]\) : \(f'' = e^x\), donc \(\max|f''| = e\). La majoration des trapezes donne \[ |I - T_n| \le \frac{(b-a)h^2}{12}\max|f''| = \frac{e}{12n^2}. \] Pour que ce majorant soit \(< 10^{-6}\) : \[ n > \sqrt{\frac{e}{12\cdot10^{-6}}} = \sqrt{226\,516} \approx 476. \] On prend \(n = 476\).

from math import exp, e, ceil, sqrt

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b) + 4*sum(f(a+k*h) for k in range(1, n, 2))
                              + 2*sum(f(a+k*h) for k in range(2, n, 2)))


I = e - 1

n_trap = ceil(sqrt(e / (12 * 1e-6)))
err = abs(trapezes(exp, 0, 1, n_trap) - I)
majorant = e / (12 * n_trap**2)
print(f"trapezes : n = {n_trap}, erreur reelle {err:.3e}, "
      f"majorant {majorant:.3e}, rapport {majorant/err:.2f}")

# 3. Simpson : f^(4) = e^x, meme majorant
n_simp = ceil((e / (180 * 1e-6))**0.25)
if n_simp % 2: n_simp += 1
err_s = abs(simpson(exp, 0, 1, n_simp) - I)
maj_s = e / (180 * n_simp**4)
print(f"Simpson  : n = {n_simp}, erreur reelle {err_s:.3e}, "
      f"majorant {maj_s:.3e}, rapport {maj_s/err_s:.2f}")

trapezes : n = 476, erreur reelle 1.000e-06, majorant 1.000e-06, rapport 1.00
Simpson  : n = 12, erreur reelle 8.815e-08, majorant 7.288e-07, rapport 8.27

2. Pour les trapezes, le majorant est atteint presque exactement (rapport 1,00) : c'est parce que \(f'' = e^x\) est monotone et que la formule \(I - T_n = -\frac{(b-a)h^2}{12}f''(\xi)\) est une egalite. La majoration est donc optimale ici.

3. Pour Simpson, \(f^{(4)} = e^x\) et \(|I - S_n| \le \frac{e}{180 n^4}\). Pour \(10^{-6}\) : \(n > \left(\frac{e}{180\cdot10^{-6}}\right)^{1/4} \approx 11{,}0\), donc \(n = 12\) (pair). Le majorant est ici 8 fois trop pessimiste, parce que \(\xi\) tombe pres du milieu de l'intervalle ou \(e^\xi \ll e\).

4. Comparaison du cout.

Pour une precision de \(10^{-12}\), l'ecart serait bien plus spectaculaire : \(n \approx 476\,000\) pour les trapezes contre \(n \approx 380\) pour Simpson, soit un facteur 1250. L'ordre de la methode compte infiniment plus que l'optimisation de son code.

import math
from math import exp, sin, sqrt, pi, log

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b)
                  + 4*sum(f(a+k*h) for k in range(1, n, 2))
                  + 2*sum(f(a+k*h) for k in range(2, n, 2)))

# 1. gaussienne
I1 = simpson(lambda t: exp(-t*t), 0, 1, 1000)
exact1 = math.erf(1) * sqrt(pi) / 2
print(f"int exp(-t^2) = {I1:.15f}, exact {exact1:.15f}, "
      f"erreur {abs(I1-exact1):.2e}")

# 2. sinus cardinal : prolongement par continuite en 0
def sinc(t):
    return 1.0 if t == 0 else sin(t)/t
I2 = simpson(sinc, 0, 1, 1000)
print(f"int sin(t)/t  = {I2:.15f}   (Si(1) = 0.946083070367183)")

# 3.
I3 = simpson(lambda t: 1/(1+t**4), 0, 1, 2000)
exact3 = (pi + 2*log(1+sqrt(2))) / (4*sqrt(2))
print(f"int 1/(1+t^4) = {I3:.15f}, exact {exact3:.15f}, "
      f"erreur {abs(I3-exact3):.2e}")

# 4.
from math import cos
I4 = simpson(lambda t: 1/(2+cos(t)), 0, pi, 2000)
exact4 = pi / sqrt(3)
print(f"int 1/(2+cos)  = {I4:.15f}, exact {exact4:.15f}, "
      f"erreur {abs(I4-exact4):.2e}")

int exp(-t^2) = 0.746824132812427, exact 0.746824132812427, erreur 1.11e-16
int sin(t)/t  = 0.946083070367183   (Si(1) = 0.946083070367183)
int 1/(1+t^4) = 0.866972987339911, exact 0.866972987339911, erreur 2.22e-16
int 1/(2+cos)  = 1.813799364234218, exact 1.813799364234218, erreur 4.44e-16

2. Le point delicat. La fonction \(t \mapsto \sin(t)/t\) n'est pas definie en 0 : Python leverait ZeroDivisionError. Or elle admet une limite finie (\(\sin t/t \to 1\)) et se prolonge en une fonction analytique. Deux solutions :

• tester explicitement \(t = 0\) et renvoyer 1, comme ci-dessus ;
• utiliser la methode du point milieu, qui n'evalue jamais aux bornes.

def point_milieu(f, a, b, n):
    h = (b-a)/n
    return h * sum(f(a + (k+0.5)*h) for k in range(n))

print(point_milieu(lambda t: sin(t)/t, 0, 1, 100000))   # aucune erreur en 0

from math import exp, sqrt, pi, erf

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b) + 4*sum(f(a+k*h) for k in range(1, n, 2))
                              + 2*sum(f(a+k*h) for k in range(2, n, 2)))


def adaptative(f, a, b, eps=1e-10, prof=50):
    """Simpson adaptative. Renvoie (valeur, evaluations, points evalues)."""
    points = []
    def g(x):
        points.append(x)
        return f(x)

    def simp(a, b):
        m = (a+b)/2
        return (b-a)/6 * (g(a) + 4*g(m) + g(b))

    def aux(a, b, eps, entier, prof):
        m = (a+b)/2
        G, D = simp(a, m), simp(m, b)
        if prof <= 0 or abs(G + D - entier) <= 15*eps:
            return G + D + (G + D - entier)/15
        return aux(a, m, eps/2, G, prof-1) + aux(m, b, eps/2, D, prof-1)

    total = simp(a, b)
    v = aux(a, b, eps, total, prof)
    return v, len(points), points

f = lambda x: exp(-100*(x-0.3)**2)
exact = sqrt(pi)/20 * (erf(10*0.7) + erf(10*0.3))     # valeur analytique

v, ev, pts = adaptative(f, 0, 1, 1e-12)
print(f"adaptative : {v:.14f}, erreur {abs(v-exact):.2e}, {ev} evaluations")

# 3. Simpson uniforme, meme precision
n = 2
while abs(simpson(f, 0, 1, n) - exact) > abs(v - exact):
    n *= 2
print(f"uniforme   : n = {n}, soit {n+1} evaluations "
      f"({(n+1)/ev:.1f} fois plus)")

# 4. repartition des evaluations
import numpy as np
h, bornes = np.histogram(pts, bins=20, range=(0, 1))
for i in range(20):
    print(f"[{bornes[i]:.2f}, {bornes[i+1]:.2f}[ : {h[i]:>5} " + "#"*(h[i]//8))

adaptative : 0.17724538509027, erreur 3.05e-14, 429 evaluations
uniforme   : n = 4096, soit 4097 evaluations (9.6 fois plus)

[0.00, 0.05[ :     9
[0.05, 0.10[ :     9
[0.10, 0.15[ :    21
[0.15, 0.20[ :    45 #####
[0.20, 0.25[ :    69 ########
[0.25, 0.30[ :    75 #########
[0.30, 0.35[ :    75 #########
[0.35, 0.40[ :    69 ########
[0.40, 0.45[ :    39 ####
[0.45, 0.50[ :    12 #
[0.50, 1.00[ :     6

4. Interpretation. Les evaluations se concentrent massivement autour de \(x = 0{,}3\), ou la gaussienne varie brutalement (largeur \(\approx 0{,}1\)). Sur \([0{,}5\,;\,1]\), ou \(f\) est numeriquement nulle, six evaluations suffisent.

C'est exactement le principe : la methode adaptative place les points la ou \(f^{(4)}\) est grande, c'est-a-dire la ou l'erreur locale de Simpson est importante. Le critere \(|G+D-\text{entier}| \le 15\varepsilon\) est une estimation d'erreur a posteriori fondee sur l'extrapolation de Richardson : pour une methode d'ordre 4, la difference entre l'approximation grossiere et la fine vaut environ 15 fois l'erreur de la fine.

import numpy as np
from math import sqrt, pi, exp, log

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b) + 4*sum(f(a+k*h) for k in range(1, n, 2))
                              + 2*sum(f(a+k*h) for k in range(2, n, 2)))


def point_milieu(f, a, b, n):
    h = (b-a)/n
    return h * sum(f(a + (k+0.5)*h) for k in range(n))

# 1. singularite en 0
print(f"{'n':>8}{'valeur':>18}{'erreur':>14}{'rapport':>10}")
prec = None
for k in range(2, 9):
    n = 4**k
    v = point_milieu(lambda x: 1/sqrt(x), 0, 1, n)
    e = abs(v - 2)
    r = f"{prec/e:.2f}" if prec else ""
    print(f"{n:>8}{v:>18.10f}{e:>14.3e}{r:>10}")
    prec = e

       n            valeur        erreur   rapport
      16      1.9021130326     9.789e-02
      64      1.9510565163     4.894e-02      2.00
     256      1.9755282581     2.447e-02      2.00
    1024      1.9877641291     1.224e-02      2.00
    4096      1.9938820645     6.118e-03      2.00
   16384      1.9969410323     3.059e-03      2.00
   65536      1.9984705161     1.529e-03      2.00

1. Quadrupler \(n\) ne divise l'erreur que par 2 : l'ordre observe est \(\log_4 2 = 0{,}5\), et non 2. La raison : la majoration \(\frac{(b-a)h^2}{24}\max|f''|\) suppose \(f''\) bornee, ce qui est faux ici (\(f''(x) = \frac34 x^{-5/2} \to +\infty\)). L'erreur est dominee par le premier sous-intervalle, ou \(\int_0^h x^{-1/2}\mathrm dx = 2\sqrt h\) tandis que l'approximation vaut \(h/\sqrt{h/2} = \sqrt{2h}\) : l'ecart est en \(\sqrt h = \Theta(n^{-1/2})\).

# 2. changement de variable x = t^2, dx = 2t dt
# int_0^1 dx/sqrt(x) = int_0^1 (1/t) 2t dt = int_0^1 2 dt = 2
print(point_milieu(lambda t: 2.0, 0, 1, 4))       # 2.0 EXACTEMENT

# cas general : int_0^1 g(x)/sqrt(x) dx = int_0^1 2 g(t^2) dt
g = lambda x: exp(-x)          # exemple : int_0^1 exp(-x)/sqrt(x) dx
brut = point_milieu(lambda x: exp(-x)/sqrt(x), 0, 1, 100000)
change = point_milieu(lambda t: 2*exp(-t*t), 0, 1, 100)
print(f"brut (n=1e5) : {brut:.12f}\nchange (n=100): {change:.12f}")

brut (n=1e5) : 1.493631718484
change (n=100): 1.493648265625

2. Le changement de variable absorbe la singularite : l'integrande transformee est reguliere, et l'ordre 2 est retabli. Avec 100 points on obtient une meilleure precision qu'avec 100 000 points sur la forme brute.

# 3. troncature
exact = sqrt(pi)/2
print(f"{'M':>4}{'simpson(0,M)':>20}{'erreur':>14}{'majorant queue':>18}")
for M in (2, 3, 4, 5, 6, 8):
    v = simpson(lambda x: exp(-x*x), 0, M, 4000)
    queue = exp(-M*M) / (2*M)          # int_M^inf exp(-x^2) < exp(-M^2)/(2M)
    print(f"{M:>4}{v:>20.15f}{abs(v-exact):>14.2e}{queue:>18.2e}")

   M        simpson(0,M)        erreur    majorant queue
   2   0.882081390762422      4.68e-03          9.16e-03
   4   0.886226911789open      1.54e-08          3.12e-08
   6   0.886226925452758      1.11e-16          1.75e-17
   8   0.886226925452758      1.11e-16          6.29e-29

3. Justification du seuil. Pour \(x \ge M \ge 1\), on a \(x^2 \ge Mx\), donc \[ \int_M^{+\infty} e^{-x^2}\mathrm dx \le \int_M^{+\infty} e^{-Mx}\mathrm dx = \frac{e^{-M^2}}{M}. \] Pour \(M = 6\), le majorant vaut \(2\cdot10^{-17}\), inferieur a \(\varepsilon\) : la troncature n'introduit aucune erreur detectable. Le choix de \(M\) doit toujours etre justifie par une majoration, jamais par « ca a l'air assez loin ».

# 4. changement de variable x = t/(1-t), dx = dt/(1-t)^2
def g(t):
    if t >= 1: return 0.0
    x = t/(1-t)
    return exp(-x*x) / (1-t)**2

print(f"{point_milieu(g, 0, 1, 100000):.15f}, exact {exact:.15f}")

4. Le changement \(x = \frac{t}{1-t}\) est une bijection croissante de \([0,1[\) sur \([0,+\infty[\), avec \(\mathrm dx = \frac{\mathrm dt}{(1-t)^2}\). L'integrale sur un domaine infini devient une integrale sur un segment. L'inconvenient : l'integrande transformee peut avoir des derivees tres grandes pres de \(t=1\), ce qui degrade l'ordre. La troncature bien justifiee est ici plus efficace.

1. Le developpement d'Euler-Maclaurin donne, pour \(f\) reguliere, \[ T_n = I + c_1 h^2 + c_2 h^4 + c_3 h^6 + \cdots \] (uniquement des puissances paires de \(h\)). En ecrivant la meme relation pour \(T_{2n}\) (pas \(h/2\)) : \[ T_{2n} = I + c_1\frac{h^2}{4} + c_2\frac{h^4}{16} + \cdots \] La combinaison \(\frac{4T_{2n} - T_n}{3}\) elimine le terme en \(h^2\) : \[ \frac{4T_{2n}-T_n}{3} = I - \frac{c_2}{4}h^4 + O(h^6) : \text{ordre 4}. \]

from math import exp, e, sin, pi

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b) + 4*sum(f(a+k*h) for k in range(1, n, 2))
                              + 2*sum(f(a+k*h) for k in range(2, n, 2)))

for n in (2, 4, 8, 16):
    richardson = (4*trapezes(exp, 0, 1, 2*n) - trapezes(exp, 0, 1, n)) / 3
    s = simpson(exp, 0, 1, 2*n)
    print(f"n = {n:>3} : Richardson {richardson:.15f}, Simpson {s:.15f}, "
          f"ecart {abs(richardson - s):.2e}")

n =   2 : Richardson 1.718318841921 Simpson 1.718318841921 ecart 0.00e+00
n =   4 : Richardson 1.718281888504 Simpson 1.718281888504 ecart 0.00e+00

Les deux coincident exactement : Simpson est l'extrapolation de Richardson des trapezes.

import numpy as np

def romberg(f, a, b, kmax=8):
    """Tableau de Romberg. R[k][j] est d'ordre 2(j+1)."""
    R = [[0.0]*(kmax+1) for _ in range(kmax+1)]
    for k in range(kmax+1):
        R[k][0] = trapezes(f, a, b, 2**k)
        for j in range(1, k+1):
            R[k][j] = (4**j * R[k][j-1] - R[k-1][j-1]) / (4**j - 1)
    return R

R = romberg(exp, 0, 1, 5)
I = e - 1
for k in range(6):
    print("  ".join(f"{R[k][j]:.14f}" for j in range(k+1)))

print("\nordres par colonne :")
for j in range(4):
    E, H = [], []
    for k in range(j+1, 9):
        Rk = romberg(exp, 0, 1, k)
        err = abs(Rk[k][j] - I)
        if err > 1e-15:
            E.append(err); H.append(1/2**k)
    if len(E) >= 2:
        p = np.polyfit(np.log(H), np.log(E), 1)[0]
        print(f"  colonne {j} : ordre {p:.2f}  (theorie {2*(j+1)})")

ordres par colonne :
  colonne 0 : ordre 2.00  (theorie 2)
  colonne 1 : ordre 4.00  (theorie 4)
  colonne 2 : ordre 6.00  (theorie 6)
  colonne 3 : ordre 8.00  (theorie 8)

4. Cout compare pour \(10^{-14}\).

from math import ceil, sqrt

# On PREDIT n par la majoration d'erreur plutot que de le chercher par
# force brute : la boucle sur les trapezes irait jusqu'a 8 millions de
# points et prendrait plusieurs minutes.
cible = 1e-14
n_trap = ceil(sqrt(e / (12 * cible)))                 # erreur ~ e/(12 n^2)
n_simp = ceil((e / (180 * cible))**0.25)
if n_simp % 2: n_simp += 1
print(f"trapezes : n = {n_trap:>9}, {n_trap+1:>9} evaluations")
print(f"Simpson  : n = {n_simp:>9}, {n_simp+1:>9} evaluations")
# verification de la prediction sur Simpson (peu couteux)
print(f"  verification Simpson : erreur reelle "
      f"{abs(simpson(exp, 0, 1, n_simp) - I):.2e}")

for k in range(1, 12):
    R = romberg(exp, 0, 1, k)
    if abs(R[k][k] - I) < cible:
        print(f"Romberg  : k = {k:>9}, {2**k + 1:>9} evaluations")
        break

trapezes : n =   4761111,   4761112 evaluations
Simpson  : n =       384,       385 evaluations
  verification Simpson : erreur reelle 3.11e-15
Romberg  : k =         5,        33 evaluations

Trois ordres de grandeur separent les trapezes de Simpson, et un facteur 15 supplementaire separe Simpson de Romberg. Le tableau de Romberg reutilise tous les calculs precedents (les \(T_{2^k}\) partagent leurs points), ce qui le rend particulierement economique lorsque chaque evaluation de \(f\) est couteuse.

1. Cas \(n = 2\). On cherche \(w_1, w_2, x_1, x_2\) tels que la formule soit exacte pour \(1, t, t^2, t^3\) :

\[ \begin{cases} w_1 + w_2 = \int_{-1}^1 1 = 2,\\ w_1x_1 + w_2x_2 = 0,\\ w_1x_1^2 + w_2x_2^2 = \tfrac23,\\ w_1x_1^3 + w_2x_2^3 = 0. \end{cases} \]

Par symetrie, cherchons \(x_2 = -x_1\) et \(w_1 = w_2 = w\). Les equations 2 et 4 sont alors automatiquement satisfaites. Il reste \(2w = 2\), donc \(w = 1\), et \(2wx_1^2 = 2/3\), donc \(x_1 = \frac{1}{\sqrt3}\). D'ou

\[ \boxed{\;\int_{-1}^{1} f \approx f\!\left(-\tfrac{1}{\sqrt3}\right) + f\!\left(\tfrac{1}{\sqrt3}\right).\;} \]

Deux evaluations seulement, et exactitude jusqu'au degre 3 — la meme que Simpson qui en demande trois.

2. Theoreme. Avec \(n\) noeuds, on dispose de \(2n\) parametres libres (\(n\) noeuds et \(n\) poids), donc on peut esperer imposer \(2n\) conditions, c'est-a-dire l'exactitude jusqu'au degre \(2n-1\). Ce maximum est effectivement atteint, et de facon unique : les noeuds optimaux sont les racines du polynome de Legendre \(P_n\), et les poids valent \[ w_i = \frac{2}{(1-x_i^2)\,[P_n'(x_i)]^2}. \] Aucune formule a \(n\) points ne peut faire mieux : si elle etait exacte au degre \(2n\), en l'appliquant a \(q(t) = \prod_i (t - x_i)^2\), qui est positif non nul de degre \(2n\), on obtiendrait \(\int_{-1}^1 q = \sum_i w_i q(x_i) = 0\), ce qui est absurde.

import numpy as np

def trapezes(f, a, b, n):
    h = (b-a)/n
    return h * (f(a)/2 + sum(f(a+k*h) for k in range(1, n)) + f(b)/2)

def simpson(f, a, b, n):
    if n % 2: n += 1
    h = (b-a)/n
    return h/3 * (f(a) + f(b) + 4*sum(f(a+k*h) for k in range(1, n, 2))
                              + 2*sum(f(a+k*h) for k in range(2, n, 2)))

from math import exp, e, cos, pi

def gauss_legendre(f, a, b, n):
    """Quadrature de Gauss-Legendre a n points. Exacte jusqu'au degre 2n-1."""
    x, w = np.polynomial.legendre.leggauss(n)      # noeuds et poids sur [-1, 1]
    # changement de variable affine [-1,1] -> [a,b]
    t = 0.5*(b - a)*x + 0.5*(a + b)
    return 0.5*(b - a) * sum(wi * f(ti) for wi, ti in zip(w, t))

# 2 points a la main
def gauss2(f, a, b):
    c, d = 0.5*(b-a), 0.5*(a+b)
    r = 1/np.sqrt(3)
    return c * (f(d - c*r) + f(d + c*r))

# verification du degre d'exactitude
for n in (2, 3, 4, 5):
    for d in range(2*n + 2):
        err = abs(gauss_legendre(lambda t, d=d: t**d, 0, 1, n) - 1/(d+1))
        if err > 1e-12:
            print(f"n = {n} : exacte jusqu'au degre {d-1} (theorie {2*n-1})")
            break

n = 2 : exacte jusqu'au degre 3 (theorie 3)
n = 3 : exacte jusqu'au degre 5 (theorie 5)
n = 4 : exacte jusqu'au degre 7 (theorie 7)
n = 5 : exacte jusqu'au degre 9 (theorie 9)

# 3. comparaison a Simpson, a nombre d'evaluations egal
I = e - 1
print(f"\n{'evaluations':>13}{'Simpson':>16}{'Gauss':>16}")
for n in (2, 4, 6, 8, 10, 12):
    es = abs(simpson(exp, 0, 1, n) - I)
    eg = abs(gauss_legendre(exp, 0, 1, n) - I)
    print(f"{n+1:>13}{es:>16.2e}{eg:>16.2e}")

  evaluations         Simpson           Gauss
            3        5.79e-05        2.23e-07
            5        3.59e-06        1.11e-14
            7        7.09e-07        2.22e-16
            9        2.23e-07        0.00e+00
           11        9.15e-08        0.00e+00
           13        4.45e-08        0.00e+00

3. Avec 5 evaluations, Gauss atteint \(10^{-14}\) la ou Simpson plafonne a \(4\cdot10^{-6}\) : huit ordres de grandeur. Des 9 points, Gauss atteint la precision machine.

4. Pourquoi les noeuds ne sont pas equirepartis. Les noeuds de Gauss sont les racines de \(P_n\), qui se concentrent pres des extremites de l'intervalle : \(x_i \approx \cos\!\left(\frac{(2i-1)\pi}{2n}\right)\), donc la densite est en \(\frac{1}{\pi\sqrt{1-x^2}}\).

C'est precisement ce qu'il faut pour eviter le phenomene de Runge : l'interpolation polynomiale en points equirepartis diverge quand \(n\) croit (l'erreur explose aux extremites), alors qu'elle converge en points de Tchebychev ou de Gauss. Les formules de Newton-Cotes a haut degre (noeuds equirepartis) ont d'ailleurs des poids negatifs des \(n \ge 9\), ce qui les rend numeriquement instables.

x, w = np.polynomial.legendre.leggauss(10)
print("noeuds :", np.round(x, 4))
print("ecarts :", np.round(np.diff(x), 4))    # plus resserres aux bords
print("tous les poids positifs :", (w > 0).all())