# Points de divergence négative dans l'équation de diffusion pour la modélisation de capture des moustiques

Aya Chaieb, Baptiste Petiot, Louis Mudarra, Lucas Bourret, Facundo Muñoz

► To cite this version:

Aya Chaieb, Baptiste Petiot, Louis Mudarra, Lucas Bourret, Facundo Muñoz. Points de divergence négative dans l'équation de diffusion pour la modélisation de capture des moustiques. Cirad. 2025, pp.21. hal-05084285

HAL Id: hal-05084285

https://hal.inrae.fr/hal-05084285v1

Submitted on 26 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Points de divergence négative dans l'équation de
diffusion pour la modélisation de capture des
moustiques

Aya Chaieb†, Baptiste Petiot†, Louis Mudarra†, Lucas Bourret†,
Facundo Muñoz‡ ID

25 mai 2025

L'étude des méthodes de réduction d'une population de moustiques, et notamment de la technique SIT (sterile insect technique), repose sur un modèle de diffusion spatiale formulé par une EDP de type équation de la chaleur qui gouverne le comportement de la densité surfacique de moustiques. L'action de pièges à moustiques est intégrée au modèle via des termes puits localisés (points de divergence négative), puis la résolution numérique de l'EDP ainsi obtenue est détaillée. Les méthodes étudiées sont la méthode des différences finies (MDF), la méthode des éléments finis (MEF) et la méthode des volumes finis (MVF), dont la formulation mathématique, la discrétisation spatiale et temporelle ainsi que l'algorithme de résolution sont précisés. Les résultats obtenus sont concluants pour la MDF et la MEF, tandis que des difficultés persistent pour la MVF. Une fonction Python capable de résoudre l'équation enrichie et de renvoyer, pour une surface donnée, le nombre de moustiques est finalement développée, et son implémentation est détaillée. L'accent est mis sur la complexité temporelle, la précision et la robustesse de cette fonction de résolution. Ce travail s'inscrit dans une étude plus large menée par le CIRAD, qui vise à estimer les paramètres de diffusion.

<figure><img src="image_2.png" /></figure>

CentraleSupélec

<footnote>† Étudiants-chercheurs du programme ingénieur de CentraleSupélec, Université Paris-Saclay.</footnote>

<footnote>‡ CIRAD, UMR ASTRE, F-34398 Montpellier, France.
ASTRE, Univ Montpellier, CIRAD, INRAE, Montpellier, France.</footnote>
Table des matières

<table><tr><td>1 Introduction</td><td>3</td></tr><tr><td>1.1 Contexte global de l'étude</td><td>3</td></tr><tr><td>1.2 Objectif et méthodologie</td><td>3</td></tr><tr><td>1.3 Cadre mathématique</td><td>4</td></tr><tr><td>2 Étude préliminaire de l'EDP</td><td>5</td></tr><tr><td>2.1 Formulation faible du problème</td><td>5</td></tr><tr><td>2.2 Caractère bien posé du problème</td><td>6</td></tr><tr><td>3 Modélisation des pièges – Choix des termes puits</td><td>7</td></tr><tr><td>3.1 Terme puits ponctuel</td><td>7</td></tr><tr><td>3.2 Terme puits constant</td><td>8</td></tr><tr><td>3.3 Terme puits gaussien</td><td>8</td></tr><tr><td>3.4 Modélisation multi-pièges</td><td>9</td></tr><tr><td>4 Méthode des différences finies</td><td>9</td></tr><tr><td>4.1 Définition des maillages</td><td>10</td></tr><tr><td>4.2 Définition du schéma de récurrence</td><td>10</td></tr><tr><td>4.3 Complexité et précision</td><td>10</td></tr><tr><td>4.4 Algorithme et résultats</td><td>10</td></tr><tr><td>5 Méthode des éléments finis</td><td>12</td></tr><tr><td>5.1 Définition des éléments finis P1</td><td>12</td></tr><tr><td>5.2 Définition du maillage</td><td>12</td></tr><tr><td>5.3 Étude de l'EDO</td><td>13</td></tr><tr><td>5.4 Complexité et précision</td><td>14</td></tr><tr><td>5.5 Algorithme et résultats</td><td>14</td></tr><tr><td>6 Méthode des volumes finis</td><td>15</td></tr><tr><td>6.1 Définition du maillage</td><td>15</td></tr><tr><td>6.2 Discrétisation en espace</td><td>15</td></tr><tr><td>6.3 Discrétisation en temps</td><td>16</td></tr><tr><td>6.4 Algorithme et résultats</td><td>17</td></tr><tr><td>7 Comparaison des algorithmes</td><td>18</td></tr><tr><td>8 Conclusion</td><td>20</td></tr></table>
# 1 Introduction

## 1.1 Contexte global de l'étude

Les moustiques peuvent jouer un rôle important dans l'équilibre de certains écosystèmes, mais dans d'autres cas, ils constituent un vecteur épidémiologique préoccupant, contribuant à la propagation de maladies graves pour l'être humain, telles que la dengue, le virus Zika ou le chikungunya. On cherche donc à réguler certaines populations, voire à les éradiquer localement dans des contextes ciblés. L'utilisation d'insecticides a toutefois montré ses limites, notamment en raison de l'apparition rapide de mutations adaptatives chez les moustiques, comme la mutation *kdr*, qui confère une résistance accrue en seulement quelques générations. Par ailleurs, dès les années 1940-1950, une alternative a émergé avec la méthode dite SIT (Sterile Insect Technique, ou technique de l'insecte stérile), qui a déjà démontré son efficacité (NURHAYATI et al., 2013). De nombreuses autres études expérimentales et campagnes de terrain à travers le monde ont confirmé le potentiel de cette méthode dans des contextes variés. Cette technique repose sur la stérilisation des moustiques via l'utilisation de rayons X ou rayons γ (voir figure 1).

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Stérilisation des moustiques (NURHAYATI et al., 2013)</figcaption></figure>

Dans le cadre de la mise en œuvre de cette technique de l'insecte stérile (SIT), la connaissance fine des dynamiques spatiales des moustiques constitue un enjeu central (voir figure 2). En particulier, la vitesse de diffusion des mâles stériles relâchés ainsi que leur survie conditionnent la fréquence et la densité des sessions de lâcher à grande échelle pour garantir l'efficacité de la stratégie. Afin d'estimer le paramètre d'intérêt de l'équation de diffusion, il est nécessaire de mener des études expérimentales de comptage des moustiques à l'aide de pièges, ce qui permet de remonter indirectement aux paramètres de l'équation de diffusion par une analyse statistique. Cependant, l'introduction de ces pièges modifie significativement le modèle. C'est pourquoi dans le cadre précis de notre étude, on s'intéresse à la modélisation de la libération des moustiques stériles en présence de pièges.

## 1.2 Objectif et méthodologie

L'objectif de cet article est de proposer une approche mathématique permettant de modéliser l'effet de pièges à moustiques dans un cadre de diffusion bidimensionnelle. Pour cela, des puits locaux affectant la densité dans une équation aux dérivées partielles de type diffusion sont introduits, représentant les effets de mortalité induits par des pièges répartis spatialement. Un code Python permettant une
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Simulation de la dispersion des moustiques stériles relâchés (MUÑOZ, 2023)</figcaption></figure>

résolution efficace de cette nouvelle EDP est développé en parallèle. Cette résolution numérique sera
utilisée par la suite pour estimer les paramètres d'intérêt dans le cadre d'une analyse statistique du
système.

Cet article se concentre sur la résolution mathématique et numérique de cette équation, en comparant plusieurs méthodes classiques : différences finies, éléments finis, et volumes finis. Nous cherchons à déterminer quelle méthode numérique offre la meilleure précision et le meilleur temps de calcul pour l'équation aux dérivées partielles issue de notre modélisation. Une attention particulière est portée au temps de calcul de chacune de ces approches, dans le but d'identifier la méthode la plus adaptée à un usage intensif en simulation inverse ou en inférence bayésienne. Ce travail s'inscrit dans une démarche plus large de modélisation multiscalaire visant à fournir des outils prédictifs robustes pour les stratégies de lutte anti-vectorielle.

1.3 Cadre mathématique

Cette sous-section vise à préciser le cadre mathématique et les principales équations de l'étude, et
aboutit à l'équation de diffusion enrichie qu'il convient de résoudre.

On se place dans le plan $\mathbb{R}^2$ et on repère la position d'un point par $\mathbf{r} = (x, y) \in \mathbb{R}^2$. On définit

$$
\rho : \mathbb{R}^2 \times \mathbb{R}_+ \rightarrow \mathbb{R}_+ \\
(\mathbf{r}, t) \mapsto \rho(\mathbf{r}, t)
$$

représentant la densité de moustiques (en moustiques · m⁻²). Le flux de moustiques est défini par

$$
\mathbf{j} : \mathbb{R}^2 \times \mathbb{R}_+ \rightarrow \mathbb{R}^2 \\
(\mathbf{r}, t) \mapsto \mathbf{j}(\mathbf{r}, t)
$$

et se mesure en moustiques · m⁻¹ · s⁻¹. Le lien avec la densité de moustique ρ est donné par la première
loi de Fick (équation 2). Le coefficient de diffusion D est une constante (en m² · s⁻¹) et le terme source
K, défini par

$$
K : \mathbb{R}^2 \to \mathbb{R}_+ \\
\mathbf{r} \mapsto K(\mathbf{r})
$$

représente ici un puits (en s⁻¹). On note $K(\mathbf{r}) = M_0 + S(\mathbf{r})$, où $M_0$ est le taux de mort naturel des moustiques, et $S(\mathbf{r})$ représente la diminution locale de la densité de moustiques relative aux pièges. L'équation de conservation de la densité de moustiques s'écrit :

$$
\partial_t \rho(\mathbf{r}, t) + \nabla \cdot \mathbf{j} = -\rho(\mathbf{r}, t) K(\mathbf{r}) \quad (1)
$$
On écrit le terme source directement comme produit de $S$ et de $\rho$ afin de simplifier les notations, de cette manière les unités apparaissent plus clairement. La loi de Fick, qui relie le flux à la densité, se formule :

$$ \mathbf{j} = -D \nabla \rho \qquad (2) $$

En injectant l'équation 2 dans l'équation 1, on obtient l'équation de diffusion :

$$ \partial_t \rho(\mathbf{r}, t)(\mathbf{r}, t) = D \Delta \rho(\mathbf{r}, t) - \rho(\mathbf{r}, t) K(\mathbf{r}) \qquad (3) $$

Partant de l'équation 3, en considérant les conditions aux bords spatiales et temporelles, on a l'équation aux dérivées partielles :

$$ \left\{ \begin{array}{ll} \partial_t \rho(\mathbf{r}, t) = D \Delta \rho(\mathbf{r}, t) - \rho(\mathbf{r}, t) K(\mathbf{r}) & , \mathbf{r} \in \Omega \\ \rho(\mathbf{r}, t) = 0 & , \mathbf{r} \in \partial \Omega \\ \rho(\mathbf{r}, 0) = \rho_0 \mathbb{1}_{B(0, \varepsilon)} & , \rho_0 \in \mathbb{R}_+^* \end{array} \right. \qquad (4) $$

Ici, $\rho_0$ désigne la densité initiale de moustiques lors de leur relâche, densité cruciale puisqu'elle constitue le point de départ de leur diffusion. Dans notre cas, la libération s'effectue au centre du domaine, soit dans la boule ouverte $B(0, \varepsilon)$, où $\varepsilon > 0$.

La section 2 se concentre sur l'étude théorie de l'EDP obtenue, en particulier, les outils d'analyse fonctionnelle nécessaires sont présentés, et le caractère bien posé du problème est démontré. La section 3 présente différents termes puits $S(\mathbf{r})$ possibles pour cette modélisation, ainsi que leurs avantages respectifs. Cette section aboutit au choix du piège gaussien (sous-section 3.3). Les sections suivantes plongent dans la résolution numérique et mathématique de l'équation 4. La section 4 détaille les schémas numériques utilisés pour la méthode des différences finies. La section 5 présente la formulation faible du problème et permet l'obtention d'un système linéaire par la méthode des éléments finis. Enfin, la section 6 décrit l'application de la méthode des volumes finis pour aboutir à un système linéaire.

## 2 Étude préliminaire de l'EDP

L'étude présente dans cette section sert à définir les éléments d'analyse fonctionnelle nécessaire pour la suite, ainsi qu'à s'assurer du caractère bien posé du problème. Il convient alors de démontrer d'abord l'existence et l'unicité d'une solution faible, positive et continue par rapport aux données initiales, puis d'établir la formulation variationnelle dans $H_0^1(\Omega)$ qui servira de base aux discrétisations numériques ultérieures (en particulier pour la méthode des éléments finis, voir section 5).

### 2.1 Formulation faible du problème

On considère le problème équation 4, modélisant la diffusion spatiale des moustiques en présence de puits locaux (pièges), où $\Omega \subset \mathbb{R}^2$ est un domaine borné sur le bord duquel on impose une condition de Dirichlet homogène. On suppose dans cette section que $K(\mathbf{r})$ est suffisamment régulier sur $\Omega$.

Pour des données initiales peu régulières, il n'est pas garanti que $\rho$ admette des dérivées classiques en chaque point. Toutefois, on peut définir $\rho$ comme solution dans un sens faible (ou distributionnel). Cette approche, fondée sur les outils d'analyse fonctionnelle (espaces de Sobolev, théorie de l'intégration, etc.), permet non seulement d'établir l'existence de solutions en toute rigueur, mais aussi de préparer la résolution numérique, notamment par la méthode des éléments finis (voir section 5).

Afin d'établir la formulation variationnelle du problème, ainsi que son caractère bien posé, nous adoptons un cadre fonctionnel approprié. Notamment, pour chaque $t$, la fonction $\mathbf{r} \mapsto \rho(\mathbf{r}, t)$ appartient à l'espace de Sobolev $H_0^1(\Omega)$, qui est un sous-espace de $H^1(\Omega)$. On définit alors ces espaces de la manière suivante (pour plus de détails, voir BREZIS, 2011) :

**Définition 1 (Espace de Sobolev $H^1(\Omega)$)** Soit $\Omega \subset \mathbb{R}^d$ un ouvert borné. On définit

$$ H^1(\Omega) := \{ u \in L^2(\Omega) \mid \partial_{x_i} u \in L^2(\Omega) \text{ pour } i = 1, \dots, d \}, $$

où $\partial_{x_i} u$ désigne la dérivée faible (au sens des distributions) de $u$ suivant $x_i$.
**Définition 2 (Espace $H_0^1(\Omega)$)** Soit $\gamma : H^1(\Omega) \to L^2(\partial\Omega)$ l'opérateur de trace. On définit

$$H_0^1(\Omega) := \{ u \in H^1(\Omega) \mid \gamma(u) = 0 \text{ sur } \partial\Omega \},$$

c'est-à-dire l'ensemble des fonctions de $H^1(\Omega)$ dont la valeur au bord s'annule.

**Définition 3 (Espace de Sobolev négatif $H^{-1}(\Omega)$)** Pour un ouvert borné $\Omega \subset \mathbb{R}^d$, on définit

$$H^{-1}(\Omega) := (H_0^1(\Omega))^*$$

c'est-à-dire le dual topologique de $H_0^1(\Omega)$. Un élément $f \in H^{-1}(\Omega)$ est donc une forme linéaire continue $f : H_0^1(\Omega) \to \mathbb{R}$.

On rappelle que $H_0^1(\Omega)$ est un espace de Hilbert pour la norme $\|u\|_{H_0^1} = \|\nabla u\|_{L^2}$, à équivalence près, d'après l'inégalité de Poincaré (BREZIS, 2011). Le terme $D\Delta\rho$ est interprété au sens des distributions, et la condition aux limites $\rho = 0$ sur $\partial\Omega$ garantit que le terme de diffusion ne produit pas de flux sortant.

Nous allons donc écrire la formulation faible du problème et en donner la définition mathématique précise, en suivant notamment EVANS, 2010 et LIONS et MAGENES, 1968.

Pour obtenir la formulation faible du problème posé par l'équation 4, on multiplie les termes de l'équation par une fonction test $\omega \in H_0^1(\Omega)$ et on intègre sur $\Omega$. Il vient :

$$\int_{\Omega} \left( \frac{\partial \rho}{\partial t}(\mathbf{r}, t) - D \Delta \rho(\mathbf{r}, t) + \rho(\mathbf{r}, t) K(\mathbf{r}) \right) \omega(\mathbf{r}) \, d\mathbf{r} = 0 \quad (5)$$

Par la formule de Green et le théorème de dérivation sous le signe intégral, on obtient :

$$\frac{d}{dt} \int_{\Omega} \rho(\mathbf{r}, t) \omega(\mathbf{r}) \, d\mathbf{r} + \int_{\Omega} D \nabla_{\mathbf{r}} \rho(\mathbf{r}, t) \cdot \nabla_{\mathbf{r}} \omega(\mathbf{r}) \, d\mathbf{r} + \int_{\Omega} \rho(\mathbf{r}, t) K(\mathbf{r}) \omega(\mathbf{r}) \, d\mathbf{r} = 0 \quad (6)$$

Posons maintenant, pour tout couple $(u, v) \in H_0^1(\Omega) \times H_0^1(\Omega)$,

$$a(u, v) = D \int_{\Omega} \nabla u \cdot \nabla v \, d\mathbf{r} + \int_{\Omega} K(\mathbf{r}) uv \, d\mathbf{r} \quad \text{et} \quad m(u, v) = \int_{\Omega} uv \, d\mathbf{r}$$

La forme bilinéaire $a(\cdot, \cdot)$ est continue et coercive sur $H_0^1(\Omega)$, tandis que $m(\cdot, \cdot)$ est le produit scalaire de $L^2(\Omega)$.

La formulation variationnelle du problème s'écrit alors :

$$\frac{d}{dt}m(\rho, \omega) + a(\rho, \omega) = 0, \quad \forall \omega \in H_0^1(\Omega) \quad (7)$$

Ce qui nous permet de définir la notion de solution faible, comme suit :

**Définition 4 (Solution faible)** On dit que $\rho$ est une solution faible du problème sur $[0, T]$ associée à la donnée initiale $\rho_0 \in L^2(\Omega)$ si $\rho$ vérifie :

— $\rho \in L^2(0, T; H_0^1(\Omega))$ et $\partial_t \rho \in L^2(0, T; H^{-1}(\Omega))$,

— $\rho(0, \mathbf{r}) = \rho_0(\mathbf{r})$ pour tout $r \in \Omega$ (au sens des traces dans $L^2$),

— l'égalité variationnelle (7) pour presque tout $t \in (0, T)$ et pour tout $v \in H_0^1(\Omega)$.

## 2.2 Caractère bien posé du problème

Soit $A : H_0^1(\Omega) \to H^{-1}(\Omega)$ défini par $Au = -D\Delta u + Ku$ où $K \in L^\infty(\Omega)$, $K \ge 0$. L'opérateur $A$ est auto-adjoint et positif sur $L^2(\Omega)$.

Le théorème de Hille-Yosida (BREZIS, 2011, Th. 7.4) s'applique alors directement :

**Théorème 1 (Corollaire du théorème de Hille-Yosida)** Pour toute donnée initiale $\rho_0 \in L^2(\Omega)$ il existe une unique fonction

$$\rho \in C([0, T]; L^2(\Omega)) \cap L^2(0, T; H_0^1(\Omega)), \quad \partial_t \rho \in L^2(0, T; H^{-1}(\Omega)), \quad (8)$$

qui vérifie l'équation 7 et la condition initiale $\rho(0) = \rho_0$. La solution dépend continûment de $\rho_0$ dans $L^2(\Omega)$. En particulier, si deux données initiales coïncident alors les solutions coïncident.

Le problème est par conséquent bien posé.
### 3 Modélisation des pièges – Choix des termes puits

L'objectif de cette section est d'étudier différents termes puits afin d'obtenir une modélisation réaliste. Le terme puit correspond à la diminution locale de la densité de moustique S telle que :

$$S : \mathbb{R}^2 \to \mathbb{R}_+ \\ \mathbf{r} \mapsto S(\mathbf{r})$$

Ce terme vient compléter le taux de mortalité naturel $M_0$ et apparaît dans l'équation de diffusion modifiée.

#### 3.1 Terme puits ponctuel

Le premier modèle que nous considérons est celui d'un terme puits ponctuel, représentant un piège dont l'effet est confiné en un point unique du plan, noté $\mathbf{a} \in \mathbb{R}^2$. Dans ce cadre, la diminution locale de la densité de moustiques est modélisée par une distribution de Dirac centrée en $\mathbf{a}$, ce qui conduit à :

$$S(\mathbf{r}) = S_0 \delta(\mathbf{r} - \mathbf{a}),$$

<figure><img src="image_5.png" /><figcaption>FIGURE 3 – Visualisation 2D du terme puits ponctuel $S(\mathbf{r}) = S_0\delta(\mathbf{r})$.</figcaption></figure>

où $S_0 > 0$ désigne l'intensité du piège. Ce modèle correspond à un piège idéalement localisé, d'efficacité maximale au point $\mathbf{a}$, sans influence alentour. En réinjectant ce terme dans l'équation de diffusion, on obtient :

$$\frac{\partial \rho}{\partial t}(\mathbf{r}, t) = D \Delta \rho(\mathbf{r}, t) - M_0 \rho(\mathbf{r}, t) - S_0 \delta(\mathbf{r} - \mathbf{a}) \rho(\mathbf{r}, t).$$

Ce dernier terme étant singulier, le produit $\delta \cdot \rho$ n'est bien défini que si $\rho$ est continue en $\mathbf{a}$, auquel cas on pose $\delta(\mathbf{r} - \mathbf{a})\rho(\mathbf{r}, t) = \rho(\mathbf{a}, t)\delta(\mathbf{r} - \mathbf{a})$. L'équation s'écrit alors, au sens faible :

$$\frac{\partial \rho}{\partial t} = D \Delta \rho - M_0 \rho - S_0 \rho(\mathbf{a}, t) \delta(\mathbf{r} - \mathbf{a}).$$

Un tel terme pose des difficultés analytiques et numériques en raison de sa nature distributionnelle. En effet, la multiplication d'une fonction potentiellement non régulière $\rho$ par une distribution comme $\delta(\mathbf{r} - \mathbf{a})$ n'est pas nécessairement bien définie. Toutefois, ces difficultés peuvent être contournées dans le cadre des méthodes de résolution numérique reposant sur une formulation variationnelle, telles que la méthode des éléments finis (MEF).
## 3.2 Terme puits constant

Un second modèle, plus régulier, consiste à considérer un terme puits constant sur une région $A \subset \mathbb{R}^2$, représentant un piège réparti spatialement avec une efficacité uniforme. On suppose que la capture s'effectue à un taux constant $S_0 > 0$ sur l'ensemble $A$, et que l'effet est nul ailleurs. On définit alors :

$$S(\mathbf{r}) = S_0 \mathbf{1}_A(\mathbf{r}),$$

où $\mathbf{1}_A$ désigne la fonction indicatrice de $A$. Ce modèle correspond à un piège étendu, tel qu'un filet, une grille ou une zone de piégeage fixe, dont la portée est connue et mesurable. La visualisation de ce puits est donnée dans la figure 4 où l'on a pris un disque pour $A$.

<figure><img src="image_6.png" /><figcaption>FIGURE 4 – Visualisation 2D du puits constant $S(\mathbf{r}) = S_0 \mathbf{1}_A(\mathbf{r})$.</figcaption></figure>

En remplaçant ce terme dans l'équation de diffusion, on obtient :

$$\frac{\partial \rho}{\partial t}(\mathbf{r}, t) = D \Delta \rho(\mathbf{r}, t) - M_0 \rho(\mathbf{r}, t) - S_0 \mathbf{1}_A(\mathbf{r}) \rho(\mathbf{r}, t).$$

Cette expression est bien définie dans le cadre classique, car $\mathbf{1}_A$ est une fonction mesurable bornée. La discontinuité introduite à la frontière de $A$ peut néanmoins engendrer des effets de bord dans les solutions numériques, mais ces effets sont maîtrisables, notamment via un raffinement local du maillage dans les méthodes de type éléments finis ou volumes finis.

L'avantage principal de ce modèle réside dans sa simplicité d'implémentation et sa compatibilité directe avec les méthodes de résolution variationnelles. Ce type de terme s'intègre parfaitement dans la formulation faible, et est donc parfaitement adapté à une implémentation efficace, même sur des domaines complexes, et permet de modéliser fidèlement des pièges répartis sur une surface donnée.

## 3.3 Terme puits gaussien

Un troisième modèle consiste à considérer un terme puits diffus, c'est-à-dire une zone de capture dont l'intensité décroît progressivement avec la distance au centre du piège. Un choix naturel est de modéliser cette décroissance par une fonction gaussienne centrée en un point $\mathbf{a} \in \mathbb{R}^2$, ce qui conduit à :

$$S(\mathbf{r}) = S_0 \exp \left( - \frac{\| \mathbf{r} - \mathbf{a} \|^2}{2\sigma^2} \right),$$

où $S_0 > 0$ est l'intensité maximale du piège et $\sigma > 0$ contrôle l'étendue spatiale de son influence. Ce modèle correspond à un piège dont l'action est plus forte au centre, mais non nulle à distance finie : il est donc particulièrement adapté à la modélisation de dispositifs à efficacité décroissante.

En introduisant ce terme dans l'équation de diffusion, on obtient :

$$\frac{\partial \rho}{\partial t}(\mathbf{r}, t) = D \Delta \rho(\mathbf{r}, t) - M_0 \rho(\mathbf{r}, t) - S_0 e^{-\frac{\|\mathbf{r}-\mathbf{a}\|^2}{2\sigma^2}} \rho(\mathbf{r}, t).$$
<figure><img src="image_7.png" /><figcaption>FIGURE 5 – Visualisation 2D du terme puits gaussien $S(\mathbf{r}) = S_0 e^{-\|\mathbf{r}\|^2/(2\sigma^2)}$.</figcaption></figure>

Cette expression est lisse sur $\mathbb{R}^2$, ce qui présente un double avantage : d'une part, l'équation reste définie de manière classique et, d'autre part, cette régularité facilite les calculs numériques, notamment dans les méthodes variationnelles.

Ce modèle gaussien peut aussi être interprété comme une régularisation naturelle du terme puits ponctuel, dans la limite $\sigma \to 0$. Il combine ainsi réalisme physique et confort d'analyse, tout en permettant un contrôle fin de l'intensité et de la portée du piège à l'aide des paramètres $S_0$ et $\sigma$.

Pour des raisons de réalisme physique et de correspondance aux données empiriques préalables, nous considérons dans cet article le terme puits gaussien.

## 3.4 Modélisation multi-pièges

Dans le cas où plusieurs pièges sont répartis dans le domaine, une approche naturelle consiste à modéliser chaque piège individuellement par un terme puits $S_j(\mathbf{r})$ de type gaussien, centré en un point $\mathbf{a}_j \in \mathbb{R}^2$. Le terme puits total est alors défini comme la somme de ces contributions locales :

$$S(\mathbf{r}) = \sum_{j=1}^{N} S_j(\mathbf{r}), \quad \text{avec} \quad S_j(\mathbf{r}) = S_0 \exp \left( - \frac{\|\mathbf{r} - \mathbf{a}_j\|^2}{2\sigma^2} \right),$$

où N est le nombre total de pièges, $S_0$ l'intensité maximale (éventuellement variable selon j), et $\sigma$ le paramètre caractéristique de dispersion de chaque piège, que l'on suppose identique pour tous les pièges.

Cette construction permet de représenter un ensemble de pièges comme une superposition d'effets localisés, tout en conservant une fonction $S(\mathbf{r})$ lisse et bien adaptée à une résolution numérique par éléments finis ou volumes finis. En pratique, cela signifie que la configuration des pièges est entièrement déterminée par la donnée des positions $(\mathbf{a}_1, \dots, \mathbf{a}_N)$.

## 4 Méthode des différences finies

L'objectif de cette section est de détailler la résolution numérique de l'équation 3 lorsqu'on applique des schémas de différences finies. Pour appliquer la méthode des différences finies, on approxime les dérivées de $\rho$ via des schémas discrets, temporels et spatiaux. Pour plus de détails, voir STRIKWERDA, 2004.
## 4.1 Définition des maillages

Pour approximer $\frac{\partial \rho}{\partial t}$ on définit le maillage temporel suivant :

$$t^n = n \Delta t, \quad \text{où} \quad n \in \{0, \dots, N_t - 1\} \qquad (9)$$

Puis on définit un maillage pour approximer les dérivées spatiales :

$$x_i = ih \quad \text{et} \quad y_j = jh \quad \text{pour} \quad i \in 0, \dots, N_x - 1 \quad \text{et} \quad j \in 0, \dots, N_y - 1 \qquad (10)$$

On note alors :

$$\rho_{i,j}^n \approx \rho(x_i, y_j, t^n) \qquad (11)$$

## 4.2 Définition du schéma de récurrence

On utilise le schéma d'Euler explicite en temps et le schéma à 5 points en espace :

$$\frac{\rho_{i,j}^{n+1} - \rho_{i,j}^n}{\Delta t} \approx \frac{\partial \rho}{\partial t} \qquad (12)$$

$$\frac{\rho_{i+1,j}^n + \rho_{i-1,j}^n + \rho_{i,j+1}^n + \rho_{i,j-1}^n - 4\rho_{i,j}^n}{h^2} \approx \Delta\rho \qquad (13)$$

En réarrangeant, on obtient une relation de récurrence qui nous permet de définir $\rho_{i,j}^{n+1}$ en fonction des valeurs de $\rho$ au temps $n$ :

$$\rho_{i,j}^{n+1} = \rho_{i,j}^n + \Delta t \left[ D \frac{\rho_{i+1,j}^n + \rho_{i-1,j}^n + \rho_{i,j+1}^n + \rho_{i,j-1}^n - 4 \rho_{i,j}^n}{h^2} - K(x_i, y_j) \rho_{i,j}^n \right] \qquad (14)$$

L'implémentation d'un tel schéma est relativement simple et très facilement adaptable pour les différentes formes prises par les pièges. Toutefois, il est nécessaire de s'assurer de la stabilité du schéma utilisé via la condition CFL (RICHARD COURANT et LEWY, 1956).

**Théorème 2 (Condition CFL)** Soit $K_{\max} = \sup_{\mathbf{r} \in \Omega} K(\mathbf{r})$. Le schéma est convergent d'ordre $O(h^2+\Delta t)$ si

$$\Delta t \le \frac{h^2}{4D + h^2 K_{\max}} \qquad (15)$$

**Preuve 1 (Condition CFL)** Voir ALLAIRE, 2012 & HUNDSDORFER et VERWER, 2003.

## 4.3 Complexité et précision

Le calcul du laplacien et du terme source en chaque point de la grille constitue une opération de complexité $O(N^2)$, en posant $N := N_x = N_y$ (on suppose une grille carrée pour simplifier l'analyse). Pour un nombre total d'itérations temporelles noté $n$, la complexité globale en temps de calcul s'élève donc à $(n \times N^2)$.

Concernant la précision du schéma, l'approximation spatiale du laplacien repose sur le schéma aux différences finies à 5 points, qui est une méthode consistante et convergente d'ordre 2 : l'erreur de discrétisation spatiale est en $O(h^2)$, avec $h$ le pas d'espace. En ce qui concerne la discrétisation temporelle, un schéma d'Euler explicite a été retenu, lequel est d'ordre 1 en temps : l'erreur temporelle est donc en $O(\Delta t)$.

## 4.4 Algorithme et résultats

L'algorithme repose sur la relation de récurrence définie par l'équation 4.2. Le code source est disponible sur le dépôt GitLab du projet.

Les figure 6 et figure 7 présentent la densité surfacique de moustiques $\rho(\mathbf{r}, t)$ simulée par la méthode
<figure><img src="image_8.png" /><figcaption>FIGURE 6 – Simulation à l'aide de la méthode des différences finies en $t_1 = n_1 \Delta t$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 7 – Simulation à l'aide de la méthode des différences finies en $t_2 = n_2 \Delta t$</figcaption></figure>

des différences finies à des instants $t_1 = n_1 \Delta t$ et $t_2 = n_2 \Delta t$ ($t_1 < t_2$). On y observe une propagation progressive du front de diffusion vers les zones périphériques, ainsi que des zones de faible densité localisées précisément là où sont placés les pièges gaussiens. Ce comportement est en excellent accord avec la théorie de la diffusion enrichie de termes puits. Le gradient de couleur reflète la valeur de $\rho$ en moustiques · m$^{-2}$ et doit être interprété comme une simple visualisation de la solution numérique.

Toutefois, les paramètres ($D, S_0, \sigma, \Delta t, h$) employés ici ne sont pas seulement des estimations : ils sont purement *ad hoc*, dépourvus de toute signification physique concrète et ne prétendent pas refléter un contexte réel.
# 5 Méthode des éléments finis

L'objectif de cette section est d'approcher numériquement une solution de l'équation 3 via la méthode des éléments finis. Pour les fondements de la méthode et les éléments d'analyse fonctionnelle nécessaires à son application, voir PIRONNEAU, 1988.

L'intérêt de la méthode réside dans la robustesse de sa formulation variationnelle et dans la possibilité d'obtenir une précision locale accrue : les zones de forte variation (autour des pièges, notamment) y sont mieux résolues, et la convergence est assurée théoriquement.

## 5.1 Définition des éléments finis $\mathbb{P}_1$

Pour l'approximation par la méthode des éléments finis, on définit l'espace discret

$$V_h = \mathrm{Vect}\{\varphi_{i,j} \mid (ih, jh) \in N_h\} \qquad (16)$$

où la fonction chapeau associée est donnée par

$$\varphi_{i,j}(x, y) = \max\left(0, 1 - \frac{|x - ih|}{h} - \frac{|y - jh|}{h}\right) \qquad (17)$$

Cette fonction, qui vaut 1 en $(ih, jh)$ et décroît linéairement à zéro en dehors du losange

$$\{(x, y) \in \mathbb{R}^2 \mid |x - ih| + |y - jh| \leq h\} \qquad (18)$$

est définie de manière linéaire par morceaux sur chaque élément $T \in \mathcal{T}_h$.

Pour approximer la solution $\rho$, on la représente par une combinaison linéaire de fonctions de base locales :

$$\rho_h(x, y, t) = \sum_{(i,j) \in I} U_{i,j}(t) \varphi_{i,j}(x, y) \qquad (19)$$

où $U_{i,j}(t)$ sont des coefficients temporels à déterminer. En réutilisant la formulation variationnelle définie en section 2, il convient de trouver $\rho_h(t) \in V_h$ tel que, pour tout $\omega_h \in V_h$,

$$\frac{d}{dt} \int_{\Omega} \rho_h \omega_h dxdy + \int_{\Omega} D \nabla \rho_h \cdot \nabla \omega_h dxdy + \int_{\Omega} K(x, y) \rho_h \omega_h dxdy = 0 \qquad (20)$$

## 5.2 Définition du maillage

On considère $h > 0$ le pas de discrétisation et $\Omega = [-N, N] \times [-M, M]$ le domaine de diffusion. On définit la grille de nœuds par :

$$N_h = \{(ih, jh) \mid i = -N, \dots, N, j = -M, \dots, M\} \qquad (21)$$

où $(x_i, y_j) = (ih, jh)$ représente un point du maillage.

Pour $i = -N, \dots, N - 1$ et $j = -M, \dots, M - 1$, on associe à chaque nœud le carré

$$S_{i,j} = [ih, (i+1)h] \times [jh, (j+1)h] \qquad (22)$$

qui constitue une subdivision élémentaire de $\Omega$.

Chaque carré $S_{i,j}$ est découpé en deux triangles (voir figure 8), définis par

$$T_{i,j}^{(1)} = \mathrm{conv}\{(ih, jh), ((i+1)h, jh), ((i+1)h, (j+1)h)\} \qquad (23)$$

$$T_{i,j}^{(2)} = \mathrm{conv}\{(ih, jh), ((i+1)h, (j+1)h), (ih, (j+1)h)\} \qquad (24)$$

où

$$\mathrm{conv} : \mathcal{P}(\mathbb{R}^d) \to \mathcal{P}(\mathbb{R}^d) \\ A \mapsto \mathrm{conv}(A)$$
<figure><img src="image_10.png" /><figcaption>FIGURE 8 – Subdivision d'un carré $S_{i,j}$ en deux triangles $T_{i,j}^{(1)}$ et $T_{i,j}^{(2)}$</figcaption></figure>

désigne l'application qui associe à tout ensemble fini $A = \{x_1, \dots, x_m\} \subset \mathbb{R}^d$ son enveloppe convexe :

$$ \mathrm{conv}(A) = \{ \lambda_1 x_1 + \dots + \lambda_m x_m \mid \lambda_i \ge 0, \sum_{i=1}^{m} \lambda_i = 1 \} \quad (25) $$

Autrement dit, $\mathrm{conv}(A)$ est le plus petit ensemble convexe contenant $A$.

Le maillage triangulaire <fnref n="1" /> complet est alors

$$ \mathcal{T}_h = \bigcup_{i=0}^{N-1} \bigcup_{j=0}^{M-1} \{T_{i,j}^{(1)}, T_{i,j}^{(2)}\} \quad (26) $$

## 5.3 Étude de l'EDO

Après projection sur la base $(\varphi_{i,j})_{i,j}$, l'assemblage conduit au système d'équations différentielles ordinaires

$$ M \dot{U}(t) + (A + R) U(t) = 0 \quad (27) $$

où les matrices sont définies par :

$$ M_{(i,j),(k,l)} = \int_{\Omega} \varphi_{i,j}(x,y) \varphi_{k,l}(x,y) \, dxdy \quad (28) $$

$$ A_{(i,j),(k,l)} = \int_{\Omega} D \nabla \varphi_{i,j}(x,y) \cdot \nabla \varphi_{k,l}(x,y) \, dxdy \quad (29) $$

$$ R_{(i,j),(k,l)} = \int_{\Omega} K(x,y) \varphi_{i,j}(x,y) \varphi_{k,l}(x,y) \, dxdy \quad (30) $$

Pour étudier l'EDO matricielle obtenue, on peut utiliser le schéma de Crank-Nicolson qui consiste à approximer la dérivée temporelle par une différence centrée et à considérer une moyenne des termes à $t^n$ et $t^{n+1}$ pour le second membre.

Le schéma s'écrit donc :

$$ M \frac{U^{n+1} - U^n}{\Delta t} + \frac{1}{2}(A + R)(U^{n+1} + U^n) = 0 \quad (31) $$

En réarrangeant, on obtient la relation de récurrence suivante :

$$ \left( M + \frac{\Delta t}{2}(A + R) \right) U^{n+1} = \left( M - \frac{\Delta t}{2}(A + R) \right) U^n \quad (32) $$

<footnote n="1">¹Le choix d'un tel maillage se justifie par sa facilité d'implémentation et son adaptation au problème.</footnote>
## 5.4 Complexité et précision

Dans le cadre de la méthode des éléments finis, la complexité dépend principalement de l'assemblage du système linéaire (matrice de raideur et vecteur second membre), ainsi que de sa résolution à chaque pas de temps. En considérant un maillage régulier composé de $N^2$ éléments triangulaires (avec $N$ le nombre de divisions par direction dans un domaine carré), la phase d'assemblage est de complexité $O(N^2)$, tandis que la résolution d'un système linéaire creux de taille $N^2 \times N^2$ est dans notre cas en $O(N^3)$. La complexité temporelle globale s'élève donc à $O(n \times N^3)$.

Sur le plan de la précision, la méthode des éléments finis linéaires ($\mathbb{P}_1$) sur un maillage régulier fournit une approximation spatiale du second ordre : l'erreur est en $O(h^2)$, où $h$ désigne la taille caractéristique des éléments du maillage. Pour le temps, l'utilisation de la méthode de Crank-Nicolson entraîne une erreur en $O(\Delta t^2)$. La précision globale dépend donc du raffinement spatial et du pas de temps choisi, avec une convergence garantie par la théorie de Céa (voir CIARLET, 1978) sous hypothèses classiques de régularité.

## 5.5 Algorithme et résultats

L'algorithme mis en œuvre repose sur la relation de récurrence formulée par l'équation 5.3. L'implémentation en langage Python, issue directement de cette modélisation mathématique, est accessible sur le dépôt GitLab associé au projet. Relativement à la méthode des différences finies, cette implémentation se révèle plus complexe en raison de la nature de la méthode.

<figure><img src="image_11.png" /><figcaption>FIGURE 9 – Simulation à l'aide de la méthode des éléments finis en $t'_1 = n'_1 \Delta t$</figcaption></figure>

Les figure 9 et figure 10 présentent la densité surfacique de moustiques $\rho(\mathbf{r}, t)$ simulée par la méthode des éléments finis aux instants $t'_1 = n'_1 \Delta t$ et $t'_2 = n'_2 \Delta t$ (avec $t'_1 < t'_2$). On observe des résultats très similaires à ceux obtenus dans la section 4, ce qui confirme la pertinence de cette approche alternative.

Tout comme dans la section 4, les paramètres ($D, S_0, \sigma, \Delta t, h$) employés ici ne sont pas seulement des estimations : ils sont purement *ad hoc*, dépourvus de toute signification physique concrète et ne prétendent pas refléter un contexte réel.
<figure><img src="image_12.png" /><figcaption>FIGURE 10 – Simulation à l'aide de la méthode des différences finies en $t'_2 = n'_2 \Delta t$</figcaption></figure>

# 6 Méthode des volumes finis

Cette section introduit la méthode des volumes finis : partant des lois de conservation intégrales, elle définit le volume de contrôle, discrétise les flux sur le maillage puis présente les principaux schémas afin de résoudre l'équation 4.

## 6.1 Définition du maillage

Un maillage $\mathcal{T}$ de $\Omega$ est défini à partir de cellules de contrôle (ou volumes de contrôle), notées $C$, selon les règles suivantes :

1. Chaque cellule $C$ est un polygone convexe, et l'ensemble des cellules constitue une partition de $\Omega$, c'est-à-dire :

$$ \bar{\Omega} = \bigcup_{C \in \mathcal{T}} \bar{C} \quad \text{et} \quad \mathring{C} \cap \mathring{L} = \emptyset \quad \text{pour tout } C \neq L \qquad (33) $$

2. À chaque cellule $C$, on associe un point intérieur $\mathbf{r}_C$, appelé centre, qui satisfait les conditions suivantes :

- Si $L$ est une cellule voisine de $C$, alors les centres $\mathbf{r}_C$ et $\mathbf{r}_L$ sont distincts et le segment les reliant est orthogonal à l'arête commune $e$ entre les deux cellules $C$ et $L$. Cette arête est notée $e = (C|L)$.

- Si une arête $e$ du bord $\partial\Omega$ appartient à la cellule $C$, alors la droite passant par le centre $\mathbf{r}_C$ et perpendiculaire à $e$ coupe nécessairement cette arête.

Un tel maillage est dit admissible dans le cadre de la méthode des Volumes Finis.

## 6.2 Discrétisation en espace

Sur une cellule $C$, on note :

$$ f_C(t) = \frac{1}{|C|} \int_C K(\mathbf{r}) \rho(\mathbf{r}, t) \, dx \qquad (34) $$

et $\mathbf{n}_C$ désigne la normale unitaire dirigée vers l'extérieur de $C$. On pose

$$ \rho_C(t) = \frac{1}{|C|} \int_C \rho(\mathbf{r}, t) \, dx \qquad (35) $$

et on a, pour $t \in [0, T]$,
<figure><img src="image_13.png" /><figcaption>FIGURE 11 – Cellules intérieures admissibles</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 12 – Cellules sur les bords</figcaption></figure>

$$ |C| \frac{d\rho_C}{dt}(t) - \nu \sum_{e \in \mathcal{E}_C} \int_e \nabla \rho \cdot \mathbf{n}_{C,e} d\Gamma = |C| f_C(t) \quad (36) $$

où $\mathcal{E}_C$ désigne l'ensemble des arêtes de la cellule $C$ et $\mathbf{n}_{C,e}$ est la normale unitaire à $e$ dirigée vers l'extérieur de $C$.

On approche le flux

$$ - \int_e \nabla \rho \cdot \mathbf{n}_{C,e} d\Gamma \approx F_{C,e} \quad (37) $$

avec

$$ F_{C,e} = \frac{(\rho_C - \rho_L)}{|\mathbf{r}_C - \mathbf{r}_L|} |e| \quad \text{si } e \notin \partial\Omega \quad (38) $$

Pour simplifier les calculs, on suppose que les cellules sont toutes des cellules intérieures, on peut justifier cette approximation en prenant une surface d'étude assez grande pour qu'il n'y ait aucun moustique sur les cellules non intérieures. Le schéma Volumes Finis de discrétisation en espace s'écrit ainsi, pour $t \in [0, T]$,

$$ |C| \frac{d\rho_C}{dt}(t) + \sum_{\substack{e \in \mathcal{E}_C \\ e \notin \partial\Omega}} \nu F_{C,e}(t) = |C| f_C(t) \quad (39) $$

## 6.3 Discrétisation en temps

Soit $n \in \mathbb{N}^*$ et $\Delta t = \frac{T}{n} > 0$ le pas de temps. On pose $t^k = k\Delta t$ pour $k = \{0, 1, \dots, n\}$. L'approximation de la moyenne de $\rho$ sur une cellule $C$ au temps $t^k$ est notée :

$$ \rho_C^k \approx \frac{1}{|C|} \int_C \rho(\mathbf{r}, t^k) \, dx \quad (40) $$

Le schéma explicite s'écrit, pour toute cellule $C \in \mathcal{T}$, et $k = \{0, \dots, n\}$,

$$ |C| \left( \frac{\rho_C^{k+1} - \rho_C^k}{\Delta t} \right) + \sum_{\substack{e \in \mathcal{E}_C \\ e \notin \partial\Omega}} \nu F_{C,e}^k = |C| f_C^k \quad (41) $$

où

$$ F_{C,e}^k = \frac{\rho_C^k - \rho_L^k}{|\mathbf{r}_C - \mathbf{r}_L|} |e|, \quad \text{pour } e = (C|L), \ e \notin \partial\Omega \quad (42) $$

On applique ces équations à chaque cellule $C$ (tous supposées comme intérieures) :

$$ \rho_C^{k+1} = \rho_C^k - \frac{\nu \Delta t}{|C|} \sum_{e \in \mathcal{E}_C, e=(C|L)} \frac{|e|}{|\mathbf{r}_C - \mathbf{r}_L|} (\rho_C^k - \rho_L^k) + \Delta t f_C^k \quad (43) $$
## Système linéaire

Soit $N$ le nombre total de cellules de contrôle intérieures. On regroupe les inconnues sous forme de vecteur :

$$ \rho^k = (\rho_1^k, \rho_2^k, \dots, \rho_N^k)^\top \qquad (44) $$

et le vecteur des sources :

$$ \mathbf{f}^k = (f_1^k, f_2^k, \dots, f_N^k)^\top \qquad (45) $$

où

Le système linéaire correspondant devient alors :

$$ \rho^{k+1} = (I_N - \nu \Delta t H^{-1} A) \rho^k + \Delta t \mathbf{f}^k \qquad (46) $$

où $I_N$ est la matrice identité d'ordre $N$, $A$ est la matrice du Laplacien discrétisé, et $H$ est la matrice diagonale contenant les volumes des cellules :

$$ H = \begin{pmatrix} |C_1| & 0 & \cdots & 0 \\ 0 & |C_2| & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & |C_N| \end{pmatrix} \qquad (47) $$

$$ A = \begin{pmatrix} \sum_{\substack{e_{C,L} \in \partial C \\ e_{C,L} \not\subset \partial \Omega}} \frac{|e_{C,L}|}{d_{C,L}} & \cdots & -\frac{|e_{C,L_1}|}{d_{C,L_1}} & \cdots & -\frac{|e_{C,L_2}|}{d_{C,L_2}} & \cdots \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots \\ -\frac{|e_{L_1,C}|}{d_{L_1,C}} & \cdots & \sum_{\substack{e_{L_1,L} \in \partial L_1 \\ e_{L_1,L} \not\subset \partial \Omega}} \frac{|e_{L_1,L}|}{d_{L_1,L}} & \cdots & -\frac{|e_{L_1,L_2}|}{d_{L_1,L_2}} & \cdots \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots \\ -\frac{|e_{L_2,C}|}{d_{L_2,C}} & \cdots & -\frac{|e_{L_2,L_1}|}{d_{L_2,L_1}} & \cdots & \sum_{\substack{e_{L_2,L} \in \partial L_2 \\ e_{L_2,L} \not\subset \partial \Omega}} \frac{|e_{L_2,L}|}{d_{L_2,L}} & \cdots \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots \end{pmatrix} \qquad (48) $$

On a noté $e_{C,L} = (C|L) \not\subset \partial \Omega$ l'arête commune aux cellules $C$ et $L$, et $d_{C,L} = |\mathbf{r}_C - \mathbf{r}_L|$. Les cellules $L_1, L_2, \dots$ sont les cellules intérieures ayant une arête en commun avec la cellule $C$.

## 6.4 Algorithme et résultats

La première étape consiste à définir le domaine de calcul et à générer un maillage adapté. Le code utilise une triangulation de Delaunay, qui correspond au maillage admissible dans le cadre de cette méthode. Le maillage est créé à partir de points aléatoires et un exemple est donné figure 13.

Une fois le maillage généré, le code calcule plusieurs propriétés géométriques essentielles :

- Pour chaque paire de triangles voisins, la longueur de leur arête commune est calculée. Cette information est cruciale pour le calcul des flux entre les volumes de contrôle.

- Les centres des triangles sont calculés, ainsi que les distances entre centres de triangles voisins. Ces distances sont utilisées pour approximer les gradients dans le schéma numérique.

- L'air des triangles.

- Les bornes d'intégration sur un triangle en 2D.

Finalement, la fonction principale calcule les différents termes de l'équation 4, et retourne la solution globale sous la forme d'une matrice où chaque ligne correspond à un instant de temps et chaque colonne à la valeur de $\rho$ sur un triangle.
<figure><img src="image_15.png" /><figcaption>FIGURE 13 – Visualisation d'un maillage par la méthode de triangulation de Delaunay</figcaption></figure>

Bien que la démarche précédente – de la triangulation admissible jusqu'au système linéaire (équation 46-équation 48) – établisse rigoureusement la consistance, la stabilité et un ordre de convergence du schéma, sa mise en œuvre numérique n'a pu être menée à terme dans le cadre de ce projet. À chaque pas de temps, l'algorithme reconstruit les longueurs d'arêtes, les centres de gravité des triangles et la matrice des volumes de contrôle ; la complexité de ces opérations croît en $O(N^3)$ avec le nombre $N$ de nœuds lorsqu'on raffine le maillage. Dans ces conditions, la validation expérimentale de la MVF a été suspendue : Les comparaisons quantitatives présentées plus loin se limitent donc aux schémas MDF et MEF, tandis que l'étude exhaustive de la MVF est remise à de futurs travaux disposant de ressources de calcul plus importantes.

# 7 Comparaison des algorithmes

L'objectif de cette section est de comparer les différentes méthodes. Comme indiqué dans la sous-section 6.4, la méthode des volumes finis, dont les résultats se sont révélés non concluants, ne sera pas approfondie. Seules la méthode des différences finies et la méthode des éléments finis sont comparées. La figure 14 repose sur le script Python accessible à l'adresse suivante : comparaison.py, hébergé dans le dépôt GitLab.

La figure 14 met en évidence que, dans le cadre expérimental considéré, la méthode des différences finies présente des performances temporelles inférieures à celles de la méthode des éléments finis. Cette différence peut s'expliquer par plusieurs facteurs. Tout d'abord, la complexité algorithmique de la méthode des différences finies est en $O(n \times N^2)$, contre $O(n \times N^3)$ pour la méthode des éléments finis, où $n$ désigne le nombre d'itérations temporelles et $N$ la taille caractéristique du maillage spatial. Cette différence d'ordre de grandeur dans la complexité explique en grande partie l'écart observé. Dans les simulations menées, la valeur $N = 100$ a été utilisée, ce qui rend cet effet particulièrement significatif. Par ailleurs, la constante implicite associée à la notation en $\mathcal{O}$ pourrait être sensiblement plus faible pour la méthode des différences finies. Cela se traduit par un avantage pratique pour un nombre modéré d'itérations, en dépit d'une complexité asymptotiquement plus favorable pour la méthode des éléments finis. Enfin, il convient de noter que la méthode des éléments finis induit un surcoût d'initialisation non négligeable, lié notamment au calcul des matrices d'assemblage et à la construction des structures topologiques associées. Dans les configurations étudiées, ce surcoût initial est de l'ordre de 0,5 seconde, même pour des cas peu itératifs. Toutefois, cet investissement initial tend à être compensé lorsque la
taille du problème augmente ou lorsqu'un raffinement spatial important est requis.

<figure><img src="image_16.png" /><figcaption>FIGURE 14 – Comparaison des performances temporelles des algorithmes</figcaption></figure>
# 8 Conclusion

La diffusion spatiale de la population de moustiques dans un milieu bidimensionnel a été modélisée par une EDP de type équation de la chaleur. L'effet de pièges à moustiques a été intégré par l'introduction de termes puits localisés (points de divergence négative), avec une préférence accordée à un modèle gaussien pour ces puits. Plusieurs approches numériques ont été comparées pour résoudre cette équation enrichie, notamment la méthode des différences finies (MDF), la méthode des éléments finis (MEF) et la méthode des volumes finis (MVF). La robustesse supérieure de la MEF a été mise en évidence, tandis que la simplicité de mise en œuvre de la MDF a été confirmée. En revanche, des limites ont été observées pour la MVF. Des fonctions Python associées ont été développées et implémentées pour la résolution de l'équation et le calcul du nombre de moustiques sur une surface donnée. Les performances de ces méthodes ont fait l'objet d'une comparaison en termes de précision et de temps de calcul. Enfin, la présente étude a été menée dans le contexte de la technique de l'insecte stérile (SIT) et de la lutte anti-vectorielle, et servira à une étude statistique plus large servant à estimer les paramètres de diffusion.
Références

ALLAIRE, G. (2012). *Analyse numérique et optimisation : une introduction avec applications* (2e éd.)
[Voir §2.2.3, équation (1.31) pour la condition CFL du schéma explicite de diffusion]. École
polytechnique.

BREZIS, H. (2011). *Functional Analysis, Sobolev Spaces and Partial Differential Equations*. Springer.

CAMARERO, R. (2024). *Triangulation de Delaunay* (rapp. tech.) [Présentation pour le cours MEC6212 :
Génération de maillages]. Polytechnique Montréal, Département de génie mécanique.

CIARLET, P. G. (1978). *The Finite Element Method for Elliptic Problems*. North-Holland.

EVANS, L. C. (2010). *Partial Differential Equations* (2e éd., T. 19). American Mathematical Society.

HUNDSDORFER, W., & VERWER, J. G. (2003). *Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations* (T. 33) [Chapitre 3 : la présence d'un terme réaction linéaire $-K\rho$ n'altère pas la borne CFL car il est strictement dissipatif]. Springer.

LIONS, J.-L., & MAGENES, E. (1968). *Problèmes aux limites non homogènes et applications. Vol. 1*. Dunod.

MUÑOZ, F. (2023). *SIT Methods – Techniques stérilisantes pour la lutte anti-vectorielle* [En cours de rédaction, non publié]. https://umr-astre.pages.mia.inra.fr/sit-methods/

NURHAYATI, S., YUNIANTO, B., RAMADHANI, T., IKAWATI, B., SANTOSO, B., & RAHAYU, A. (2013).
Controlling Aedes aegypti population as DHF vector with radiation based-sterile insect tech-
nique in Banjarnegara Regency, Central Java. *Jurnal Sains dan Teknologi Nuklir Indonesia*,
14 (1), 1-10. Récupérée 9 mai 2025, à partir de https://core.ac.uk/outputs/231115268

PIRONNEAU, O. (1988). *Finite Element Methods for Fluids* [Notes de cours données à l'Université Pierre et Marie Curie & INRIA; à paraître en 1988, Wiley]. Wiley. Récupérée 9 mai 2025, à partir de https://www.ljll.fr/pironneau/publi/publications/OPfemInFluids.pdf

RICHARD COURANT, K. F., & LEWY, H. (1956). *On the Partial Difference Equations of Mathematical Physics* (AEC Research and Development Report NYO-7689) [Translated by Phyllis Fox]. New York University. https://archive.org/details/onpartialdiffere00cour

SCHEID, J.-F. (2017). *Volumes Finis* (rapp. tech.) [Master 2 IMOL, Méthodes numériques avancées pour la résolution des EDP]. Lorraine, France, Université de Lorraine.

STRIKWERDA, J. C. (2004). *Finite Difference Schemes and Partial Differential Equations* (rapp. tech.)
[Livre électronique, deuxième édition, Classics in Applied Mathematics]. SIAM — Society for
Industrial and Applied Mathematics. https://doi.org/10.1137/1.9780898717938