DEVOIR ¶

Choix des points d'interpolation¶

Dans ce devoir, 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'
color_marron = 'xkcd:burnt sienna'
color_violet = 'xkcd:blue violet'

Rappel : une fonction qui évalue le polynôme interpolateur de Lagrange¶

In [2]:
def interp_vdm(x, y, xx):
    """
    calcule et évalue le polynôme interpolateur de Lagrange
    en utilisant la méthode de la matrice de Vandermonde
    
    Parameters
    ----------
    
    x: ndarray
        abscisses des points d'interpolation
        
    y: ndarray
        ordonnées des points d'interpolation
        
    xx: ndarray
        abscisses des points d'évaluation
        
    Returns
    -------
    
    yy: ndarray
        ordonnées des points d'évaluation
    """
    if x.size != y.size:
        msg = "Erreur dans interp_vdm : x et y doivent avoir la même taille\n"
        msg += f"taille de x : {x.size}\n"
        msg += f"taille de y : {y.size}\n"
        raise ValueError(msg)
    n = x.size
    A = np.fromfunction(
        lambda i, j: x[i]**j, (n, n), dtype='int'
    )
    P = np.linalg.solve(A, y)
    yy = 0.
    for ak in P[::-1]:
        yy *= xx
        yy += ak
    return yy
In [3]:
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):
        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])
        yy += y[i] * li
    return yy

Exercice 1 : interpolation et approximation¶

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$.

Indication : vous pourrez utiliser les fonctions vectorielles de numpy(par exemple np.linspace, np.cos) et la valeur de $\pi$ np.pi.

In [4]:
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

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 ?

Indication : vous pouvez utiliser une des fonctions d'interpolation codée au TP01.

In [5]:
def f(x):
    return 1./(1+x**2)

a, b = -5, 5
In [6]:
interp = interp_vdm
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(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)

EXERCICE 2 : une autre manière d'interpoler¶

Dans cet exercice, nous étudions le comportement du polynôme interpolateur de Lagrange lorsque les points d'interpolation se rapprochent. Pour simplifier, nous supposerons ici que le nombre de points d'interpolation est égal à $4$.

Nous rappelons que, étant donnée une fonction $f:[a, b]\to\mathbb{R}$, le polynôme interpolateur aux points $a\leq x_0 < x_1 < x_2 < x_3 \leq b$ est défini comme l'unique polynôme $P$ de $\mathbb{R}_3[X]$ tel que $P(x_i)=f(x_i)$, $0\leq i<4$.

Dans cette définition, nous avons obligatoirement des points $x_i$ deux à deux disjoints.

Nous posons à présent $$ x_0=0, \qquad x_1 = h, \qquad x_2=1-h, \qquad x_3=1, $$ et nous considérons $$ f(x) = 1 - e^{-10x}. $$

Question

Dans une fenêtre graphique,

  • tracez la fonction $f$ sur l'intervalle $[0,1]$,
  • les polynômes interpolateurs de Lagrange obtenus avec les 4 points pour $h=0.3$,
  • ajoutez sur la même figure les courbes obtenues pour $h=0.1$, $h=0.01$ et $h=0.001$.

Que constatez-vous ?

Vous pourrez essayer de produire une image ressemblant à celle-ci : DM_Hermite

In [7]:
def f(x):
    return 1 - np.exp(-10*x)

def points_x(h):
    return np.array([0, h, 1-h, 1])
In [9]:
taille = 4

xx = np.linspace(0, 1, 1025)

fig = plt.figure(figsize=(taille, taille))
ax = fig.add_subplot(1, 1, 1)
ax.plot(
    xx, f(xx),
    linewidth=2, alpha=0.5, color=color_marron, linestyle='dashed',
    label='fonction'
)

for h in [0.3, 0.1, 0.01, 0.001]:
    x = points_x(h)
    y = f(x)
    yy = interp_Lagrange(x, y, xx)

    ax.scatter(x, y, s=50, color=color_bleu, alpha=1)
    ax.plot(
        xx, yy,
        linewidth=2, alpha=0.75,
        label=f'$h={h}$'
    )
ax.set_title(f"Interpolation à 4 points", fontsize=16, color=color_bleu)
ax.grid(True, color=color_bleu, linestyle='dotted', alpha=0.5)
ax.legend()
fig.savefig("DM_Hermite.png")
In [ ]: