FEUILLE DE TP 1

Construction du polynôme interpolateur de Lagrange¶

Dans ce TP, nous construisons le polynôme interpolateur de Lagrange par deux méthodes (écritures dans la base canonique et dans la base des polynômes 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'

Nous rappelons le résultat d'existence et d'unicité vu en cours

(Polynôme interpolateur de Lagrange)
Etant donnés $N$ un entier strictement positif, $x_1,\ldots,x_N$ des réels deux à deux distincts et $y_1,\ldots,y_N$ des réels, il existe un unique polynôme $P$ qui vérifie $$ P\in\mathbb{R}_{N-1}[X], \qquad P(x_i)=y_i, \quad 1\leq i\leq N.$$

Construction par la méthode de la matrice de Vandermonde¶

Nous rappelons que le vecteur $a=(a_1,\ldots,a_N)$ des coefficients du polynôme interpolateur de Lagrange dans la base canonique de $\mathbb{R}_{N-1}[X]$ est solution du système linéaire $$ M a = y$$ où $M$ est la matrice de Vandermonde associée aux points $(x_1,\ldots, y_N)$ et $y=(y_1,\ldots,y_N)$. Nous avons $$ M = \begin{pmatrix} 1&x_1&\ldots&x_1^{N-1}\\ 1&x_2&\ldots&x_2^{N-1}\\ \vdots&\vdots&&\vdots\\ 1&x_N&\ldots&x_N^{N-1} \end{pmatrix}. $$

Question

  1. Proposez une fonction interp_vdm_build qui prend en argument un ndarray noté x et qui retourne la matrice de Vandermonde associée.
  2. Proposez une fonction interp_vdm_poly qui
    • prend en argument deux ndarray notés x et y ;
    • vérifie que les deux arguments ont bien la même taille ;
    • calcule la matrice de Vandermonde M associée au vecteur x ;
    • retourne la solution du système Ma=y.
  3. Vérifiez que vos fonctions fonctionnent correctement en affichant les résultats pour des points bien choisis (vous pouvez prendre ce que vous voulez...)

Indication : vous pourrez utiliser la fonction fromfunction du module numpy, la fonction solve du module numpy.linalg.

In [2]:
def interp_vdm_build(x):
    """
    retourne la matrice de Vandermonde associée à x
    
    Parameters
    ----------
    
    x: ndarray
        les valeurs des points
        
    Returns
    -------
    
    out: ndarray
        une matrice carrée de taille NxN (N est la taille de x)
        out[i, j] = x[i]**j
    
    """
    n = x.size
    return np.fromfunction(
        lambda i, j: x[i]**j, (n, n), dtype='int'
    )


def interp_vdm_poly(x, y):
    """
    retourne le polynôme interpolateur de Lagrange
    
    Parameters
    ----------
    
    x: ndarray
        les abscisses des points d'interpolation
    y: ndarray
        les ordonnées des points d'interpolation
        
    Returns
    -------
    
    P: ndarray
        les coefficients dans la base canonique du polynôme interpolateur
        P = P_0 + P_1 X + P_2 X^2 + ... + P_n X^n
        
    Remark
    ------
    
    Les deux ndarray x et y doivent avoir la même taille
    """
    if x.size != y.size:
        msg = "Erreur dans interp_vdm_poly : 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)
    return np.linalg.solve(interp_vdm_build(x), y)
In [3]:
# test avec le polynôme X^2-1
x = np.array([0., 1, 2])
y = x**2-1
print(interp_vdm_build(x))
print(interp_vdm_poly(x, y))  # [-1, 0, 1]

[[1. 0. 0.]
 [1. 1. 1.]
 [1. 2. 4.]]
[-1.  0.  1.]

Question

  1. Proposez une fonction horner qui prend en argument un ndarray noté P et un ndarray noté xx et qui retourne l'évaluation aux points du vecteur xx par l'algorithme de Hörner du polynôme dont les coefficients dans la base canonique sont stockés dans le vecteur P.
  2. Proposez une fonction interp_vdm qui
    • prend en argument deux ndarray x et y de taille N et un ndarray xx de taille M ;
    • calcule le polynôme interpolateur de Lagrange aux points donnés par les vecteurs x et y en utilisant la fonction interp_vdm_poly ;
    • retourne l'évaluation de ce polynôme aux points du vecteur xx en utilisant la fonction horner.
In [4]:
def horner(P, xx):
    """
    évalue par l'algorithme de Hörner le polynôme
    dont les coefficients dans la base canonique sont stockés dans le vecteur P
    aux points stockés dans le vecteur xx
    
    Parameters
    ----------
    
    P: ndarray
        Les coefficients dans la base canonique du polynôme P
    xx: ndarray
        Les abscisses des points où l'on évalue le polynôme P
    """
    yy = 0.
    for ak in P[::-1]:
        yy *= xx
        yy += ak
    return yy


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)
    a = interp_vdm_poly(x, y)
    return horner(a, xx)

Question

Afin de tester votre fonction : en prenant $N=5$,

  1. prenez $N$ points équirépartis entre 0 et 1 (ce sera notre vecteur x) ;
  2. générez $N$ valeurs aléatoires entre 0 et 1 (ce sera notre vecteur y) ;
  3. tracez dans une fenêtre graphique le nuage de points d'abscisses x et d'ordonnées y à l'aide d'une commande scatter ;
  4. ajoutez le tracé du polynôme interpolateur en prenant xx un vecteur de taille grande (plutôt 100 ou 200 points équi-répartis entre 0 et 1) ;
  5. vérifiez que le polynôme interpolateur passe bien par les points d'interpolation.
  6. recommencez avec des valeurs de $N$ plus grandes. Que constatez-vous ?
In [21]:
for N in range(5, 30, 5):
    x = np.linspace(0, 1, N)
    y = np.random.rand(N)
    xx = np.linspace(0, 1, 1025)
    yy = interp_vdm(x, y, xx)

    fig = plt.figure(figsize=(12, 4))
    ax = fig.add_subplot(1, 2, 1)
    ax.scatter(x, y, s=50, color=color_bleu, alpha=1)
    ax.plot(xx, yy, linewidth=2, alpha=0.75, color=color_orange)
    ax.set_title(f"Polynôme interpolateur à {N} points", fontsize=16, color=color_bleu)
    ax.grid(True, color=color_bleu, linestyle='dotted', alpha=0.5)
    ax = fig.add_subplot(1, 2, 2)
    ax.scatter(x, y, s=50, color=color_bleu, alpha=1)
    ax.plot(xx, yy, linewidth=2, alpha=0.75, color=color_orange)
    ax.set_ylim(-1.5, 1.5)
    ax.set_title("zoom", fontsize=16, color=color_bleu)
    ax.grid(True, color=color_bleu, linestyle='dotted', alpha=0.5)

Construction en utilisant la base duale des polynômes de Lagrange¶

Nous rappelons que les polynômes de Lagrange associés aux points $x_1,\ldots,x_N$ sont définis par $$ L_i = \prod_{j\neq i} \frac{X-x_j}{x_i-x_j}, \qquad 1\leq i\leq N. $$ Le polynôme interpolateur de Lagrange s'exprime alors dans la base des polynômes de Lagrange $$ P = \sum_{i=1}^N y_i L_i. $$

Question

Proposez une fonction interp_Lagrange qui

  • prend en argument deux ndarray x et y de taille N et un ndarray xx de taille M ;
  • retourne l'évaluation du polynôme interpolateur aux points du vecteur xx en utilisant la décomposition dans la base des polynômes de Lagrange. Testez votre nouvelle fonction en reprenant le test fait avec la méthode de Vandermonde. Voyez-vous une différence en particulier pour les grandes valeurs de $N$ ?
In [22]:
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
    -------
    
    yy: ndarray
        ordonnées des points d'évaluation
    """
    if x.size != y.size:
        msg = "Erreur dans interp_Lagrange : 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
    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
In [23]:
for N in range(5, 30, 5):
    x = np.linspace(0, 1, N)
    y = np.random.rand(N)
    xx = np.linspace(0, 1, 1025)
    yy = interp_Lagrange(x, y, xx)

    fig = plt.figure(figsize=(12, 4))
    ax = fig.add_subplot(1, 2, 1)
    ax.scatter(x, y, s=50, color=color_bleu, alpha=1)
    ax.plot(xx, yy, linewidth=2, alpha=0.75, color=color_orange)
    ax.set_title(f"Polynôme interpolateur à {N} points", fontsize=16, color=color_bleu)
    ax.grid(True, color=color_bleu, linestyle='dotted', alpha=0.5)
    ax = fig.add_subplot(1, 2, 2)
    ax.scatter(x, y, s=50, color=color_bleu, alpha=1)
    ax.plot(xx, yy, linewidth=2, alpha=0.75, color=color_orange)
    ax.set_ylim(-1.5, 1.5)
    ax.set_title("zoom", fontsize=16, color=color_bleu)
    ax.grid(True, color=color_bleu, linestyle='dotted', alpha=0.5)
In [ ]: