-- Simulation numérique en mécanique des fluides --

TP n°4

Dérivation numérique

simon.marie@lecnam.net
In [1]:
import numpy as np
import matplotlib.pyplot as plt
# la librairie misc de scipy permet de convertir les images en objets numpy
import scipy.misc as scm
# la librairie time nous sera utile pour calculer le temps d'execution des scripts
import time 
#Option pour afficher les figures dans le notebook et eviter le plt.show():
%matplotlib inline   

1 - Derivation numérique d'une fonction analytique

On se propose ici de tester un schéma aux différences finies pour calculer la dérivée d'une fonction sinus $f(x)=\sin(kx)$ soit (comme chacun sait) $f'(x)=k\cos(kx)$:

In [13]:
# On définie une première fonction qui va calculer une approximation numérique de la dérivée avec un schéma
# décentré d'ordre 1
def diff(f,dx):
    df=np.zeros(f.shape)
    # on peut l'écrire avec une boucle for
    for i in range(len(f)-2): 
        df[i]= (f[i+1]-f[i])/dx # ordre 1 avec boucle
    #df[0:-1]=(f[1:]-f[0:-1])/dx # ordre 1 avec slicing
    #df[0:-2]=... # ordre 2
    return df

# On définie une autre fonction qui va calculer une approximation numérique de la dérivée avec un schéma
# centré d'ordre 2
def diffc(f,dx):
    df=np.zeros(f.shape)
    for i in range(1,len(f)-1):
        df[i]= (f[i+1]-f[i-1])/(2.*dx) # ordre 1
    #df[1:-1]=... # ordre 2
    return df
In [14]:
dx=0.05
x=np.arange(-1,1,dx)
y=np.sin(10.*x)
dy_th=10.*np.cos(10.*x)

dy=diff(y,dx)
dy2=diffc(y,dx)
In [16]:
plt.figure(figsize=(10,5))
plt.plot(x,y,'k',label='$f(x)$')
plt.plot(x,dy,'ro',label='ordre $1$ décentré')
plt.plot(x,dy2,'gs',label='ordre $2$ centré')
plt.plot(x,dy_th,'b',label='$f\'(x)$')
plt.legend()
Out[16]:
<matplotlib.legend.Legend at 0x7ff455396390>

Le schéma d'ordre 2 est plus précis lorsque $\Delta x$ est faible. On pourra essayer de changer les valeurs de $\Delta x$.

2 - Dérivation 2D

On se propose dans cette partie de définir des fonctions permetant de calculer le gradient et le laplacien d'une quantité $f$ donnée, dépendant de deux variables $x$ et $y$. On rapelle:

$$ \mathbf{grad}.x= \left( \begin{array}{c} \dfrac{\partial f}{\partial x} \\ \dfrac{\partial f}{\partial y} \\ \end{array} \right) $$

$$ \Delta f = \dfrac{\partial^2 f}{\partial x^2}+\dfrac{\partial^2 f}{\partial y^2} $$

On utilisera dans ces fonctions les différents schémas d'espace vues en cours.

In [46]:
# On construit la fonction grad avec 3 argument, 1 obligatoire: f
# et deux facultatifs o et dx qui sont l'ordre souhaité et la taille de maille.
def grad(f,o=1,dx=1.):
    #On note dfx la composante suivant x du gradient:
    dfx=np.zeros((f.shape[0],f.shape[1])) #on initialise dfx a 0
    #On note dfy la composante suivant y du gradient:
    dfy=np.zeros((f.shape[0],f.shape[1])) #on initialise dfy a 0
    # Calcul de dfx et dfy en fonction de l'ordre:
    if o==1: # ordre 1 decentre
        dfx[0:-1,:]=(f[1:,:]-f[0:-1,:])/dx
        dfy[:,0:-1]=(f[:,1:]-f[:,0:-1])/dx
    elif o==2: # ordre 2 centre
        dfx[1:-1,:]=(f[2:,:]-f[0:-2,:])/dx
        dfy[:,1:-1]=(f[:,2:]-f[:,0:-2])/dx
    else:
        print('o < 2 ')
    return dfx,dfy


def lap(f,dx=1.):
    Lx=np.zeros((f.shape[0],f.shape[1])) #on initialise Lx a 0
    Ly=np.zeros((f.shape[0],f.shape[1])) #on initialise Ly a 0
    Lx[1:-1,:]=(f[0:-2,:]-2.*f[1:-1,:]+f[2:,:])/dx
    Ly[:,1:-1]=(f[:,0:-2]-2.*f[:,1:-1]+f[:,2:])/dx
    return Lx+Ly

Interprétation visuelle

Pour bien illustrer la notion de dérivation numérique, on peut utiliser une image en niveau de gris et appliquer les opérateur de dérivation sur cette image. On pourra utiliser l' image suivante (attention à bien l'enregistrer en local si vous executez le notebook):

In [52]:
f=scm.imread('CNAM.jpg')
In [53]:
plt.imshow(f,cmap='gray')
print('min/max:',f.min(),f.max())
print('taille:',f.shape)
min/max: 0 255
taille: (1099, 736)
In [54]:
fig=plt.figure(figsize=(15,10))
dfx,dfy=grad(f,2)
plt.imshow(dfx+dfy,cmap='Greys');plt.colorbar()
Out[54]:
<matplotlib.colorbar.Colorbar at 0x7ff44ea79048>
In [55]:
fig=plt.figure(figsize=(10,5))
fig.add_subplot(121)
dfx,dfy=grad(f,2)
plt.imshow((dfx+dfy)[200:550,300:500],cmap='Greys',vmin=-100,vmax=100);plt.colorbar()
fig.add_subplot(122)
Lap=lap(f)
plt.imshow(Lap[200:550,300:500],cmap='Greys',vmin=-100,vmax=100);plt.colorbar()
Out[55]:
<matplotlib.colorbar.Colorbar at 0x7ff44c466390>