Python 21-22 for dummies : problème 5

Auteurs : Loïc Beuken, Olivier Cheffert

De nouveaux polynômes magiques ukrainiens !

Dans ce devoir, on doit réaliser une approximation polynomiale intégrale d’une fonction $u(x)$ aux sens des moindres carrés. Le but sera de calculer les coefficients $a_i$ approximation polynomiale de degré $p$

$$ u^h(x) = \sum_{i=0}^{p} a_i Z_i(x)$$

où les fonctions $Z_i(x)$ sont les célèbres polynômes ukrainiens de Zelensky.

Comme vous pouvez vous en douter, ces polynômes portent un autre nom. Il s'agit en réalité des polynômes de Legendre! Et oui, ces polynômes servent aussi à déterminer les points d'intégrations utilisé par la méthode de Gauss-Legendre. Dans la suite de Notebook, on remplecara les $Z$ par $\mathcal{P}$.

Les polynômes de Legendre sont définis sur l'intervalle $[-1, 1]$ et peuvent être obtenus par une formule de récurrence :

$$ \mathcal{P_0}(x) = 1 $$$$ \mathcal{P_1}(x) = x $$$$ ... $$$$ n \mathcal{P_n}(x) = (2n-1)x\mathcal{P_{n-1}}(x) - (n-1)\mathcal{P_{n-2}}(x) $$

Graphiques des polynômes de Legendre

import numpy as np

def Legendre(x,p):
    P = np.ones((p+1,len(x)))
    P[1] = x
    for n in range(2,p+1):
        P[n] = (2 - 1/n) *x*P[n-1] - (1-1/n)*P[n-2]
    return P

p = 5
m = 100
x = np.linspace(-1,1,m)

from matplotlib import pyplot as plt
plt.rcParams['toolbar'] = 'None'
plt.rcParams['figure.facecolor'] = 'moccasin'
  
fig = plt.figure("Polynomes de Legendre : p = %d :-)" % p)
plt.plot(x,Legendre(x,p).T)  
plt.show()

Propriétés des polynômes de Legendre

Dans ce devoir, on va tirer profit du fait que les polynôme $\mathcal{P_i}(x)$ sont des polynômes orthogonaux pour le produit scalaire de l'espace $L^2(\Omega)$ avec $\Omega = ] − 1, 1[$. En d’autres mots, les polynômes de Legendre sont construits afin de satisfaire la propriété suivante :

$$\int_{-1}^{1} \mathcal{P_i}(x) \mathcal{P_j} (x) dx = 0 \text{ si } i \neq j \;\; (1)$$

Ceci va nous permettre de ne pas devoir résoudre de système pour déterminer les coefficients $a_i$. Mais que se passe-t-il quand $i = j$ ? En cherchant un peu sur internet, on peut trouver que les polynômes de Legendre respectent aussi la propriété suivante :

$$ \int_{-1}^{1} \mathcal{P_n}(x) \mathcal{P_n} (x) dx = \frac{2}{2n + 1} \; \; (2)$$

ce résultat nous sera utile pour la suite du devoir!

Déterminer les coefficients $a_i$

Les coefficients $a_i$ sont obtenus afin de minimiser l'intégrale du carré de l'écart $u(x) − u^h(x)$ sur l'intervalle $[−1, 1]$. Il s’agit donc de calculer les paramètres $a_j$ de l’approximation par moindres carrés en minimisant l’intégrale du carré de l'écart $u(x)−u^h(x)$ en tout point de l’intervalle

$$ J(a_0, a_1, \cdots, a_n) = \int_{-1}^{1} \left(u(x) - \sum_{j=0}^{p} a_j \mathcal{P_j}(x) \right)^2 dx$$

Afin de minimiser la fonction objective, on doit calculer le gradient de $J$. Avant de faire ce calculs, simplifions notre fonction objective en dévelopant le carré. Pour cela prenont un exemple un peu plus simple pour commencer : $(z - (a+b+c+d))^2$ où $a$, $b$, $c$ et $d$ correspond aux éléments de la somme dans $J$ et $z = u(x)$. En dévelopant, on obtient:

$$a^2 + 2 a b + 2 a c + 2 a d - 2 a z + b^2 + 2 b c + 2 b d - 2 b z + c^2 + 2 c d - 2 c z + d^2 - 2 d z + z^2$$

En utilisant la propriété $(1)$ des polynômes de Legendre, on constate que la plus part des termes seront nuls. On obtiens donc:

$$ J = \int_{-1}^{1} \sum_{j=0}^p a_j^2 \mathcal{P_j}^2(x) - 2 \sum_{j=0}^{p} a_j \mathcal{P_j}(x) u(x) dx$$$$ J = \sum_{j=0}^p a_j^2 \int_{-1}^{1} \mathcal{P_j}^2(x) dx - 2 \sum_{j=0}^{p} a_j \int_{-1}^{1} \mathcal{P_j}(x) u(x) dx $$

On calcule ensuite le gradient : $$ \frac{\partial J}{\partial a_j} = 2 a_j \int_{-1}^{1} \mathcal{P_j}^2(x) dx - 2 \int_{-1}^{1} \mathcal{P_j}(x) u(x) dx$$ Dont la condition a imposé est $\frac{\partial J}{\partial a_j} = 0$. En utilisant la propriété $(2)$ des polynômes de Legendre, on trouve donc le résultat final :

$$ a_j = \frac{2j+1}{2} \int_{-1}^{1} \mathcal{P_j}(x) u(x) dx $$

On peut déterminer ce résultat en regardant directement la matrice des équations normales à résoudre pour l'approximation polynomiale dont on connait la fonction $u(x)$ (page 27 du syllabus).

Pour obtenir une estimation de l'intgrale se trouvant dans l'expression des coefficients $a_j$, on utilisera uniquement les points et les poids de Gauss-Legendre qui seront fournis en argument $Xint$, $Wint$.

Le code

def legendreApproximation(u, p, x, Xint, Wint):
        
    P = Legendre(Xint,p)
    a = np.zeros(p+1)
    Uint = u(Xint)
    for j in range(0,p+1):
        a[j] = (2*j+1)/2 * (P[j]*Uint @ Wint)
    
    P = Legendre(x,p)
    uh = a @ P
    return uh,a

Le test

u = lambda x : np.sin(1.5*x)+0.5*np.sin(5*1.5*x)

from scipy.special import roots_legendre
Xint,Wint = roots_legendre(p+1)
uh,a  = legendreApproximation(u,p,x,Xint,Wint)

fig = plt.figure("Approximation de Zelensky : p = %d :-)" % p)
plt.plot(x,u(x),'-b',label='Fonction à approximer !')  
plt.plot(x,uh,'-r',label='Approximation :-)')
plt.legend(loc='upper right')
plt.show()