FEUILLE DE TP 3

Intégration numérique (Méthodes élémentaires)¶

Dans ce TP, nous programmons et nous étudions les méthodes élémentaires d'intégration numérique basées sur la méthode d'interpolation de Lagrange.

In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np                       # pour les numpy array
import matplotlib.pyplot as plt          # librairie graphique

color_bleu = 'xkcd:dusky blue'
color_vert = 'xkcd:grassy green'
color_orange = 'xkcd:orangish'
color_rouge = 'xkcd:blood orange'

Dans ce TP, nous allons programmer et étudier différentes formules de quadrature élémentaires. Nous les testerons sur les fonctions suivantes : \begin{align} f_0(x) &= \exp(-x), & f_1(x) &= \frac{\exp(-x^2)}{2},& f_2(x) &= \sqrt{\sin x},\\ %f_2(x) &= 0.1 \,(1+x/3)^3 \log(1+x),& %f_3(x) &= \sin(2\pix) / (1.1 - \sin(\pix)). \end{align}

Nous allons programmer les méthodes élémentaires autour du point 0. Notons (pour l'analyse) que $f_0'(0) = -1$, $f_1'(0) = 0$ et $f_2$ est non dérivable en 0.

Question 1.

Définissez ces 3 fonctions et tracez les sur l'intervalle $[0, 3]$.

In [2]:
def f_0(x):
    return np.exp(-x)

def f_1(x):
    return .5 * np.exp(-x*x)

def f_2(x):
    return np.sqrt(np.sin(x))

liste_f = [f_0, f_1, f_2]
In [13]:
x = np.linspace(0, 3, 200)
fig = plt.figure(figsize=(4*len(liste_f), 4))
fig.suptitle("Les fonctions que l'on va utiliser", fontsize=16, color=color_bleu)
for k, f in enumerate(liste_f):
    ax = fig.add_subplot(1, len(liste_f), k+1)
    ax.plot(x, f(x), color=color_orange, linewidth=2)
    ax.set_title(f"${f.__name__}$")

Intégration en utilisant le module scipy¶

La fonction quad de scipy.integrate permet de calculer la valeur approchée d'une intégrale à une précision donnée. Lisez attentivement la documentation https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html.

Question 2.

Calculer les intégrales de $f_i$, $ i = 0, 1, 2 $, sur $[0, 0.1]$ en utilisant la fonction quad de scipy.integrate. On choissira une erreur absolue epsabs=1e-14. Pour $ i = 0 $ comparer la valeur obtenue avec la valeur de l'intégrale calculée à la main.

In [4]:
from scipy.integrate import quad
x = 0.1

liste_Ie = [quad(f, 0, x, epsabs=1.e-14)[0] for f in liste_f]
I0, I1, I2 = liste_Ie
print(I0, I0 - (1.-np.exp(-x)))
print(I1)
print(I2)

0.09516258196404043 -5.551115123125783e-17
0.04983383214516818
0.0210743222343681

Formule du rectangle¶

La formule du rectangle consiste à faire l'approximation

\begin{equation} \int_{a}^{b} f \approx J^R( f ) = (b-a) f(a). \end{equation}

Question 3.

  1. Programmez une fonction quad_rect(f,a,b) qui prend en argument une fonction f et deux réels a et b, et qui retourne la valeur approchée de $ \int_a^b f(x)dx $ obtenue par la formule du rectangle.
  2. Calculez les valeurs approchées des intégrales de $ f_i $, $ i = 0, 1, 2 $, sur $ [ 0, 0.1 ] $, comparez les valeurs obtenues avec celles de la question précédente.
In [5]:
def quad_rect(f, a, b):
    """
    méthode élémentaire du rectangle à gauche
    
    Parameters
    ----------
    
    f: function
        la fonction que l'on intègre
    a: float
        le bord gauche de l'intégrale
    b: float
        le bord droit de l'intégrale
    """
    return (b - a) * f(a)
In [6]:
for f, Ie in zip(liste_f, liste_Ie):
    Ia = quad_rect(f, 0, x)
    print(f"Ia = {Ia:9.7f}, Ie = {Ie:9.7f}, E = {abs(Ie-Ia):9.3e}")

Ia = 0.1000000, Ie = 0.0951626, E = 4.837e-03
Ia = 0.0500000, Ie = 0.0498338, E = 1.662e-04
Ia = 0.0000000, Ie = 0.0210743, E = 2.107e-02

Etude du comportement d'erreur¶

On supposera que $ a = 0 $ et $ b = h $. On utilisera la notation $ J^R_h( f ) = J^R( f ) $ dans ce cas.

Question 4.

  1. Calculez les valeurs approchée de $I_h(f_0) = \int_{0}^{h}f_0(x)dx $ à l'aide de la formule du rectangle pour $h=\lbrace 1, 1/2, 1/4, \ldots, 2^{-10}\rbrace$.
  2. Affichez alors en respectant le format ci-dessous au maximum les valeurs de $h$, $E_h(f_0)$, $E_h(f_0)/h$ et $E_h(f_0)/h^2$, où $E_h(f_0)=\vert I_h( f_0 )-J^R_h( f_0 )\vert$ :
-----------------------------------------------
  h         Eh(f0)/h^0  Eh(f0)/h^1  Eh(f0)/h^2
-----------------------------------------------
 1.000e+00  3.6788e-01  3.6788e-01  3.6788e-01
 5.000e-01  1.0653e-01  2.1306e-01  4.2612e-01
 2.500e-01  2.8801e-02  1.1520e-01  4.6081e-01
...
-----------------------------------------------

Rappel : voici un exemple pour écrire le nombre 0.14159 sous la forme 1.4159e-01

print(f"{0.14159:11.4e}")
  1. Comparez la convergence de $ E_h( f_0 ) $ vers $ 0 $ à celle d'une puissance de $ h $ vers 0. Les résultats obtenus sont-ils conformes aux résultats vus en cours?
  2. Recommencez avec $f_1$ et $f_2$. Que remarquez-vous ? Pouvez-vous l'expliquer ?
In [7]:
hs = 1/2.**np.arange(0, 11)

lst_puiss = [1, 1.5, 2, 3]
ndigits = 7
tc = ndigits+8
taille = (1+len(lst_puiss))*tc

for f in liste_f:
    print('*'*taille)
    print(f"Comportement pour {f.__name__}")
    print('-'*taille)
    print(f"{'h':^{tc}}", end='')
    for p in lst_puiss:
        print(f"{f'Eh({f.__name__})/h^{p}':^{tc}}", end='')
    print()
    print('-'*taille)
    for h in hs:
        Iex = quad(f, 0, h, epsabs = 1e-14)[0]
        Iappr = quad_rect(f, 0, h)
        E = np.abs(Iex - Iappr)
        msg = f"{h:^{ndigits+6}.{ndigits}e}"
        print(f"{msg:^{tc}}", end='')
        for p in lst_puiss:
            msg = f"{E/h**p:^{ndigits+6}.{ndigits}e}"
            print(f"{msg:^{tc}}", end='')
        print()
    print('-'*taille)

***************************************************************************
Comportement pour f_0
---------------------------------------------------------------------------
       h         Eh(f_0)/h^1   Eh(f_0)/h^1.5   Eh(f_0)/h^2    Eh(f_0)/h^3  
---------------------------------------------------------------------------
 1.0000000e+00  3.6787944e-01  3.6787944e-01  3.6787944e-01  3.6787944e-01 
 5.0000000e-01  2.1306132e-01  3.0131421e-01  4.2612264e-01  8.5224528e-01 
 2.5000000e-01  1.1520313e-01  2.3040626e-01  4.6081253e-01  1.8432501e+00 
 1.2500000e-01  5.9975221e-02  1.6963554e-01  4.7980177e-01  3.8384141e+00 
 6.2500000e-02  3.0609005e-02  1.2243602e-01  4.8974408e-01  7.8359053e+00 
 3.1250000e-02  1.5463503e-02  8.7474784e-02  4.9483210e-01  1.5834627e+01 
 1.5625000e-02  7.7719683e-03  6.2175747e-02  4.9740597e-01  3.1833982e+01 
 7.8125000e-03  3.8960973e-03  4.4079309e-02  4.9870046e-01  6.3833658e+01 
 3.9062500e-03  1.9505844e-03  3.1209350e-02  4.9934959e-01  1.2783350e+02 
 1.9531250e-03  9.7592703e-04  2.2082708e-02  4.9967464e-01  2.5583341e+02 
 9.7656250e-04  4.8812234e-04  1.5619915e-02  4.9983728e-01  5.1183337e+02 
---------------------------------------------------------------------------
***************************************************************************
Comportement pour f_1
---------------------------------------------------------------------------
       h         Eh(f_1)/h^1   Eh(f_1)/h^1.5   Eh(f_1)/h^2    Eh(f_1)/h^3  
---------------------------------------------------------------------------
 1.0000000e+00  1.2658793e-01  1.2658793e-01  1.2658793e-01  1.2658793e-01 
 5.0000000e-01  3.8718994e-02  5.4756926e-02  7.7437987e-02  1.5487597e-01 
 2.5000000e-01  1.0224226e-02  2.0448451e-02  4.0896903e-02  1.6358761e-01 
 1.2500000e-01  2.5920049e-03  7.3312970e-03  2.0736039e-02  1.6588831e-01 
 6.2500000e-02  6.5027944e-04  2.6011177e-03  1.0404471e-02  1.6647154e-01 
 3.1250000e-02  1.6271274e-04  9.2044228e-04  5.2068078e-03  1.6661785e-01 
 1.5625000e-02  4.0687124e-05  3.2549699e-04  2.6039759e-03  1.6665446e-01 
 7.8125000e-03  1.0172340e-05  1.1508689e-04  1.3020595e-03  1.6666361e-01 
 3.9062500e-03  2.5431199e-06  4.0689918e-05  6.5103869e-04  1.6666590e-01 
 1.9531250e-03  6.3578215e-07  1.4386108e-05  3.2552046e-04  1.6666648e-01 
 9.7656250e-04  1.5894567e-07  5.0862616e-06  1.6276037e-04  1.6666662e-01 
---------------------------------------------------------------------------
***************************************************************************
Comportement pour f_2
---------------------------------------------------------------------------
       h         Eh(f_2)/h^1   Eh(f_2)/h^1.5   Eh(f_2)/h^2    Eh(f_2)/h^3  
---------------------------------------------------------------------------
 1.0000000e+00  6.4297763e-01  6.4297763e-01  6.4297763e-01  6.4297763e-01 
 5.0000000e-01  4.6720107e-01  6.6072209e-01  9.3440214e-01  1.8688043e+00 
 2.5000000e-01  3.3258953e-01  6.6517906e-01  1.3303581e+00  5.3214325e+00 
 1.2500000e-01  2.3557074e-01  6.6629467e-01  1.8845659e+00  1.5076527e+01 
 6.2500000e-02  1.6664342e-01  6.6657366e-01  2.6662947e+00  4.2660714e+01 
 3.1250000e-02  1.1784702e-01  6.6664342e-01  3.7711046e+00  1.2067535e+02 
 1.5625000e-02  8.3332607e-02  6.6666085e-01  5.3332868e+00  3.4133036e+02 
 7.8125000e-03  5.8925437e-02  6.6666521e-01  7.5424559e+00  9.6543435e+02 
 3.9062500e-03  4.1666644e-02  6.6666630e-01  1.0666661e+01  2.7306652e+03 
 1.9531250e-03  2.9462779e-02  6.6666658e-01  1.5084943e+01  7.7234906e+03 
 9.7656250e-04  2.0833333e-02  6.6666664e-01  2.1333333e+01  2.1845333e+04 
---------------------------------------------------------------------------

Formules de Newton-Cotes¶

Supposons que les points $ x_1, x_2, \ldots, x_N $ sont définies par l'une des formules suivantes: $$ x_j=a+(j-1)\frac{b-a}{N-1}, \quad j = 1, \ldots, N, \quad (1) $$ ou $$ x_j=a+j\frac{b-a}{N+1}, \quad j = 1, \ldots, N. \quad (2) $$

On rappelle que la formule de Newton-Cotes à $N$ points consiste à approcher $\int_a^b f(x)\,dx$ par
$$ J^{NC}( f ) = \sum\limits_{k=1}^{N} \lambda_k f(x_k) \quad \text{ avec } \quad \lambda_k= \int_a^b L_i( x ) dx\,, $$ où $ L_i $ est le $i$-ème polynome de Lagrange associé aux points $ x_1, x_2, \ldots, x_N $.

Si les points $ x_1, \ldots, x_N $ sont donnés par la relation (1), on parle de la formule de Newton-Cotes fermée; s'ils sont donnés par la relation (2), il s'agit de la formule de Newton-Cotes ouverte.

Question 5.

  1. Proposez une fonction points_NC_fermee(a, b, N) qui prend en arguments 2 réels, les valeurs des bornes de l'intervalle a et b et 1 entier N et qui retourne un ndarray contenant les points d'interpolation pour la formule de Newton-Cotes fermée.
  2. Proposez une fonction points_NC_ouverte(a, b, N) qui prend en arguments 2 réels, les valeurs des bornes de l'intervalle a et b et 1 entier N et qui retourne un ndarray contenant les points d'interpolation pour la formule de Newton-Cotes ouverte.
  3. Proposez une fonction quad_NC(f, N, a, b, lam, ouv) qui prend en arguments une fonction f, un entier N, deux réels a, b, un ndarray lam et un booléen ouv et qui retourne la valeur approchée de $ \int_a^b f(x)dx $ obtenue par la formule de Newton-Cotes à $N$ points avec les coefficients lam, ouverte si ouv vaut True et fermée si ouv est False.
In [8]:
def points_NC_fermee(a, b, N):
    """retourne les points de Newton-Cotes fermés"""
    return np.linspace(a, b, N)

def points_NC_ouverte(a, b, N):
    """retourne les points de Newton-Cotes ouverts"""
    return np.linspace(a, b, N+2)[1:-1]

def quad_NC(f, N, a, b, lam, ouv):
    """
    Formule de quadrature de Newton-Cotes
    
    Parameters
    ----------
    
    f: function
        fonction que l'on souhaite intégrer
    N: int
        nombre de points pour l'interpolation
    a: float
        borne de gauche pour l'intégrale
    b: float
        borne de droite pour l'intégrale
    lam: ndarray
        poids de la formule de quadrature
    ouv: bool
        formule ouverte ou fermée
    """
    assert len(lam) == N, "Erreur dans quad_NC: lam n'a pas la bonne taille"
    if ouv:
        x = points_NC_ouverte(a, b, N)
    else:
        x = points_NC_fermee(a, b, N)
    return np.sum(lam * f(x))

Cas particuliers: formules du point milieu, trapèze et Simpson¶

Question 6.

  1. Rappelez les valeurs des coefficients $ \lambda_k $ pour les formules du point milieu, trapèze et Simpson (toutes les trois sont des cas particuliers de la formule de Newton-Cotes).
  2. A l'aide de la fonction quad_NC programmez les fonctions quad_milieu, quad_trapeze, quad_simpson qui prend en argument une fonction f et deux réels a et b et qui retournent la valeur approchée de $ \int_a^b f(x)dx $ obtenue par la formule du point milieu, du trapèze, de Simpson respectivement.
  3. Testez vos fonctions en calculant une valeur approchée de l'intégrale de $ f_i $, $ i = 0, 1, 2 $; comparez avec la valeur de l'intégrale obtenue à l'aide de la fonction quad.
In [9]:
def quad_milieu(f, a, b):
    """méthode élémentaire du point milieu"""
    N = 1
    lam_milieu = (b - a) * np.array([1])
    return quad_NC(f, N, a, b, lam_milieu, True)

def quad_trapeze(f, a, b):
    """méthode élémentaire du trapèze"""
    N = 2
    lam_trapeze = (b - a) * np.array([0.5, 0.5])
    return quad_NC(f, N, a, b, lam_trapeze, False)

def quad_simpson(f, a, b):
    """méthode élémentaire de Simpson"""
    N = 3
    lam_simpson = (b - a) * np.array([1/6, 4/6, 1/6])
    return quad_NC(f, N, a, b, lam_simpson, False)

h = 0.1
for meth, nom in zip(
    [quad_milieu, quad_trapeze, quad_simpson],
    ["Point milieu", "Trapèze", "Simpson"]
):
    print("-"*80)
    print(f"{nom:^80s}")
    print("-"*80)
    for f in liste_f:
        Iapp = meth(f, 0, h)
        Iex = quad(f, 0, h, epsabs = 1e-14)[0]
        print(f"{f.__name__} : J = {Iapp:9.7f},   I = {Iex:9.7f},   erreur = {abs(Iapp-Iex):9.3e}")

--------------------------------------------------------------------------------
                                  Point milieu                                  
--------------------------------------------------------------------------------
f_0 : J = 0.0951229,   I = 0.0951626,   erreur = 3.964e-05
f_1 : J = 0.0498752,   I = 0.0498338,   erreur = 4.132e-05
f_2 : J = 0.0223560,   I = 0.0210743,   erreur = 1.282e-03
--------------------------------------------------------------------------------
                                    Trapèze                                     
--------------------------------------------------------------------------------
f_0 : J = 0.0952419,   I = 0.0951626,   erreur = 7.929e-05
f_1 : J = 0.0497512,   I = 0.0498338,   erreur = 8.259e-05
f_2 : J = 0.0157982,   I = 0.0210743,   erreur = 5.276e-03
--------------------------------------------------------------------------------
                                    Simpson                                     
--------------------------------------------------------------------------------
f_0 : J = 0.0951626,   I = 0.0951626,   erreur = 3.303e-09
f_1 : J = 0.0498339,   I = 0.0498338,   erreur = 2.055e-08
f_2 : J = 0.0201701,   I = 0.0210743,   erreur = 9.042e-04
  1. Reprenez la question 4 pour chacune de ces trois méthodes pour étudier le comportement de l'erreur d'intégration en fonction de la longuer de l'intervalle d'intégration $ h $ (on étudiera $ E_h(f_0)/h^r $ pour quelques valeurs $ r \in \mathbb{N} $).
In [10]:
hs = 1/2.**np.arange(0, 11)

lst_puiss = [1.5, 2, 3, 5, 6]
ndigits = 5
tc = ndigits+8
taille = (1+len(lst_puiss))*tc

for meth, nom in zip(
    [quad_milieu, quad_trapeze, quad_simpson],
    ["Point milieu", "Trapèze", "Simpson"]
):
    print("*"*taille)
    print(f"{nom:^{taille}}")
    print("*"*taille)
    for f in liste_f:
        print('='*taille)
        print(f"Comportement pour {f.__name__} (méthode: {nom})")
        print('-'*taille)
        print(f"{'h':^{tc}}", end='')
        for p in lst_puiss:
            print(f"{f'Eh({f.__name__})/h^{p}':^{tc}}", end='')
        print()
        print('-'*taille)
        for h in hs:
            Iex = quad(f, 0, h, epsabs = 1e-14)[0]
            Iappr = meth(f, 0, h)
            E = np.abs(Iex - Iappr)
            msg = f"{h:^{ndigits+6}.{ndigits}e}"
            print(f"{msg:^{tc}}", end='')
            for p in lst_puiss:
                msg = f"{E/h**p:^{ndigits+6}.{ndigits}e}"
                print(f"{msg:^{tc}}", end='')
            print()
        print('-'*taille)

******************************************************************************
                                 Point milieu                                 
******************************************************************************
==============================================================================
Comportement pour f_0 (méthode: Point milieu)
------------------------------------------------------------------------------
      h      Eh(f_0)/h^1.5 Eh(f_0)/h^2  Eh(f_0)/h^3  Eh(f_0)/h^5  Eh(f_0)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  2.55899e-02  2.55899e-02  2.55899e-02  2.55899e-02  2.55899e-02 
 5.00000e-01  1.15087e-02  1.62758e-02  3.25516e-02  1.30206e-01  2.60413e-01 
 2.50000e-01  4.59993e-03  9.19986e-03  3.67994e-02  5.88791e-01  2.35516e+00 
 1.25000e-01  1.73020e-03  4.89373e-03  3.91499e-02  2.50559e+00  2.00447e+01 
 6.25000e-02  6.31042e-04  2.52417e-03  4.03867e-02  1.03390e+01  1.65424e+02 
 3.12500e-02  2.26612e-04  1.28191e-03  4.10212e-02  4.20057e+01  1.34418e+03 
 1.56250e-02  8.07471e-05  6.45977e-04  4.13425e-02  1.69339e+02  1.08377e+04 
 7.81250e-03  2.86601e-05  3.24252e-04  4.15043e-02  6.80006e+02  8.70407e+04 
 3.90625e-03  1.01527e-05  1.62443e-04  4.15854e-02  2.72534e+03  6.97687e+05 
 1.95312e-03  3.59302e-06  8.13008e-05  4.16260e-02  1.09120e+04  5.58695e+06 
 9.76562e-04  1.27095e-06  4.06702e-05  4.16463e-02  4.36693e+04  4.47174e+07 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_1 (méthode: Point milieu)
------------------------------------------------------------------------------
      h      Eh(f_1)/h^1.5 Eh(f_1)/h^2  Eh(f_1)/h^3  Eh(f_1)/h^5  Eh(f_1)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  1.59883e-02  1.59883e-02  1.59883e-02  1.59883e-02  1.59883e-02 
 5.00000e-01  1.19155e-02  1.68510e-02  3.37021e-02  1.34808e-01  2.69617e-01 
 2.50000e-01  4.94489e-03  9.88978e-03  3.95591e-02  6.32946e-01  2.53178e+00 
 1.25000e-01  1.81780e-03  5.14152e-03  4.11321e-02  2.63246e+00  2.10597e+01 
 6.25000e-02  6.48946e-04  2.59578e-03  4.15326e-02  1.06323e+01  1.70117e+02 
 3.12500e-02  2.29993e-04  1.30103e-03  4.16331e-02  4.26323e+01  1.36423e+03 
 1.56250e-02  8.13638e-05  6.50911e-04  4.16583e-02  1.70632e+02  1.09205e+04 
 7.81250e-03  2.87708e-05  3.25504e-04  4.16646e-02  6.82632e+02  8.73769e+04 
 3.90625e-03  1.01724e-05  1.62758e-04  4.16661e-02  2.73063e+03  6.99042e+05 
 1.95312e-03  3.59652e-06  8.13800e-05  4.16665e-02  1.09226e+04  5.59239e+06 
 9.76562e-04  1.27156e-06  4.06901e-05  4.16666e-02  4.36906e+04  4.47392e+07 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_2 (méthode: Point milieu)
------------------------------------------------------------------------------
      h      Eh(f_2)/h^1.5 Eh(f_2)/h^2  Eh(f_2)/h^3  Eh(f_2)/h^5  Eh(f_2)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  4.94280e-02  4.94280e-02  4.94280e-02  4.94280e-02  4.94280e-02 
 5.00000e-01  4.27038e-02  6.03922e-02  1.20784e-01  4.83138e-01  9.66276e-01 
 2.50000e-01  4.10071e-02  8.20143e-02  3.28057e-01  5.24891e+00  2.09956e+01 
 1.25000e-01  4.05819e-02  1.14783e-01  9.18264e-01  5.87689e+01  4.70151e+02 
 6.25000e-02  4.04756e-02  1.61902e-01  2.59044e+00  6.63152e+02  1.06104e+04 
 3.12500e-02  4.04490e-02  2.28814e-01  7.32205e+00  7.49778e+03  2.39929e+05 
 1.56250e-02  4.04423e-02  3.23539e-01  2.07065e+01  8.48137e+04  5.42808e+06 
 7.81250e-03  4.04407e-02  4.57534e-01  5.85643e+01  9.59518e+05  1.22818e+08 
 3.90625e-03  4.04403e-02  6.47044e-01  1.65643e+02  1.08556e+07  2.77903e+09 
 1.95312e-03  4.04401e-02  9.15056e-01  4.68509e+02  1.22817e+08  6.28822e+10 
 9.76562e-04  4.04401e-02  1.29408e+00  1.32514e+03  1.38951e+09  1.42286e+12 
------------------------------------------------------------------------------
******************************************************************************
                                   Trapèze                                    
******************************************************************************
==============================================================================
Comportement pour f_0 (méthode: Trapèze)
------------------------------------------------------------------------------
      h      Eh(f_0)/h^1.5 Eh(f_0)/h^2  Eh(f_0)/h^3  Eh(f_0)/h^5  Eh(f_0)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  5.18192e-02  5.18192e-02  5.18192e-02  5.18192e-02  5.18192e-02 
 5.00000e-01  2.30894e-02  3.26533e-02  6.53066e-02  2.61226e-01  5.22453e-01 
 2.50000e-01  9.20705e-03  1.84141e-02  7.36564e-02  1.17850e+00  4.71401e+00 
 1.25000e-01  3.46107e-03  9.78938e-03  7.83150e-02  5.01216e+00  4.00973e+01 
 6.25000e-02  1.26215e-03  5.04858e-03  8.07773e-02  2.06790e+01  3.30864e+02 
 3.12500e-02  4.53230e-04  2.56386e-03  8.20434e-02  8.40124e+01  2.68840e+03 
 1.56250e-02  1.61495e-04  1.29196e-03  8.26853e-02  3.38679e+02  2.16755e+04 
 7.81250e-03  5.73202e-05  6.48504e-04  8.30086e-02  1.36001e+03  1.74082e+05 
 3.90625e-03  2.03054e-05  3.24886e-04  8.31708e-02  5.45068e+03  1.39537e+06 
 1.95312e-03  7.18604e-06  1.62602e-04  8.32520e-02  2.18240e+04  1.11739e+07 
 9.76562e-04  2.54189e-06  8.13405e-05  8.32927e-02  8.73387e+04  8.94348e+07 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_1 (méthode: Trapèze)
------------------------------------------------------------------------------
      h      Eh(f_1)/h^1.5 Eh(f_1)/h^2  Eh(f_1)/h^3  Eh(f_1)/h^5  Eh(f_1)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  3.14422e-02  3.14422e-02  3.14422e-02  3.14422e-02  3.14422e-02 
 5.00000e-01  2.34488e-02  3.31616e-02  6.63232e-02  2.65293e-01  5.30586e-01 
 2.50000e-01  9.84502e-03  1.96900e-02  7.87601e-02  1.26016e+00  5.04065e+00 
 1.25000e-01  3.63138e-03  1.02711e-02  8.21687e-02  5.25880e+00  4.20704e+01 
 6.25000e-02  1.29751e-03  5.19005e-03  8.30408e-02  2.12584e+01  3.40135e+02 
 3.12500e-02  4.59952e-04  2.60188e-03  8.32601e-02  8.52584e+01  2.72827e+03 
 1.56250e-02  1.62725e-04  1.30180e-03  8.33150e-02  3.41258e+02  2.18405e+04 
 7.81250e-03  5.75413e-05  6.51006e-04  8.33288e-02  1.36526e+03  1.74753e+05 
 3.90625e-03  2.03448e-05  3.25516e-04  8.33322e-02  5.46126e+03  1.39808e+06 
 1.95312e-03  7.19304e-06  1.62760e-04  8.33330e-02  2.18453e+04  1.11848e+07 
 9.76562e-04  2.54313e-06  8.13801e-05  8.33333e-02  8.73813e+04  8.94784e+07 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_2 (méthode: Trapèze)
------------------------------------------------------------------------------
      h      Eh(f_2)/h^1.5 Eh(f_2)/h^2  Eh(f_2)/h^3  Eh(f_2)/h^5  Eh(f_2)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  1.84319e-01  1.84319e-01  1.84319e-01  1.84319e-01  1.84319e-01 
 5.00000e-01  1.71117e-01  2.41997e-01  4.83993e-01  1.93597e+00  3.87194e+00 
 2.50000e-01  1.67782e-01  3.35564e-01  1.34226e+00  2.14761e+01  8.59043e+01 
 1.25000e-01  1.66946e-01  4.72194e-01  3.77755e+00  2.41763e+02  1.93410e+03 
 6.25000e-02  1.66736e-01  6.66946e-01  1.06711e+01  2.73181e+03  4.37090e+04 
 3.12500e-02  1.66684e-01  9.42908e-01  3.01730e+01  3.08972e+04  9.88710e+05 
 1.56250e-02  1.66671e-01  1.33337e+00  8.53356e+01  3.49534e+05  2.23702e+07 
 7.81250e-03  1.66668e-01  1.88563e+00  2.41361e+02  3.95445e+06  5.06170e+08 
 3.90625e-03  1.66667e-01  2.66667e+00  6.82668e+02  4.47393e+07  1.14533e+10 
 1.95312e-03  1.66667e-01  3.77124e+00  1.93087e+03  5.06167e+08  2.59157e+11 
 9.76562e-04  1.66667e-01  5.33333e+00  5.46133e+03  5.72662e+09  5.86406e+12 
------------------------------------------------------------------------------
******************************************************************************
                                   Simpson                                    
******************************************************************************
==============================================================================
Comportement pour f_0 (méthode: Simpson)
------------------------------------------------------------------------------
      h      Eh(f_0)/h^1.5 Eh(f_0)/h^2  Eh(f_0)/h^3  Eh(f_0)/h^5  Eh(f_0)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  2.13121e-04  2.13121e-04  2.13121e-04  2.13121e-04  2.13121e-04 
 5.00000e-01  2.39729e-05  3.39028e-05  6.78057e-05  2.71223e-04  5.42446e-04 
 2.50000e-01  2.39571e-06  4.79142e-06  1.91657e-05  3.06651e-04  1.22660e-03 
 1.25000e-01  2.25284e-07  6.37199e-07  5.09759e-06  3.26246e-04  2.60997e-03 
 6.25000e-02  2.05417e-08  8.21667e-08  1.31467e-06  3.36555e-04  5.38488e-03 
 3.12500e-02  1.84417e-09  1.04322e-08  3.33831e-07  3.41843e-04  1.09390e-02 
 1.56250e-02  1.64278e-10  1.31423e-09  8.41105e-08  3.44517e-04  2.20491e-02 
 7.81250e-03  1.45755e-11  1.64903e-10  2.11076e-08  3.45826e-04  4.42657e-02 
 3.90625e-03  1.28964e-12  2.06342e-11  5.28235e-09  3.46184e-04  8.86230e-02 
 1.95312e-03  1.13047e-13  2.55795e-12  1.30967e-09  3.43323e-04  1.75781e-01 
 9.76562e-04  3.55271e-15  1.13687e-13  1.16415e-10  1.22070e-04  1.25000e-01 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_1 (méthode: Simpson)
------------------------------------------------------------------------------
      h      Eh(f_1)/h^1.5 Eh(f_1)/h^2  Eh(f_1)/h^3  Eh(f_1)/h^5  Eh(f_1)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  1.78148e-04  1.78148e-04  1.78148e-04  1.78148e-04  1.78148e-04 
 5.00000e-01  1.27392e-04  1.80160e-04  3.60319e-04  1.44128e-03  2.88255e-03 
 2.50000e-01  1.49198e-05  2.98395e-05  1.19358e-04  1.90973e-03  7.63891e-03 
 1.25000e-01  1.40804e-06  3.98254e-06  3.18604e-05  2.03906e-03  1.63125e-02 
 6.25000e-02  1.26478e-07  5.05911e-07  8.09457e-06  2.07221e-03  3.31554e-02 
 3.12500e-02  1.12241e-08  6.34933e-08  2.03179e-06  2.08055e-03  6.65776e-02 
 1.56250e-02  9.93079e-10  7.94463e-09  5.08456e-07  2.08264e-03  1.33289e-01 
 7.81250e-03  8.77983e-11  9.93325e-10  1.27146e-07  2.08315e-03  2.66644e-01 
 3.90625e-03  7.76090e-12  1.24174e-10  3.17887e-08  2.08330e-03  5.33325e-01 
 1.95312e-03  6.85816e-13  1.55183e-11  7.94535e-09  2.08282e-03  1.06641e+00 
 9.76562e-04  6.03961e-14  1.93268e-12  1.97906e-09  2.07520e-03  2.12500e+00 
------------------------------------------------------------------------------
==============================================================================
Comportement pour f_2 (méthode: Simpson)
------------------------------------------------------------------------------
      h      Eh(f_2)/h^1.5 Eh(f_2)/h^2  Eh(f_2)/h^3  Eh(f_2)/h^5  Eh(f_2)/h^6 
------------------------------------------------------------------------------
 1.00000e+00  2.84877e-02  2.84877e-02  2.84877e-02  2.84877e-02  2.84877e-02 
 5.00000e-01  2.85700e-02  4.04040e-02  8.08080e-02  3.23232e-01  6.46464e-01 
 2.50000e-01  2.85892e-02  5.71784e-02  2.28714e-01  3.65942e+00  1.46377e+01 
 1.25000e-01  2.85939e-02  8.08758e-02  6.47007e-01  4.14084e+01  3.31267e+02 
 6.25000e-02  2.85951e-02  1.14380e-01  1.83009e+00  4.68502e+02  7.49603e+03 
 3.12500e-02  2.85954e-02  1.61760e-01  5.17632e+00  5.30055e+03  1.69618e+05 
 1.56250e-02  2.85955e-02  2.28764e-01  1.46409e+01  5.99690e+04  3.83802e+06 
 7.81250e-03  2.85955e-02  3.23521e-01  4.14107e+01  6.78472e+05  8.68445e+07 
 3.90625e-03  2.85955e-02  4.57528e-01  1.17127e+02  7.67604e+06  1.96507e+09 
 1.95312e-03  2.85955e-02  6.47042e-01  3.31285e+02  8.68445e+07  4.44644e+10 
 9.76562e-04  2.85955e-02  9.15055e-01  9.37017e+02  9.82533e+08  1.00611e+12 
------------------------------------------------------------------------------

Méthodes d'ordre plus élevé¶

On souhaite obtenir les formules de Newton-Cotes à $ N $ points avec $ N = 4 $ et $ N = 5 $.

Question 7.

  1. Proposez une fonction Li(i, x) qui prend en arguments un entier i et un ndarray x et qui retourne la fonction $ f: x \mapsto L_i( x ) $, où $ L_i $ est le $i$-ème polynome de Lagrange associé aux points x (vous pouvez utiliser des fonctions programmées en TP06).
  2. Proposez une fonction NC_poids qui prend en arguments un entier N, 2 réels a et b et un ndarray x et qui retourne un ndarray contenant les poids $ \lambda_k $ de la formule de Newton-Cotes à $ N $ points x pour intégration sur $[a,b]$. On utilisera la fonction Li et la fonction quad du module scipy.integrate calculer les coéficients $ \lambda_k $ dans la formule de Newton-Cotes
  3. Testez votre code en calculant les poids $ \lambda_k $ dans la formule de Newton-Cotes fermée à $ N = 2 $ et $ N = 3 $. Obtenez-vous les poids de la formule du trapèze et de Simpson?
  4. Calculer les poids $ \lambda_k $ dans la formule de Newton-Cotes fermée à $ N = 4 $ et $ N = 5 $ points.
  5. Pouvez-vous deviner pourquoi la formule de Newton-Cotes à $ N = 4 $ points fermée est appelée "méthode de 3/8"?
In [11]:
def Li(i, x):
    """
    retourne une fonction qui évalue le polynôme Li construit avec les points x
    
    Parameters
    ----------
    
    i: int
        le numéro du polynôme de Lagrange
    x: ndarray
        les points de construction du polynôme Li
        
    Returns
    -------
    
    Li: function
        la fonction t -> Li(t) définie par
        Li(t) = \prod_{j\neq i} (t-x_j)/(x_i-x_j)
        
    Warning
    -------
    
    l'indice i doit être compris entre 0 (inclus) et taille de x (exclus)
    
    """
    n = x.size
    if i < 0 or i >= n:
        msg = "Problème dans le calcul de Li\n"
        msg += f"x.size = {n} et i = {i}"
        raise ValueError(msg)
    def P(t):
        Pt = 1.
        for j in range(i):
            Pt *= (t-x[j])/(x[i]-x[j])
        for j in range(i+1, x.size):
            Pt *= (t-x[j])/(x[i]-x[j])
        return Pt
    return P

def NC_poids(a, b, x):
    """
    retourne les poids de la formule d'intégration de Newton-Cotes 
    associée aux points du vecteur x
    sur l'intervalle [a, b]
    
    Parameters
    ----------
    
    a: float
        le bord gauche de l'intégrale
    b: float
        le bord droit de l'intégrale
    x: ndarray
        le vecteur des points d'intégration
        
    Returns
    -------
    
    lamda: ndarray
        le vecteur des poids de la formule d'intégration
        lamda_i = \int_a^b Li(t) dt
        
    Remarks
    -------
    
    La formule d'intégration associée s'écrira
    
    \int_a^b f(t) dt \sim \sum_{i=0}^{x.size-1} lamda_i f(x_i).
    """
    lamda = np.zeros(x.shape)
    for i in range(x.size):
        lamda[i], _ = quad(Li(i, x), a, b, epsabs=1.e-14)
    return lamda
In [12]:
a, b = 0, 1
Nmax = 10

for N in range(2, Nmax+1):
    x = points_NC_fermee(a, b, N)
    print("-"*80)
    print(f"poids pour la méthode fermée à {N} points :")
    print(*NC_poids(a, b, x), sep='\n')
    
for N in range(1, Nmax+1):
    x = points_NC_ouverte(a, b, N)
    print("-"*80)
    print(f"poids pour la méthode ouverte à {N} points :")
    print(*NC_poids(a, b, x), sep='\n')

--------------------------------------------------------------------------------
poids pour la méthode fermée à 2 points :
0.5
0.5
--------------------------------------------------------------------------------
poids pour la méthode fermée à 3 points :
0.16666666666666666
0.6666666666666666
0.16666666666666669
--------------------------------------------------------------------------------
poids pour la méthode fermée à 4 points :
0.12499999999999999
0.3749999999999999
0.375
0.12500000000000003
--------------------------------------------------------------------------------
poids pour la méthode fermée à 5 points :
0.07777777777777778
0.35555555555555557
0.13333333333333333
0.35555555555555557
0.07777777777777779
--------------------------------------------------------------------------------
poids pour la méthode fermée à 6 points :
0.06597222222222222
0.2604166666666667
0.1736111111111111
0.1736111111111112
0.2604166666666665
0.06597222222222224
--------------------------------------------------------------------------------
poids pour la méthode fermée à 7 points :
0.04880952380952381
0.25714285714285706
0.03214285714285711
0.32380952380952394
0.032142857142856876
0.2571428571428573
0.04880952380952384
--------------------------------------------------------------------------------
poids pour la méthode fermée à 8 points :
0.04346064814814816
0.20700231481481474
0.07656250000000009
0.17297453703703686
0.17297453703703722
0.07656249999999981
0.20700231481481493
0.043460648148148186
--------------------------------------------------------------------------------
poids pour la méthode fermée à 9 points :
0.034885361552028225
0.20768959435626105
-0.032733686067019464
0.3702292768959436
-0.16014109347442693
0.37022927689594365
-0.032733686067019506
0.20768959435626103
0.03488536155202823
--------------------------------------------------------------------------------
poids pour la méthode fermée à 10 points :
0.031886160714285725
0.1756808035714286
0.012053571428571518
0.2158928571428567
0.0644866071428578
0.06448660714285649
0.2158928571428577
0.012053571428571053
0.17568080357142876
0.03188616071428573
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 1 points :
1.0
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 2 points :
0.4999999999999999
0.5000000000000001
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 3 points :
0.6666666666666666
-0.3333333333333334
0.6666666666666666
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 4 points :
0.45833333333333337
0.041666666666666435
0.04166666666666691
0.45833333333333315
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 5 points :
0.5499999999999999
-0.7000000000000003
1.3000000000000012
-0.7000000000000013
0.5500000000000005
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 6 points :
0.42430555555555544
-0.31458333333333294
0.3902777777777766
0.3902777777777792
-0.3145833333333345
0.424305555555556
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 7 points :
0.48677248677248686
-1.00952380952381
2.323809523809525
-2.602116402116403
2.323809523809525
-1.00952380952381
0.4867724867724869
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 8 points :
0.39888392857142857
-0.6256696428571423
1.1087053571428553
-0.381919642857139
-0.3819196428571475
1.1087053571428616
-0.6256696428571451
0.3988839285714291
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 9 points :
0.4458774250440916
-1.28858024691358
3.6750440917107587
-6.070326278659613
7.475970017636678
-6.070326278659606
3.6750440917107574
-1.2885802469135796
0.4458774250440918
--------------------------------------------------------------------------------
poids pour la méthode ouverte à 10 points :
0.37925443672839537
-0.9098323137125229
2.159650573192242
-2.3541688712522113
1.2250961750441025
1.225096175044077
-2.354168871252192
2.159650573192231
-0.9098323137125186
0.3792544367283944
In [ ]:
In [ ]: