Stencil (analyse numérique)

Le stencil de la méthode de Crank-Nicolson en une dimension.

En mathématiques, un stencil est une représentation géométrique d'un réseau nodal illustrant les points d'intérêt utilisés dans un schéma de discrétisation pour la résolution numérique des équations différentielles, notamment des équations aux dérivées partielles combinant variables temporelles et spatiales.

Les stencils peuvent être compacts ou non, selon les niveaux utilisés autour du point d'intérêt.

Étymologie

Les représentations graphiques d'arrangements nodaux et les coefficients associés sont vite apparus dans l’étude des équations aux dérivées partielles. Certains auteurs utilisent des termes différents comme « schémas de relaxation », « instructions opérantes », « logenzes » ou « schémas de points »Modèle:,. Le mot « stencil » a été retenu pour de tels schémas pour refléter l'utilisation d'un stencil sur une grille de calcul afin de repérer uniquement les points utiles.

Calcul de coefficients

Les coefficients de différences finies pour un stencil donné dépendent du choix des points. Ils peuvent s'obtenir par le calcul des dérivées des polynômes de Lagrange interpolant ces points, par résolution d'un système linéaire obtenu à partir d'un développement de Taylor autour de chaque point ou en imposant que le stencil est exact pour chaque monôme jusqu'au degré du schéma. Dans le cas des points équirépartis, on peut également utiliser efficacement l'approximant de Padé de Modèle:Math, avec Modèle:Mvar l'ordre de dérivation recherché et Modèle:Mvar la « distance » entre les nœuds de calcul les plus à gauche sur les deux étages du stencil.

Avec les polynômes interpolateurs de Lagrange

Un moyen d'obtenir les poids des différences finies est de passer par la dérivation du polynôme interpolateur de Lagrange aux nœuds de calcul du stencil choisi. Pour un stencil de longueur n, on note les polynômes d'interpolation :

${\ell }_j(\xi )=\prod _{i=0,i\ne j}^n\frac{\xi -x_i}{x_j-x_i},$

On obtient ainsi un polynôme Modèle:Math de degré n interpolant Modèle:Math en ces n points

$p_n(x)=\sum _{j=0}^{n}f(x_j){\ell }_j(x)$

et de son polynôme dérivé

${p_n}^{\prime }(x)=\sum _{j=0}^{n}f(x_j)\ell {\text{'}}_j(x).$

Ainsi, on a l'approximation de la dérivée première au point Modèle:Mvar avec :

$f^{\prime }(x_c)=\sum _{j=0}^{n}f(x_j)\ell {\text{'}}_j(x_c)+O(h^n).$

Le calcul des coefficients $({{\ell }_j}^{\prime }(x_c))$ permet de retrouver les coefficients obtenus auparavant.

On remarque que la formule est totalement définie sans un pas de discrétisation, ce qui la rend directement applicable aux grilles non uniformes.

Avec les approximants de Padé

Les coefficients de pondération des schémas peuvent être également déterminés par un approximant de Padé, dans le cas où le schéma est semi-implicite.

On considère un stencil de la forme :

$\sum _{k=0}^m{b}_k{f}^{(d)}(x_k)\approx \sum _{j=-s}^n{c}_{j}f(x_j)$

où les points Modèle:Mvar sont sur une grille régulière de pas Modèle:Mvar, soit donc Modèle:Math. On remplace Modèle:Math :

$\sum _{k=0}^m{b}_k(\mathrm{i}\omega {)}^d{\mathrm{e}}^{\mathrm{i}\omega x_k}=\sum _{j=-s}^n{c}_j{\mathrm{e}}^{\mathrm{i}\omega x_j}⟺{\mathrm{e}}^{\mathrm{i}\omega x}(\mathrm{i}\omega {)}^d\sum _{k=0}^m{b}_k{\mathrm{e}}^{\mathrm{i}\omega kh}={\mathrm{e}}^{\mathrm{i}\omega x}\sum _{j=-s}^n{c}_j{\mathrm{e}}^{\mathrm{i}\omega jh}={\mathrm{e}}^{\mathrm{i}\omega x}{\mathrm{e}}^{-\mathrm{i}\omega sh}\sum _{j=0}^{n+s}{c}_{j-s}{\mathrm{e}}^{\mathrm{i}\omega jh}$

On pose Modèle:Math. Le schéma se réécrit alors, par abus de notation :

${\xi }^s{\left(\frac{\ln \xi }{h}\right)}^d=\frac{\sum _{j=0}^n{c}_j{\xi }^j}{\sum _{k=0}^m{b}_k{\xi }^k}.$

Les coefficients optimaux s'obtiennent donc en calculant l'approximant de Padé d'ordre Modèle:Math autour de Modèle:Math.

Compacité du stencil

Un stencil compact 2D utilisant le nœud central (en rouge) et les 8 nœuds adjacents.

Un stencil compact n'utilise que neuf points pour la discrétisation en 2D, en l'occurrence le nœud central et les nœuds adjacents. Ainsi, sur une grille structurée utilisant un stencil compact en dimension 1, 2 ou 3, le nombre maximum de nœuds utilisés est de 3, 9 ou 27 respectivement.

Les stencils compacts sont couramment utilisés et sont au cœur de plusieurs solveurs basant leurs modèles sur des équations aux dérivées partiellesModèle:,

Un stencil non compact 2D.

À l'inverse, un stencil non compact peut faire intervenir tout nœud du voisinage du nœud d'intérêt. Le principal désavantage est que l’utilisation de points éloignés va augmenter le temps de calcul nécessaire à la résolution.

Stencils en une dimension

À deux points

Le stencil à deux points pour la dérivée première d'une fonction est donné par :

$f^{\prime }(x_0)=\frac{f(x_0+h)-f(x_0-h)}{2h}+O\left(h^2\right)$.

qui est obtenu classiquement par un développement limité autour de Modèle:Math :

$\begin{array}{cc}{f}^{\prime }(x_0)& =\frac{f(x_0+h)-f(x_0)}{h}-\frac{f^{(2)}(x_0)}{2!}h-\frac{f^{(3)}(x_0)}{3!}{h}^2-\frac{f^{(4)}(x_0)}{4!}{h}^3+\cdots \\ f^{\prime }(x_0)& =-\frac{f(x_0-h)-f(x_0)}{h}+\frac{f^{(2)}(x_0)}{2!}h-\frac{f^{(3)}(x_0)}{3!}{h}^2+\frac{f^{(4)}(x_0)}{4!}{h}^3+\cdots \end{array}$

À trois points

Le stencil à trois points pour la dérivée seconde d'une fonction est donné par :

$f^{″}(x_0)=\frac{f(x_0+h)-2f(x_0)+f(x_0-h)}{h^2}+O(h^2).$

Ce stencil compact est le plus usuel pour cette dérivée.

À cinq points

Le stencil à 5 points en une dimension s'étend sur la grille

$\{x_0-2h,x_0-h,x_0,x_0+h,x_0+2h\}.$

Le stencil pour la dérivée première est donné par :

$f^{\prime }(x_0)\approx \frac{-f(x_0+2h)+8f(x_0+h)-8f(x_0-h)+f(x_0-2h)}{12h}$

Il est obtenu à partir des développements limités autour de Modèle:Math de Modèle:Math et Modèle:Math jusqu'à l'ordre 3 (5 si on cherche également une estimation de l'erreur) et en recombinant les quatre équations obtenues. On remarque ainsi qu'on a :

$\begin{array}{cc}f(x+h)-f(x-h)& =2h{f}^{\prime }(x)+\frac{h^3}{3}{f}^{(3)}(x)+O(h^4)\\ f(x+2h)-f(x-2h)& =4h{f}^{\prime }(x)+\frac{8h^3}{3}{f}^{(3)}(x)+O(h^4)\end{array}$

dont on tire :

$8(f(x+h)-f(x-h))-(f(x+2h)-f(x-2h))=12h{f}^{\prime }(x)+O(h^4)$

Formellement, on devrait avoir des O(h ⁵) mais on peut l'écrire ainsi sans perte de généralités dans ce cadre.

Le jeu de coefficients (8, -8,-1,1) obtenu ici est un cas spécifique de l'algorithme de Savitzky-Golay.

Estimation de l'erreur

En utilisant les termes d'ordre supérieur (négligés pour l'obtention du stencil), on a une estimation de l'erreur de ce schéma  :

$\frac{-f(x+2h)+8f(x+h)-8f(x-h)+f(x-2h)}{12h}=f^{\prime }(x)-\frac{1}{30}{f}^{(5)}(x)h^4+O(h^5).$

Dérivées d'ordres supérieurs

Sur un schéma à cinq points, on peut obtenir des formules centrées (avec une symétrie dans les coefficients autour du nœud central) pour les dérivées d'ordres supérieurs :

$\begin{array}{cc}{f}^{″}(x)& \approx \frac{-f(x+2h)+16f(x+h)-30f(x)+16f(x-h)-f(x-2h)}{12h^2}+O(h^4)\\ f^{(3)}(x)& \approx \frac{f(x+2h)-2f(x+h)+2f(x-h)-f(x-2h)}{2h^3}+O(h^2)\\ f^{(4)}(x)& \approx \frac{f(x+2h)-4f(x+h)+6f(x)-4f(x-h)+f(x-2h)}{h^4}+O(h^2)\end{array}$

Deux dimensions

Sur une grille régulière 2D, le stencil à 5 points autour du nœud (x, y) est donc donné par :

$\{(x-h,y),(x,y),(x+h,y),(x,y-h),(x,y+h)\},$

Il forme donc un quinconce.

Il apparait dans le calcul du laplacien discret d'une fonction à deux variables:

$\mathrm{\Delta }f(x,y)\approx \frac{f(x-h,y)+f(x+h,y)+f(x,y-h)+f(x,y+h)-4f(x,y)}{h^2}+O(h^2).$

Stencils des schémas numériques usuelles

Les stencils sont représentés graphiquement sur une grille : le nœud Modèle:Math correspond au point Modèle:Math utilisé dans le schéma.

• 
  Schéma d'Euler 1D (décentrage aval)
• 
  Schéma d'Euler 1D (décentrage amont)
• 
  Schéma d'Euler 1D (centré)
• 
  Schéma d'Euler implicite 1D
• 
  Schéma de Crank-Nicholson 1D

Voir aussi

Bibliographie

• Modèle:Ouvrage.
• Communications in Numerical Methods in Engineering, Copyright © 2008 John Wiley & Sons, Ltd.
• Modèle:Handbook of Mathematical Functions (Abramowitz et Stegun). Ninth printing. Table 25.2.

Articles connexes

• Dérivation numérique