Dérivation de la méthode FDTD.

1. Introduction

Cette page présente la démonstration mathématique détaillée des principales équations de mise à jour utilisées dans la méthode des différences finies dans le domaine temporel (FDTD) pour résoudre les équations de Maxwell dans le domaine temporel. L’objectif de cette page est différent de celui de la page consacrée à la théorie générale. Ici, il ne s’agit pas de fournir un aperçu général, mais de montrer explicitement comment les équations continues du rotationnel sont transformées en expressions algébriques discrètes utilisées dans un programme informatique. C’est également la partie de la méthode FDTD qui est la plus susceptible de contenir des erreurs dans le code : des erreurs de signe, des indices mal placés ou un décalage temporel incohérent entraîneront immédiatement des résultats instables ou irréalistes.

Nous partons donc de la loi d’Ampère et de la loi de Faraday sous forme différentielle, nous les développons en composantes cartésiennes, puis nous nous concentrons sur un problème optique à deux dimensions. Les dérivées spatiales sont remplacées par des différences finies centrées, les dérivées temporelles sont discrétisées à l’aide du schéma standard « leapfrog », et les expressions résultantes sont réorganisées afin d’obtenir des formules explicites de mise à jour pour les composantes du champ électrique et du champ magnétique.

Afin de faciliter la compréhension de la dérivation, nous réduisons ici les équations à un cas bidimensionnel dans lequel \(\partial/\partial y = 0\) et la conductivité magnétique sont fixées à zéro, \(\sigma_m = 0\). Cela simplifie considérablement les calculs et montre explicitement comment la réduction dimensionnelle est appliquée en pratique. Cependant, dans le solveur OghmaNano FDTD complet, les équations électromagnétiques sont résolues en trois dimensions. La dérivation bidimensionnelle présentée ici est donc conçue comme un exemple clair et détaillé de la méthode, plutôt que comme une affirmation selon laquelle le solveur sous-jacent est limité à deux dimensions.

Dans tout ce document, les composantes du champ électrique sont évaluées à des intervalles de temps entiers \(n\Delta t\), tandis que les composantes du champ magnétique sont évaluées à des intervalles de temps semi-entiers \(\left(n+\tfrac{1}{2}\right)\Delta t\). Il s’agit de la configuration standard de la méthode de Yee, utilisée dans la méthode FDTD. Contrairement à la page de présentation plus concise, l’indexation est ici explicitement définie afin que la dérivation puisse être comprise à la fois comme un guide d’implémentation et comme une note mathématique.

2. Équations de départ

Nous partons des équations de Maxwell sous forme de rotation dans un milieu linéaire doté d’une conductivité électrique \(\sigma\), d’une permittivité électrique \(\epsilon\), d’une perméabilité magnétique \(\mu\) et, pour généraliser, d’une conductivité magnétique \(\sigma_m\) :

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

Sous forme déterminée, les boucles sont

\[ \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix}, \qquad \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \dfrac{\partial}{\partial x} & \dfrac{\partial}{\partial y} & \dfrac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix}. \]

Il s’agit des équations continues dans le temps et dans l’espace. Le but du reste de cette page est de les transformer en un schéma explicite de différences finies qui peut être utilisé pour effectuer une simulation dans le temps.

3. Forme vectorielle des équations de rotationnel

En développant la loi d’Ampère en ses composantes cartésiennes, on obtient :

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = \frac{\partial H_x}{\partial z} - \frac{\partial H_z}{\partial x}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}. \]

L’extension de la loi de Faraday donne :

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z}, \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}, \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y}. \]

À ce stade, aucune hypothèse n’a encore été formulée concernant la dimensionnalité. Les équations présentées ci-dessus ne sont que la forme composante des équations complètes du rotationnel en trois dimensions.

4. Réduction bidimensionnelle

Afin de simplifier le processus de dérivation, nous réduisons maintenant le problème à deux dimensions en supposant qu’il n’y a aucune variation dans la direction \(y\), ce qui implique \(\partial / \partial y = 0\). Cette étape est effectuée uniquement pour simplifier les calculs et faciliter la compréhension de la structure des équations de mise à jour. La même procédure peut être directement étendue à trois dimensions, et dans la pratique, OghmaNano résout la forme complète en 3D des équations de Maxwell.

En partant de cette hypothèse, les équations du champ électrique se simplifient en :

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}, \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}, \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

Si nous fixons également la conductivité magnétique à zéro, \(\sigma_m = 0\), les équations du champ magnétique deviennent :

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}, \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}, \]

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

Voici les équations continues qui seront désormais discrétisées. Il est important de préciser qu’il s’agit d’une réduction en 2D du problème complet en 3D. L’objectif n’est pas d’affirmer que toutes les méthodes FDTD sont intrinsèquement bidimensionnelles, mais de dériver des équations de mise à jour claires pour le cas courant où une direction est supprimée.

5. Grille de Yee, indexation et notation.

Pour résoudre numériquement les équations ci-dessus, l’espace et le temps sont discrétisés. Soit \(\Delta x\) et \(\Delta z\) les pas de la grille de calcul dans les directions \(x\) et \(z\), respectivement, et soit \(\Delta t\) le pas de temps. Nous utilisons les indices entiers \(i\) et \(k\) pour les directions \(x\) et \(z\) :

\[ x_i = i \Delta x, \qquad z_k = k \Delta z. \]

Le temps est indexé par l’entier \(n\), de sorte que

\[ t^n = n \Delta t. \]

Dans le schéma de Yee, les champs électriques et magnétiques ne sont pas stockés au même endroit, ni dans l’espace, ni dans le temps. Au lieu de cela, ils sont décalés. Les composantes du champ électrique sont calculées à des instants discrets \(n\Delta t\), tandis que les composantes du champ magnétique sont calculées à des instants semi-discrets \(\left(n+\tfrac{1}{2}\right)\Delta t\). C’est là l’origine de la structure en « pas de grenouille » de la méthode FDTD.

Sur le plan spatial, les composantes du champ sont également décalées de manière à ce que chaque opérateur rotationnel aux différences finies, centré, soit naturellement défini entre les points de grille voisins. Dans une implémentation complète du code, chaque composante du champ est associée à sa propre paire d’indices spatiaux, par exemple \(E_x^n(i,k)\), \(H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2})\), et ainsi de suite. Afin de faciliter la lecture des équations, certains des indices « spectateurs » répétés ne sont pas toujours écrits explicitement dans les étapes suivantes. Cependant, l’interprétation reste la même :

\(i\) indexe la direction \(x\).
\(k\) indexe la direction \(z\).
\(n\) représente les niveaux temporels du champ électrique, exprimés en valeurs entières, et
\(n+\tfrac{1}{2}\) indique les niveaux temporels du champ magnétique, avec des incréments d’un demi-pas.

Par exemple, l’approximation centrée de la méthode des différences finies pour \(\partial H_y / \partial z\) est la suivante :

\[ \frac{\partial H_y}{\partial z} \approx \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}, \]

et de même, l’approximation centrée de \(\partial H_z / \partial x\) est :

\[ \frac{\partial H_z}{\partial x} \approx \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Une fois cette notation clairement définie, la dérivation qui suit peut être interprétée sans ambiguïté comme un schéma de mise à jour correct pour une grille de Yee.

6. Dérivation des équations de mise à jour du champ électrique.

Le champ électrique est calculé à des intervalles de temps entiers \(n\Delta t\), tandis que le champ magnétique est calculé à des intervalles de temps semi-entiers \(\left(n+\tfrac{1}{2}\right)\Delta t\). Pour établir les équations de mise à jour, la dérivée temporelle est approchée à l’aide d’une différence finie centrée, et le terme de conductivité est calculé en utilisant la moyenne des champs électriques anciens et nouveaux. Cela permet d’obtenir une méthode numérique stable, de type trapézoïdale, pour le calcul de la conductivité.

Cette étape de calcul de la moyenne est importante. Si le terme de conductivité était traité de manière incohérente, les coefficients de mise à jour ne présenteraient plus la forme amortie standard, et la correspondance avec les équations FDTD habituelles serait perdue.

6.1 Mettre à jour l’équation de \(E_x\).

Commencez par

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = - \frac{\partial H_y}{\partial z}. \]

Évaluez le terme de conductivité en utilisant la moyenne de \(E_x^{n+1}\) et \(E_x^n\), et évaluez la dérivée temporelle en utilisant une méthode aux différences finies centrée :

\[ \sigma \frac{E_x^{n+1}(i,k) + E_x^n(i,k)}{2} + \epsilon \frac{E_x^{n+1}(i,k) - E_x^n(i,k)}{\Delta t} = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Regroupez les termes \(E_x^{n+1}\) à gauche :

\[ \left( \frac{\sigma}{2} + \frac{\epsilon}{\Delta t} \right) E_x^{n+1}(i,k) = - \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} + \left( \frac{\epsilon}{\Delta t} - \frac{\sigma}{2} \right) E_x^n(i,k). \]

Par conséquent.

\[ E_x^{n+1}(i,k) = \frac{ \left( \dfrac{\epsilon}{\Delta t} - \dfrac{\sigma}{2} \right) E_x^n(i,k) - \dfrac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} }{ \dfrac{\epsilon}{\Delta t} + \dfrac{\sigma}{2} }. \]

Si l’on multiplie le numérateur et le dénominateur par \(2\Delta t\), on obtient la forme du coefficient que l’on connaît bien :

\[ E_x^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_x^n(i,k) - \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

6.2 Mettre à jour l’équation pour \(E_y\).

Commencez par

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = - \frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z}. \]

Discrétisez les deux termes de dérivée sur la grille de Yee :

\[ \sigma \frac{E_y^{n+1}(i,k) + E_y^n(i,k)}{2} + \epsilon \frac{E_y^{n+1}(i,k) - E_y^n(i,k)}{\Delta t} = - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z}. \]

Réorganiser permet de

\[ E_y^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_y^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \left[ - \frac{ H_z^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_z^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x} + \frac{ H_x^{n+\frac{1}{2}}\!\left(i,k+\frac{1}{2}\right) - H_x^{n+\frac{1}{2}}\!\left(i,k-\frac{1}{2}\right) }{\Delta z} \right]. \]

Il s’agit de la seule composante du champ électrique dans la réduction 2D qui contient des contributions provenant de deux dérivées distinctes du champ magnétique, et c’est donc la plus susceptible de provoquer des erreurs dans le code. Plus précisément, le signe relatif entre les termes de dérivée \(x\) et \(z\) découle directement du développement du rotationnel et ne doit pas être modifié.

6.3 Mettre à jour l’équation pour \(E_z\).

Commencez par

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x}. \]

La discrétisation permet de

\[ \sigma \frac{E_z^{n+1}(i,k) + E_z^n(i,k)}{2} + \epsilon \frac{E_z^{n+1}(i,k) - E_z^n(i,k)}{\Delta t} = \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

Par conséquent.

\[ E_z^{n+1}(i,k) = \frac{ 2\epsilon - \sigma \Delta t }{ 2\epsilon + \sigma \Delta t } E_z^n(i,k) + \frac{ 2 \Delta t }{ 2\epsilon + \sigma \Delta t } \frac{ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2},k\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2},k\right) }{\Delta x}. \]

À ce stade, les équations discrètes mises à jour pour le champ électrique ont été obtenues sous leur forme standard, avec amortissement. L’étape suivante consiste à dériver les équations correspondantes pour le champ magnétique, en utilisant des pas de temps semi-entiers.

7. Dérivation des équations de mise à jour du champ magnétique.

Le champ magnétique est mis à jour à des intervalles de temps semi-entiers. Avec \(\sigma_m = 0\), les équations du champ magnétique ne contiennent aucun terme d’amortissement, de sorte que les formules de mise à jour prennent une forme plus simple, de type « pas de temps avant ». Cette partie de la dérivation est algébriquement plus courte que dans le cas du champ électrique, car il n’y a pas de terme de conductivité trapézoïdale à prendre en compte.

7.1 Mettre à jour l’équation de \(H_x\).

Commencez par

\[ \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_y}{\partial z}. \]

En effectuant une discrétisation temporelle entre les étapes intermédiaires \(n-\tfrac{1}{2}\) et \(n+\tfrac{1}{2}\), et une discrétisation spatiale en utilisant une différence centrée dans la direction \(z\), on obtient :

\[ \mu \frac{ H_x^{n+\frac{1}{2}}(i,k+\frac{1}{2}) - H_x^{n-\frac{1}{2}}(i,k+\frac{1}{2}) }{\Delta t} = \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}. \]

Dans la notation compacte à demi-indice, souvent utilisée dans les démonstrations, le même résultat peut s’écrire comme suit :

\[ \mu \frac{ H_x^{n+\frac{1}{2}} - H_x^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Par conséquent.

\[ H_x^{n+\frac{1}{2}} = H_x^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(k+\frac{1}{2}\right) - E_y^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

7.2 Mettre à jour l’équation pour \(H_y\).

Commencez par

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} - \frac{\partial E_x}{\partial z}. \]

La discrétisation permet de

\[ \mu \frac{ H_y^{n+\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) - H_y^{n-\frac{1}{2}}(i+\frac{1}{2},k+\frac{1}{2}) }{\Delta t} = \frac{ E_z^n(i+1,k+\frac{1}{2}) - E_z^n(i,k+\frac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\frac{1}{2},k+1) - E_x^n(i+\frac{1}{2},k) }{\Delta z}. \]

Encore une fois, en utilisant la notation abrégée de l’indice partiel,

\[ \mu \frac{ H_y^{n+\frac{1}{2}} - H_y^{n-\frac{1}{2}} }{\Delta t} = \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z}. \]

Par conséquent.

\[ H_y^{n+\frac{1}{2}} = H_y^{n-\frac{1}{2}} + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n\!\left(i+\frac{1}{2}\right) - E_z^n\!\left(i-\frac{1}{2}\right) }{\Delta x} - \frac{ E_x^n\!\left(k+\frac{1}{2}\right) - E_x^n\!\left(k-\frac{1}{2}\right) }{\Delta z} \right]. \]

7.3 Mettre à jour l’équation pour \(H_z\).

Commencez par

\[ \mu \frac{\partial H_z}{\partial t} = - \frac{\partial E_y}{\partial x}. \]

La discrétisation permet de

\[ \mu \frac{ H_z^{n+\frac{1}{2}}(i+\frac{1}{2},k) - H_z^{n-\frac{1}{2}}(i+\frac{1}{2},k) }{\Delta t} = - \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

De manière équivalente, en utilisant une notation plus concise,

\[ \mu \frac{ H_z^{n+\frac{1}{2}} - H_z^{n-\frac{1}{2}} }{\Delta t} = - \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Par conséquent.

\[ H_z^{n+\frac{1}{2}} = H_z^{n-\frac{1}{2}} - \frac{\Delta t}{\mu} \frac{ E_y^n\!\left(i+\frac{1}{2}\right) - E_y^n\!\left(i-\frac{1}{2}\right) }{\Delta x}. \]

Une fois ces substitutions effectuées, les équations standard de mise à jour de la méthode FDTD (méthode de différences finies dans le domaine temporel) peuvent être appliquées directement.

8. Équations finales de mise à jour.

Une fois les résultats recueillis, les équations de mise à jour bidimensionnelles de la méthode FDTD sont les suivantes :

\[ H_x^{n+\frac{1}{2}}(i,k+\tfrac{1}{2}) = H_x^{n-\frac{1}{2}}(i,k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \frac{ E_y^n(i,k+1) - E_y^n(i,k) }{\Delta z}, \]

\[ H_y^{n+\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) = H_y^{n-\frac{1}{2}}(i+\tfrac{1}{2},k+\tfrac{1}{2}) + \frac{\Delta t}{\mu} \left[ \frac{ E_z^n(i+1,k+\tfrac{1}{2}) - E_z^n(i,k+\tfrac{1}{2}) }{\Delta x} - \frac{ E_x^n(i+\tfrac{1}{2},k+1) - E_x^n(i+\tfrac{1}{2},k) }{\Delta z} \right], \]

\[ H_z^{n+\frac{1}{2}}(i+\tfrac{1}{2},k) = H_z^{n-\frac{1}{2}}(i+\tfrac{1}{2},k) - \frac{\Delta t}{\mu} \frac{ E_y^n(i+1,k) - E_y^n(i,k) }{\Delta x}. \]

Sous cette forme, la signification du schéma de type « saut de puce » est clairement définie. Le champ électrique au niveau de temps \(n+1\) est obtenu à partir du champ magnétique au niveau de temps \(n+\tfrac{1}{2}\), tandis que le champ magnétique au niveau de temps \(n+\tfrac{1}{2}\) est obtenu à partir du champ électrique au niveau de temps \(n\). Les champs évoluent donc selon une séquence alternée dans le temps.

9. Observations sur la mise en œuvre.

Dans un code informatique concret, les coefficients qui multiplient les anciennes valeurs du champ et les termes de rotation sont souvent précalculés et stockés, en particulier lorsque \(\epsilon\), \(\mu\) ou \(\sigma\) varient d’une cellule à l’autre. Par exemple, la mise à jour du champ électrique prend généralement la forme suivante :

\[ E^{n+1} = C_a E^n + C_b (\nabla \times H), \]

où les coefficients locaux \(C_a\) et \(C_b\) représentent les paramètres du matériau et le pas de temps. Cela permet d’obtenir une implémentation plus claire et d’éviter de recalculer les mêmes facteurs à chaque pas de temps.

Il convient également de souligner que la formule présentée ci-dessus ne donne que le schéma de base de la méthode FDTD. Une simulation utilisable nécessite en outre des sources, des conditions aux limites, des couches absorbantes telles que les couches parfaitement adaptées (PML), et un pas de temps qui satisfait la condition de stabilité de Courant. Aucune de ces additions ne modifie la structure fondamentale du schéma « leapfrog » que nous avons dérivée, mais chacune d’entre elles modifie la manière dont les équations de mise à jour sont appliquées à des points spécifiques de la grille.

Enfin, bien que la notation ci-dessus ait été rédigée avec suffisamment de soin pour être pertinente en termes d’implémentation, la disposition spatiale exacte de chaque composant dépend de la convention de grille utilisée dans le code. L’exigence essentielle n’est pas la typographie elle-même, mais la cohérence : l’opérateur de rotation, le décalage et les signes doivent tous correspondre à la même disposition de cellules de Yee dans tout le programme.