FEUILLE DE TP 2

Choix des points d'interpolation¶

Dans ce TP, nous étudions l'importance du choix des points d'interpolation en particulier lorsque l'interpolation est utilisée pour approcher une fonction.

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'

Etant donné un intervalle $[a, b]$, nous rappelons la définition des familles de points classiques pour l'interpolation :

  • les points équi-répartis $$ x_{i, N} = a + (b-a) \frac{i-1}{N-1}, \qquad 1\leq i\leq N,$$
  • les points de Tchebychev $$ x_{i, N} = \frac{a+b}{2} + \frac{b-a}{2} \cos \Bigl( \frac{2i-1}{2N}\pi \Bigr), \qquad 1\leq i\leq N,$$

Question 1

  1. Proposez une fonction equirepartis(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 la famille des $N$ points équi-répartis entre $a$ et $b$.
  2. Proposez une fonction tchebychev(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 la famille des $N$ points de Tchebychev entre $a$ et $b$.
In [2]:
def equirepartis(a, b, N):
    """
    retourne les N points équirépartis entre a et b
    
    Parameters
    ----------
    
    a: float
        le bord gauche de l'intervalle
    b: float
        le bord droit de l'intervalle
    N: int
        le nombre de points
        
    Returns
    -------
    
    out: ndarray
        les N valeurs des points équi-répartis
    """
    i = np.arange(1, N+1)
    return a + (b-a)/(N-1) * (i-1)


def tchebychev(a, b, N):
    """
    retourne les N points de Tchebychev entre a et b

    Parameters
    ----------
    
    a: float
        le bord gauche de l'intervalle
    b: float
        le bord droit de l'intervalle
    N: int
        le nombre de points
        
    Returns
    -------
    
    out: ndarray
        les N valeurs des points de Tchebychev
    """
    i = np.arange(N, 0, -1)
    return .5*(a+b) + .5*(b-a)*np.cos((2*i-1)/(2*N)*np.pi)

Question 2

  1. Proposez une fonction Li(i, x, xx) qui prend en arguments un entier $i$ et deux vecteurs $x$ et $xx$ et qui retourne (si c'est possible) l'évaluation du polynôme $L_i$ aux points $xx$ avec pour rappel $$ L_i(X) = \prod_{j\neq i} \frac{X-x_j}{x_i-x_j}.$$
  2. Tracez dans des fenêtres graphiques séparés l'allure des différents polynômes pour les points équi-répartis et les points de Tchebychev (pour différentes valeurs de $N$).
  3. Que remarquez-vous ?
In [3]:
def Li(i, x, xx):
    """
    retourne l'évaluation en xx du polynôme Li construit avec les points x
    
    Parameters
    ----------
    
    i: int
        le ieme polynôme de la base duale de Lagrange
    x: ndarray
        les points d'interpolation
    xx: ndarray
        les points d'évaluation
        
    Returns
    -------
    
    yy: ndarray
        les évaluations yy_j = L_i(xx_j)
    
    Remark
    ------
    
    i doit être un entier positif ou nul et strictement inférieur à la taille de x
    """
    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)
    li = 1.
    for j in range(i):
        li *= (xx-x[j]) / (x[i] - x[j])
    for j in range(i+1, n):
        li *= (xx-x[j]) / (x[i] - x[j])
    return li
In [4]:
a, b = 0, 1
vN = [3, 5, 11, 21]
xx = np.linspace(a, b, 1025)

fig = plt.figure(figsize=(8, 4*len(vN)))

k = 0
for N in vN:
    for title, points in zip(
        ["Points équi-répartis", "Points de Tchebychev"],
        [equirepartis, tchebychev]
    ):
        k += 1
        ax = fig.add_subplot(len(vN), 2, k)
        ax.set_title(title + f" ($N={N}$)")
        x = points(a, b, N)
        ax.scatter(x, 0*x, s=50, color='black')
        ax.scatter(x, 0*x+1, s=50, color='black')
        for i in range(N):
            ax.plot(xx, Li(i, x, xx), color=(i/N, 0, 1-i/N), alpha=0.5)

Question 3

Afin de visualiser la convergence ou la non convergence du polynôme interpolateur vers une fonction régulière, nous choisissons de prendre la fonction $f$ définie sur $[-5, 5]$ par $$f(x) = \frac{1}{1+x^2}.$$

  1. Tracez dans une même fenêtre graphique
    • la fonction $f$
    • les polynômes interpolateurs de Lagrange pour $N$ points équi-répartis avec $N\in\lbrace 5, 10, 20 \rbrace$. Ajoutez un titre et une légende.
  2. Recommencez la question précédentes avec des points de Tchebychev.
  3. Que remarquez-vous ?
In [5]:
def interp_Lagrange(x, y, xx):
    """
    calcule et évalue le polynôme interpolateur de Lagrange
    en utilisant la base duale des polynôme de Lagrange
    
    Parameters
    ----------
    
    x: ndarray
        abscisses des points d'interpolation
        
    y: ndarray
        ordonnées des points d'interpolation
        
    xx: ndarray
        abscisses des points d'évaluation
        
    Returns
    -------
    
    ndarray
        ordonnées des points d'évaluation
    """
    n = x.size
    if n != y.size:
        msg = "Problème de taille dans interp_Lagrange\n"
        msg += f"x.size = {x.size}, y.size={y.size}"
        raise ValueError(msg)
    yy = 0.
    for i in range(n):
        yy += y[i] * Li(i, x, xx)
    return yy
In [6]:
def f(x):
    return 1./(1+x**2)

a, b = -5, 5
In [7]:
vect_N = [5, 10, 20]
vect_color = [color_bleu, color_vert, color_rouge]
xx = np.linspace(a, b, 1025)

fig = plt.figure(figsize=(12, 6))
ax_equi = fig.add_subplot(1, 2, 1)
ax_tche = fig.add_subplot(1, 2, 2)

for N, color in zip(vect_N, vect_color):
    xequi = equirepartis(a, b, N)
    xtche = tchebychev(a, b, N)
    for x, ax in zip([xequi, xtche], [ax_equi, ax_tche]):
        yy = interp_Lagrange(x, f(x), xx)
        ax.plot(xx, yy, label=f'$N={N}$', color=color)

for title, ax in zip(
    ['Points équi-répartis', 'Points de Tchebychev'],
    [ax_equi, ax_tche]
):
    ax.plot(xx, f(xx), color='black', label=r'$f$', lw=2, linestyle='dashed')
    ax.legend()
    ax.set_xlim(a-.5, b+.5)
    ax.set_ylim(-1, 2)
    ax.set_title(title, fontsize=20, color=color_bleu)

Question 4

Pour mesurer l'écart entre la fonction $f$ et un polynôme interpolateur $P$, nous choisissons la norme infinie.

  1. Tracez dans une fenêtre graphique la norme infinie de la différence entre la fonction et son polynôme interpolateur à $N$ points équi-répartis en fonction de $N$. Vous prendrez $N$ entre $5$ et $100$ par pas de $5$ et vous choisirez une échelle logarithmique pour les ordonnées.
  2. Ajoutez sur la même figure la courbe obtenue pour les points de Techbychev.
  3. Que remarquez-vous ?
In [8]:
vect_N = np.arange(5, 100, 5)
xx = np.linspace(a, b, 1025)
yy = f(xx)

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1)

err_equi = np.zeros(vect_N.shape)
err_tche = np.zeros(vect_N.shape)

for k, N in enumerate(vect_N):
    xequi = equirepartis(a, b, N)
    xtche = tchebychev(a, b, N)
    yyequi = interp_Lagrange(xequi, f(xequi), xx)
    yytche = interp_Lagrange(xtche, f(xtche), xx)
    err_equi[k] = np.linalg.norm(yy - yyequi, np.inf)
    err_tche[k] = np.linalg.norm(yy - yytche, np.inf)

ax.scatter(vect_N, err_equi, color=color_orange, marker='s', label='équi-répartis')
ax.scatter(vect_N, err_tche, color=color_vert, marker='o', label='Tchebychev')

ax.set_xlabel(r'$N$')
ax.set_ylabel(r'$\Vert f-P_N\Vert_\infty$')
ax.set_xscale('linear')
ax.set_yscale('log')
ax.set_title("Evolution de l'erreur", fontsize=16, color=color_bleu)
ax.legend();
In [ ]: