De nouvelles méthodes numériques

De nouvelles méthodes numériques#

Dérivée numérique décentrée#

Pour un ensemble de points de coordonnées \(\{(x_i,y_i)\}\) où \(y\) est une fonction de \(x\), nous avons vu qu’une première approximation pour calculer la dérivée \(y' = \frac{\mathrm{d}y}{\mathrm{d}x}\) est de calculer la pente entre deux points consécutifs.

\[ y'_i = \frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}} \]

Il est cependant arbitraire d’associer la dérivée en un point avec la pente par rapport au point suivant, plutôt que la pente par rapport au point précédent. Une deuxième approximation peut donc s’écrire :

\[ y'_i = \frac{y_{i}-y_{i-1}}{x_{i}-x_{i-1}} \]

Ces deux approximations sont appelées respectivement dérivée décentrée à droite et dérivée décentrée à gauche.

Comme on peut le voir sur la figure ci-dessus, ces deux approximations (en vert et en bleu) ne sont pas toujours satisfaisantes car elles sont arbitrairement biaisées d’un côté ou de l’autre du point considéré. Une troisième approximation consiste à prendre la pente entre les deux points adjacents au point considéré. Cette approximation est appelée dérivée centrale :

\[ y'_i = \frac{y_{i+1}-y_{i-1}}{x_{i+1}-x_{i-1}} \]

Programmez une fonction pour ces trois différentes approximations. La fonction pourra prendre en arguments des points \(\{(x_i,y_i)\}\) sous forme de listes ou de tableaux numpy (plusieurs solutions sont possibles).

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

# Dans toutes ces fonctions on sera particulierement attentif aux indices, surtout en début et fin de tableau

def derivdroite(x, y): # On attend en entrée 2 listes ou 2 tableaux numpy de taille n
    n = len(x)
    # On commence par initialiser les tableaux de sorties
    # Attention : pour la dérivée à droite on ne peut pas calculer le dernier point !
    # Les tableaux de sorties seront donc de taille n - 1
    xoutput = np.zeros(n - 1)
    youtput = np.zeros(n - 1)
    for i in range(n - 1):
        xoutput[i] = x[i]
        youtput[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
    return [xoutput, youtput] # La sortie est une liste de deux tableaux numpy de tailles n - 1

def derivgauche(x, y): # On attend en entrée 2 listes ou 2 tableaux numpy de taille n
    n = len(x)
    # On commence par initialiser les tableaux de sorties
    # Attention : pour la dérivée à gauche on ne peut pas calculer le premier point !
    # Les tableaux de sorties seront donc de taille n - 1
    xoutput = np.zeros(n - 1)
    youtput = np.zeros(n - 1)
    for i in range(n - 1):
        # Ici le premier point du tableau de sortie est le deuxieme point du tableau d'entrée
        xoutput[i] = x[i + 1] 
        youtput[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
    return [xoutput, youtput] # La sortie est une liste de deux tableaux numpy

def derivcentre(x, y): # On attend en entrée 2 listes ou 2 tableaux numpy de taille n - 1
    n = len(x)
    # On commence par initialiser les tableaux de sorties
    # Normalement on ne peut pas calculer le premier ni le dernier point !
    # Mais pour les bords de l'intervalle, on utilise une dérivée gauche ou droite.
    # Ce choix n'est pas sans conséquences.
    # Les tableaux de sortie auront donc la même taille que les tableaux d'entrée...
    xoutput = np.zeros(n)
    youtput = np.zeros(n)
    # Premier point
    xoutput[0] = x[0]
    youtput[0] = (y[1] - y[0]) / (x[1] - x[0])
    # Dernier point
    xoutput[n - 1] = x[n - 1]
    youtput[n - 1] = (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])
    # Les points au milieu...
    for i in range(n - 2):
        xoutput[i + 1] = x[i + 1]
        youtput[i + 1] = (y[i + 2] - y[i]) / (x[i + 2] - x[i])
    return [xoutput, youtput] # La sortie est une liste de deux tableaux numpy de tailles n

Comparez les résultats obtenus pour les trois fonctions mathématiques connues que vous avez dérivées mathématiquement dans le calepin à préparer pour aujourd’hui, et comparez au résultat théorique. On prendra pour la liste de points en abscisse un intervalle contenu dans l’intervalle de définition de la fonction et un pas de 0,1 pour commencer.

../../../_images/fb67fdfc0eb693ad46fe3b5ca8c49cd6d0d7642ad3a8b7bbcce6a0536bafe5c1.png

../../../_images/46e33ec252faceedf93a1fdc9f363bc617022e8ab15bd24757fa475de3d3621c.png

../../../_images/77ecbfd442f369ad7e3ae17f9195b4996918e2bf3edfdfd9200f73b9dd4a4394.png

Quantifiez et affichez l’erreur entre les 3 méthodes numériques, et le résultat théorique. Faites varier le pas de la dérivée et observez son influence sur la précision obtenue.

../../../_images/26ffa97381cd74820d1a79fe858ad24a5351fb190a180081f7a325080871d5ae.png

../../../_images/fb8e925254ebca769b29f26da46cdb320c61b203713177f25c0e00a2c1679003.png

../../../_images/4a2afd14c12b51d2f075622e1f92b509cbc4cd1adf14d71875bec0513a712330.png

Dérivée d’une courbe expérimentale#

En physique expérimentale, il est rare d’analyser des points bien sagement ordonnés avec une précision à 12 chiffres. De plus, toute mesure expérimentale contient du bruit.

Dans la suite, vous allez essayer de dériver un exemple de courbe expérimentale à partir des données brutes.

Importez le fichier de données dataBSCCO.csv (par exemple avec le module pandas) qui contient une mesure de résistance d’un échantillon 2D de supraconducteur à haute température critique (le \(\mathrm{Bi}_2\mathrm{Sr}_2\mathrm{CaCu}_2\mathrm{O}_{8+\delta}\)) en fonction de la température. Affichez quelques lignes de son contenus, ainsi qu’un graphe des données Rs (résistance 2D en ohms par carré - \(\Omega/\square\) abbrégé Ohmsq) en fonction de T (température en kelvins).

         T(K)   Rs(Ohmsq)
    29.527     1.85548
    29.565     7.28151
    29.693     3.79314
    29.881   -10.66002
    30.110     2.55768
...       ...         ...
193.020  5059.90546
193.000  5057.80715
193.020  5058.66410
193.010  5059.54701
193.020  5058.83260

[8583 rows x 2 columns]

../../../_images/3b7618560f7228c58f66e0fea9252b58eb6137c3bb35894595e4cd06128646c4.png

Vous pouvez essayer d’utiliser les méthodes de dérivation numériques précédentes sur cette courbe. A priori ça ne marchera pas car les données T(K) ne sont pas uniformément réparties et ne varient pas de façon strictement monotone.

Avant de dériver cette courbe, il faut donc mettre en forme les données. On peut par exemple :

Interpoler les données sur un tableau régulier avec la fonction numpy.interp.
Moyenner localement la courbe. Créer un tableau régulier en abscisse {\(x_i\)} avec un pas \(p\), sélectionner les points de données tels que \(|T-x_i|<p\) et en faire une moyenne pour chaque \(i\).

Créez une nouvelle fonction qui permet de créer, à partir des données brutes, des données qui pourront être utilisées par les méthodes de dérivation créées précédemment.
Appliquez cette nouvelle fonction aux données brutes et étudiez le comportement de la dérivée en fonction du pas choisi.

Show code cell source Hide code cell source

# Avant de dériver, il faut mettre en forme les données
# Ici, on utilise une méthode de moyenne locale pour obtenir des points R(T) sur une grille régulière
# D'autres solutions sont possibles

def lavg(x, y, step): # Fonction moyenne locale 
    
    # L'idée est de créer un tableau régulier en abscisse avec un pas 'step'
    # puis de moyenner tous les points d'abscisses proches des points du nouveau tableau
    
    # Bornes des abcisses
    xmin = np.min(x) 
    xmax = np.max(x)
    
    # Premièrement on vérifie que le pas n'est pas trop petit pour les données
    # Si le pas est trop petit, on pourrait trouver un point sans données à moyenner
    # On cherche donc le pas le plus petit acceptable, 
    # c'est-a-dire l'écart maximum entre deux abscisses consécutives
    stepmin = 0 
    for i in range(len(x) - 1):
        if abs(x[i] - x[i + 1]) > stepmin: # La valeur absolue permet d'avoir un tableau x classé en ordre décroissant si on veut
            stepmin = abs(x[i] - x[i + 1])
    
    if step < stepmin: # Si le pas est trop petit pour la moyenne, on envoie un message d'erreur
        print("step too small")
        return 0
    
    # On crée le nouveau tableau d'abscisses
    nx = np.arange(xmin, xmax, step)
    # On initialise le nouveau tableau d'ordonnées
    ny = np.zeros(len(nx))
    
    for i in range(len(nx)): # Pour chaque point de nx
        xresult = [] # Initialisation de xresult
        for j in range(len(x)): # On cherche tous les points de x contenus dans l'intervalle autour du nx[i]
            if x[j] >= (nx[i] - step / 2) and x[j] < (nx[i] + step / 2):
                xresult.append(j) # On liste ces points dans xresult
        
        for j in range(len(xresult)): # On moyenne ensuite les ordonnées correspondantes
            ny[i] = ny[i] + y[xresult[j]]
        ny[i] = ny[i] / len(xresult)
        
    return [nx, ny] # On retourne une liste de 2 tableaux numpy [abscisses, ordonnées]



def derivclean(x, y, step): # Fonction combinant moyenne locale et dérivée
    cleaned = lavg(x, y, step)
    
    if cleaned == 0: # Si le pas est trop petit pour la moyenne, on envoie un message d'erreur
        print("step too small")
        return 0
    
    nx = cleaned[0]
    ny = cleaned[1]
    deriv = derivcentre(nx, ny) # On apelle la fonction de dérivation 
    
    return [nx, ny, deriv[1]] # On retourne à la fois la nouvelle abscisse, la courbe moyennée et sa dérivée

../../../_images/ec7137a12383ad260ced915ff17d77c73870c33ec8a3babba0bef9994423b7bf.png

../../../_images/c17234cbe463d81e9acba756a8f994c5dc779957a29a4d0fb150684f3ad8d53f.png

Intégration numérique#

Nous allons étudier plusieurs méthodes d’intégration numérique. Nous commençons par rappeler la méthode des rectangles vue dans la partie cours.

Méthode des rectangles#

Pour une fonction intégrable \(f\), on cherche à calculer l’intégrale suivante :

\[ I(f) = \int_{a}^{b} f(x)\mathrm{d}x \]

On choisit de discrétiser l’intervalle \([a,b]\) en \(n\) intervalles de même longueur, repérés par les points \(x_i = a+i\times \Delta x\) avec \(\Delta x=\frac{b-a}{n}\) . Une approximation de l’intégrale \(I(f)\) est donnée par :

\[ I(f) \approx I_n(f) = \sum_{i=0}^{n-1} f(x_i)\Delta x \]

Méthode du point milieu#

Une approximation plus précise de l’intégrale \(I(f)\) peut être calculée en considérant les valeurs de \(f\) au centre des \(n\) intervalles considérés. C’est-à-dire les points de coordonnées \(x_i = a + (i+1/2)\times\Delta x\). L’intervalle est alors calculé par l’approximation :

\[ J_n(f) = \sum_{i=0}^{n-1} f(x_i)\Delta x \]

Méthode des trapèzes#

Enfin, la méthode des trapèzes consiste à approximer la fonction \(f\) par une fonction affine par morceaux entre les points de coordonnées \(x_i = a + i\times\Delta x\). On peut alors calculer l’intégrale :

\[ K_n(f) = \sum_{i=0}^{n-1} \frac{f(x_i)+f(x_{i+1})}{2}\Delta x \]

Implémentez ces trois méthodes d’intégration numérique sous la forme de fonctions. On pourra mettre en arguments des fonctions : la fonction à intégrer, les bornes d’intégrations et le nombre de points d’intégration.

Appliquez ces méthodes d’intégration aux trois fonctions données en introduction et comparez aux valeurs analytiques que vous aurez calculées en choisissant des bornes appropriées.

exact :		 0.8999999999999999
rectangle :	 0.9017649037043426
milieu :	 0.9000175467586048
trapeze :	 0.8999649037043423

exact :		 0.6156264703860141
rectangle :	 0.6078596436599304
milieu :	 0.6156163647222771
trapeze :	 0.6156466822832047

exact :		 0.6201145069582775
rectangle :	 0.6189587691633259
milieu :	 0.6201147294096454
trapeze :	 0.620114062056476

Étudiez le comportement de l’erreur d’approximation avec le nombre de points. On pourra par exemple tracer un graphe de l’erreur en fonction du nombre de points (10, 100, 1000, etc.)

../../../_images/f9f69a78875c4116af2dde0c758907083833e2b9b385e7b4c802f01b07df331d.png

../../../_images/8d4f2489c449183031858548919dd22c886d0e502a04390e764c7d2b8571d692.png

../../../_images/ae8fbca102d5e78b6f318a34a02cec4e157d4fd0700a952147faea20165f0b30.png

Pour aller plus loin#

Dérivation et intégration étant incontournables en mathématiques et en physique, il serait étonnant de ne pas trouver ces outils dans certains modules. Le module scipy propose en effet les fonctions derivative et integrate, dont vous pouvez consulter la documentation :

On peut également mentionner le module sympy qui permet de différencier et intégrer de manière formelle.

import sympy as sp
x, y, z = sp.symbols('x y z')

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[21], line 1
----> 1 import sympy as sp
      2 x, y, z = sp.symbols('x y z')

ModuleNotFoundError: No module named 'sympy'

sp.diff(sp.tan(x), x)

sp.integrate(sp.tan(x), x)