- Email: [email protected]
Une méthode particulaire déterministe pour des équations diffusives non linéaires
C. R. Acad. Sci. Paris, t. 332, Série I, p. 369–376, 2001 Analyse numérique/Numerical Analysis (Équations aux dérivées partielles/Partial Differential Equations)
Une méthode particulaire déterministe pour des équations diffusives non linéaires Pierre-Louis LIONS a , Sylvie MAS-GALLIC b,c a
Ceremade, UMR CNRS 7534, Université Paris-9–Dauphine, place du Maréchal-de-Lattre-de-Tassigny, 75775 Paris cedex 16, France b CMAP, École polytechnique, 91128 Palaiseau cedex, France c Département de mathématiques, Université Evry-Val-d’Essonne, 91025 Evry cedex, France (Reçu le 27 novembre 2000, accepté le 4 décembre 2000)
Résumé.
Nous étudions dans cette Note une méthode particulaire déterministe pour des équations de type chaleur (ou Fokker–Planck) ou de type milieux poreux. Cette méthode est basée sur une approximation de ces équations par des équations de transport non linéaires et nous montrons la convergence de cette approximation. Enfin, nous présentons quelques expériences numériques pour l’équation de la chaleur et pour un exemple d’équations de milieux poreux. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
A deterministic particle method for nonlinear diffusion equations Abstract.
We study in this Note a deterministic particle method for heat (or Fokker–Planck) equations or for porous media equations. This method is based upon an approximation of these equations by nonlinear transport equations and we prove the convergence of that approximation. Finally, we present some numerical experiments for the heat equation and for an example of porous media equations. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Abridged English version We study in this Note a numerical approximation method for heat type equations (possibly nonlinear). We present this approach on a porous media equation ∂u 1 − ∆u2 = 0 ∂t 2
(1)
∂u − ∆u = 0. ∂t
(2)
and on the usual heat equation
Note présentée par Pierre-Louis L IONS. S0764-4442(00)01795-X/FLA 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés.
369
P.-L. Lions, S. Mas-Gallic
We detail here the case of (1) and refer to the text below (in French) for the case of (2). In order to simplify the presentation and notation, we only consider the case of nonnegative, periodic (of period Ti in each xi , 1 i N ) solutions of (1). The idea of the approximation is to write 1 ∆u2 = div(V u) with V = ∇u, 2 to consider V as a “transport velocity” and finally to smooth V for exampleby convolution, i.e., to consider Vε = V ∗ ζε , where ζε = ε1N ζ ε· , ε > 0, ζ ∈ S(RN ) (for instance. . .) and RN ζ dx = 1. This leads to ∂u = div(Vε u), ∂t
Vε = ∇u ∗ ζε .
(3)
If we prescribe an initial condition u|t=0 = u0 , where u0 is a nonnegative bounded measure, we explain below that there exists a unique solution of (3) which is smooth in (x, t) when u0 is smooth. In addition, (3) is a conservative (nonlinear) transport equation which can be solved by particle methods and we refer the reader to the French version below for a discussion of the origin of this approach. We present below numerical experiments in order to compare the numerically obtained solutions with the solutions of (1) and we show below a rigorous proof of the convergence of solutions uεof (3) to solutions of (1) as ε goes to 0. More precisely, if uε solves (3) with uε |t=0 = uε0 ∈L1 such that uε0 |Log uε0 | < ∞ and uε0 converges weakly in L1 to u0 and if ζε = ρˇε ∗ ρε with ρε = ε1N ρ ε· , ρ ∈ S(RN ), ρ > 0, RN ρ dx = 1, ρˇ = ρ(−·), then uε converges weakly in L1 to the solution of (1) satisfying: u|t=0 = u0 , u ∈ L2t (H1) ∩ Ct (L1 ), supt0 u |Log u| < ∞. In addition, if uε0 converges strongly to u0 and for uε0 Log uε0 → u0 Log u0 , then uε converges to u in L1 , uniformly in t ∈ [0, T ] (for all T ∈ (0, ∞)).
ε
1. Introduction Nous étudions dans cette Note une méthode d’approximation numérique d’équations de type chaleur (non linéaires ou non). Nous présentons cette méthode sur deux exemples à savoir une équation de type milieux poreux ∂u 1 − ∆u2 = 0 ∂t 2
(1)
et l’équation de la chaleur ∂u − ∆u = 0. (2) ∂t Pour simplifier la présentation et les notations, nous ne considérons ici que le cas d’équations posées avec des conditions périodiques (u est périodique de période Ti > 0 en xi pour 1 i N ) et de solutions positives de (1) ou (2). L’idée est d’écrire 1 ∆u2 = div(V u) avec V = ∇u, 2
∆u = div(V u) avec V =
∇u = ∇ Log u, u
de considérer V comme une « vitesse de transport » et enfin de régulariser V par convolution, i.e. de remplacer V par Vε = V ∗ ζε , où ε > 0, ζε = ε1N ζ ε· , ζ ∈ S(RN ) (par exemple. . .) et RN ζ dx = 1. Ceci conduit donc à l’approximation suivante de (1), respectivement (2) ∂u = div(Vε u), ∂t
370
Vε = ∇u ∗ ζε
(3)
Une méthode particulaire déterministe pour des équations diffusives non linéaires
respectivement ∂u = div(Vε u), Vε = (∇ Log u) ∗ ζε . ∂t On peut également réécrire (4) en introduisant v = Log u et on obtient ∂v = Vε · ∇v + div Vε , ∂t
Vε = (∇v) ∗ ζε .
(4)
(5)
ε Remarquons enfin que, dans le cas de (2), on peut également régulariser V en introduisant Vε = ∇u∗ζ u∗ζε ce qui donne à l’évidence un modèle très proche de (4) (avec essentiellement la même précision que (4) par rapport à ε . . .). Nous montrons dans ce qui suit que les solutions de (3) et (4) (ou (5)) existent et convergent, quand ε tend vers 0, vers les solutions de (1) et (2) respectivement. D’autre part, les équations de transport (3) ou (4) peuvent être résolues numériquement par des méthodes particulaires classiques et nous présentons quelques expériences numériques illustrant notre approche. Signalons pour conclure cette introduction que l’approche introduite ici est très naturelle du point de vue physique et a été mise en évidence au début des années 80 par des physiciens des plasmas (voir [3] et [4] par exemple) pour la résolution d’équations cinétiques collisionnelles telle l’équation de Fokker– Planck. Réécrivant l’équation sous forme d’une équation de transport dont la vitesse, déterministe, prend en compte les effets de collision et de diffusion, ils ont ouvert la voie aux méthodes numériques telles que les méthodes particulaires. Depuis son introduction, cette méthode a été régulièrement utilisée, pour des calculs sur les semiconducteurs ([1] et [9], par exemple) mais le plus souvent pour des applications proches de la mécanique des fluides ([5], [6], [13], [12], [14], . . .). Pourtant, excepté quelques résultats partiels récents ([7] et [8]), aucune analyse mathématique n’en avait encore été fournie. Comme il a déjà été vu, le point de départ de la méthode consiste à remplacer l’équation initiale (1) ou (2) par une équation de transport non linéaire (3) ou (4), dans laquelle la vitesse de diffusion est régularisée. La méthode numérique choisie ensuite ici pour résoudre les équations de transport est une méthode particulaire mais il est tout à fait possible d’utiliser d’autres méthodes numériques.
2. Équation de type milieux poreux Commençons par la résolution de (3) en imposant une condition initiale : u|t=0 = u0 , où u0 0 est une mesure bornée. La classe de solutions de (3) considérée correspond à une mesure ( 0) u telle que u(t) soit pour tout t 0 une mesure bornée et u(t) est continue pour la convergence faible des mesures et telle que (3) ait lieu au sens des distributions. P ROPOSITION 1. – Il existe une unique solution u de (3) qui est régulière en (x, t) si u0 est régulière. Esquisse de la démonstration. – Comme on a d’après (3) la conservation, pour tout t 0, de u(t) dx k (et que u 0), on déduit aisément une estimation sur Vε dans L∞ t (Cx ) (pour tout k 0) et l’existence d’une solution s’obtient en considérant pour deux solutions u et v la sans difficulté. L’unicité se démontre quantité M (t) = sup (u − v)ϕ | ϕ ∈ W1,∞ , ∇ϕ ∞ 1 . ✷ La convergence des solutions uε de (3) vers la solution de (1) est assurée par le résultat suivant pour lequel on suppose que ζε = ρˇε ∗ ρε avec ρε = ε1N ρ ε· , ρˇε = ρε (−·), ρ ∈ S(RN ), ρ 0, ρ dx = 1. T HÉORÈME 2. – Soit uε0 ∈ L1 , uε0 0 telle que supε uε0 |Log uε0 | < ∞ ; on note par uε la solution faiblement dans L1 vers la solution de (1) vérifiant : de (3) telle que uε |t=0 ≡ uε0 . Alors, uε converge 2 1 1 plus, la convergence de uε vers u est forte u|t=0 = u0 , u ∈ Lt (H ) ∩ Ct (L ), supt0 u |Log u| < ∞. De 1 ε 1 ε dans Ct (L ) si u0 converge fortement dans L vers u0 et u0 Log uε0 dx converge vers u0 Log u0 dx.
371
P.-L. Lions, S. Mas-Gallic
Remarque 1. – Le résultat qui précède et la démonstration s’adaptent facilement au cas de solutions de signe quelconque : il faut alors remplacer u2 par |u|u et Vε par (∇|u|) ∗ ζε . Remarque 2. – La démonstration du théorème 2 montre que u2 ∗ ρε converge vers u dans L2 (par exemple. . .). Démonstration du théorème 2. – On observe tout d’abord que l’on a, pour tout t 0, 2 d ε ε u Log u dx + ∇(uε ∗ ρε ) dx = 0, dt
(6)
1 d’où l’on déduit des bornes sur uε ∗ ρε dans L2t (H1 ) et sur uε (1 + |Log uε |) dans L∞ t (L ). On remarque ensuite que, pour toute fonction ϕ régulière, on a ε ∇(uε ∗ ζε )u · ∇ϕ dx = ∇(uε ∗ ρε )(uε · ∇ϕ) ∗ ρε dx 1 ε 2 = ∇(u ∗ ρε ) · ∇ϕ dx + ∇(uε ∗ ρε ) (uε · ∇ϕ) ∗ ρε − ∇ϕ(uε ∗ ρε ) . 2 De plus, on vérifie aisément que l’on a (uε · ∇ϕ) ∗ ρε − ∇ϕ(uε ∗ ρε ) C D2 ϕ ∞ ε(uε ∗ ρε ). Ceci permet alors de montrer d’une part que (uε ∗ ρε ) est relativement compacte dans L2 et que l’équation (1) a lieu en faisant tendre ε vers 0. Enfin, la convergence forte se déduit aisément de l’identité (6).
3. Équation de la chaleur Tout d’abord, nous signalons que les estimations a priori données ci-dessous impliquent aisément l’existence et l’unicité d’une solution v de (5) telle que v ∈ Ct (L1 ), ev ∈ Ct (L1 ) et v|t=0 ≡ v0 , où v0 ∈ L1 et ev0 ∈ L1 . De plus, cette solution est régulière dès que v0 l’est. Nous considérons alors v0ε et supposons ε que v0ε converge dans L1 vers v0 et que uε0 = ev0 converge dans L1 vers u0 . Enfin, nous notons par v ε la solution de (5) telle que vε |t=0 ≡ v0ε et par uε = evε . Nous aurons besoin de l’hypothèse suivante sur ζε : ζε = ρˇε ∗ ρε , ρε = ε1N ρ ε· , ρˇε = ρε (−·), ρ ∈ S(RN ), ρ dx = 1, ρ 0 et ∇ ρ = 0 sur { ρ = 0} où ϕ désigne la transformée de Fourier de ϕ. T HÉORÈME 3. – Sous les hypothèses précédentes, v ε et uε convergent dans Ct (L1 ), quand ε tend vers 0, vers, respectivement v et u où u = ev , u est la solution de (2) telle que u|t=0 ≡ u0 et v est la solution de ∂v = |∇v|2 + ∆v, ∂t
v|t=0 = v0 ,
v ∈ Ct L1 ∩ L2t H1 .
(7)
De plus, v ε ∗ ρε converge vers v dans L2t (H1 ). Esquisse de la démonstration du théorème 3. – Dans une première étape, nous établissons des bornes a priori (indépendantes de ε). Comme uε > 0 résout (3), nous avons bien sûr ε ε u (t) dx = u0 dx −→ u0 dx. (8) ε
D’autre part, on déduit de (5) d dt
−v dx + ε
|∇v ε ∗ ρε |2 dx = 0.
(9)
1 ε 2 1 On déduit alors aisément de (8)–(9) des bornes sur uε et v ε dans L∞ t (L ) et sur v ∗ ρε dans Lt (H )) ε + ∞ q (remarquons également que (v ) est bornée dans Lt (L ) pour tout q < ∞). Quitte à extraire des
372
Une méthode particulaire déterministe pour des équations diffusives non linéaires
sous-suites, on peut donc supposer que uε et v ε convergent au sens des mesures vers u faiblement et v qui vérifient d’après (8) et les bornes précédentes : u dx = u0 pour tout t 0, u ev et 1 2 1 v ∈ L∞ t (L ) ∩ Lt (H ). La deuxième étape consiste à montrer que v vérifie l’inéquation suivante ∂v |∇v|2 + ∆v, ∂t
v|t=0 = v0 .
(10)
Pour ce faire, on remarque tout d’abord que les bornes précédentes sur (v ε )+ et v ε ∗ ρε impliquent que v ε ∗ ϕε est bornée dans L2t (Lq ) pour 2 q N2N −2 si N = 3, 2 q < ∞ si N = 2, q = +∞ si N = 1, où · 1 ∞ ϕε = εN ϕ ε , ϕ ∈ L et il existe α > N tel que supz (|z|α |ϕ(z)|) < ∞. À l’aide de ces bornes et des hypothèses sur ρε , on peut établir le fait suivant (∇v ε ) ∗ ζε · ∇v − |∇v ε ∗ ρε |2 −→ 0 ε
dans L1x,t .
L’inégalité (10) s’en déduit alors facilement. Enfin, la conclusion de la démonstration s’obtient en multipliant (10) par ev /(1 + δev ) pour δ > 0, d’où
ev ev δe2v 1 d v 2 2 Log 1 + δe dx = |∇v| − ∇v ∇ dx 0. dx = |∇v| dt δ 1 + δev 1 + δev (1 + δev )2 En faisant tendre δ vers 0, on en déduit pour tout t > 0 u0 dx = u(t) dx ev(t) dx ev0 dx = u0 dx. Ceci implique donc que u ≡ ev et on peut alors aisément conclure. ✷ 4. Expériences numériques Nous avons choisi ici d’utiliser comme méthode numérique, une méthode particulaire. Ce choix a été motivé par le fait que les mesures de Dirac dont les centres évoluent à la vitesse de diffusion, sont solutions exactes au sens des distributions de l’équation de transport non linéaire obtenue. En principe, il apparaît dans l’analyse de convergence de ce genre de méthode (voir [10] par exemple) que l’erreur se décompose en deux parties, dont la première, due à l’approximation du modèle de départ équation (1) (resp. (2)) par l’équation (3) (resp. (4)) est de type asymptotique et la seconde, plus spécifiquement numérique, est due à l’approximation de la condition initiale par une combinaison linéaire de mesures de Dirac. La convergence de l’erreur asymptotique, c’est-à-dire la convergence de la solution de l’équation (3) (resp. (4)) vers la solution de équation (1) (resp. (2)) a été démontrée dans les paragraphes précédents mais l’ordre de convergence en ε n’est pas clair et la forte nonlinéarité du problème rend cette analyse particulièrement difficile. En ce qui concerne l’erreur numérique, elle fait intervenir généralement le nombre, N , de particules et le paramètre ε, la fonction ζε jouant le rôle de fonction de forme de particules. Sans préciser d’avantage les hypothèses, disons que si la fonction ζ et la condition initiale u0 sont assez régulières, on obtient une erreur du type
1 u − uh = O u − uh = O 1 , + O εr , ε m m N (εN )
où u est la solution de l’équation (1) (resp. (2)), uh = p∈N αp δ(x − Xp (t)) son approximation par des
mesures de Dirac, uhε = p∈N αp ζε (x − Xp (t)) sa version régularisée et les positions Xp (t) sont les
373
P.-L. Lions, S. Mas-Gallic
positions des particules. Les normes sont des normes d’espaces de Sobolev habituels et les exposants m et r dépendent, en particulier des fonctions ζ et u0 . La seconde estimation fait apparaître, pour un nombre N de particules donné, l’existence d’une valeur optimale de ε, qui peut être recherchée lors de tests numériques. Les illustrations numériques monodimensionnelles qui suivent confirment, d’une part dans le cas de l’équation des milieux poreux et d’autre part dans le cas de l’équation de la chaleur, les résultats de l’analyse théorique. Les calculs ont été fait avec un nombre fixé de 250 particules, 100 pas de temps et une fonction cut-off P2 par morceaux de classe C1 ; les erreurs sont calculées dans les normes pour lesquelles la convergence a été obtenue. 4.1. Équation de type milieu poreux. – La fonction u0 étant donnée, on résout numériquement l’équation ∂u = div(Vε u), ∂t
Vε = ∇u ∗ ζε
où Vε , vitesse de diffusion définie précédemment a été réécrite en tenant compte des propriétés du produit de convolution. Ceci nous permet de l’appliquer à l’approximation particulaire uh de u et ainsi de définir la vitesse des particules. La position des particules au temps t est donnée par la résolution du système d’équations différentielles dXp (t) = − αq ∇ζε Xp (t) − Xq (t) , dt
Xp (0) = x0p ,
q∈N
auquel un schéma de Runge–Kutta d’ordre deux a été appliqué. Les courbes de la figure 1 montrent que, pour ε = 0.20, la solution autosimilaire dite de Barenblatt est bien calculée ; on constate néanmoins que les points d’intersection de la parabole avec l’axe n’ont pas tout à fait la bonne vitesse sans doute parce que les particules initialement placées en ces deux points n’y sont plus. Pour résoudre cette difficulté, plusieurs solutions sont envisageables. La courbe d’erreur en fonction du paramètre ε, figure 2, montre d’une part sur la courbe de gauche, l’existence d’une valeur optimale de ε recherchée, elle est ici obtenue pour ε = 0.10 et donne d’autre part, sur la courbe de droite, une estimation de l’ordre de convergence en ε qui ici à peu près 1.6. Précisons enfin que le pas de temps ∆t est, comme cela été établi dans [2], lié à la valeur de ε et non pas au nombre N de particules et qu’une étude heuristique nous à fait choisir ici ∆t = ε. 4.2. Équation de la chaleur. – Soit u0 une fonction donnée. Rappelons tout d’abord l’équation à résoudre respectivement ∂u = div Vε u , ∂t
∇u ∗ ζε Vε = , u ∗ ζε
Figure 1. – À gauche, condition initiale ; à droite, solutions approchée (trait plein) et exacte (trait pointillé). Figure 1. – Left, initial condition; right, approximate (plain line) and exact (doted line) solutions.
374
(11)
Une méthode particulaire déterministe pour des équations diffusives non linéaires
Figure 2. – À gauche, erreur en norme L1 en x, L2 en t en fonction du paramètre ε ; à droite, zoom et droite de plus proche pente donnant l’ordre de convergence. Figure 2. – Left, L1 in x, L2 in t error with respect to ε; right, closer view plus nearest line with order of convergence.
Figure 3. – À gauche, condition initiale ; à droite, solutions approchée (trait plein) et exacte (trait pointillé). Figure 3. – Left, initial condition; right, approximate (plain line) and exact (doted line) solutions.
Figure 4. – Erreur en norme L1 ln L1 en x, L2 en t en fonction du paramètre ε. Figure 4. – L1 in x, L2 in t error with respect to ε.
où Vε est la vitesse de diffusion. Cette écriture nous permet de définir la vitesse des particules dont la position au temps t est donnée par la résolution du système d’équations différentielles
dXp q∈N αq ∇ζε (Xp (t) − Xq (t)) (t) = − , Xp (0) = x0p , dt q∈N αq ζε (Xp (t) − Xq (t)) qui est résolu par un schéma de Runge–Kutta d’ordre deux.
375
P.-L. Lions, S. Mas-Gallic
Les courbes, figure 3, montrent que la Gaussienne, solution de l’équation de la chaleur, est bien calculée, la courbe d’erreur en fonction de ε, figure 4, montrant une convergence moins rapide que dans le cas de l’équation précédente. Précisons également, que le pas de temps, ∆t, qui est calculé en fonction de la valeur de ε a été fixé, après une étude heuristique, à ∆t = ε2 /2. Remarque 3. – On peut avoir des problèmes de calcul si la solution est proche de zéro. Pour pallier cette difficulté, un seuil numérique η > 0 est introduit et le système couplé d’équations différentielles ordinaires résolu est en fait le suivant dXp q∈N αq ∇ζε (Xp (t) − Xq (t)) q∈N αq ζε (Xp (t) − Xq (t)) (t) = − , Xp (0) = x0p . 2 2 dt αq ζε (Xp (t) − Xq (t)) + η q∈N
Conclusion. – La méthode particulaire qui a été utilisée ici, est, comme le sont les méthodes de marche aléatoire, une méthode particulaire à poids fixes ce qui la rend beaucoup moins coûteuse en temps de calcul que, par exemple, la méthode à poids variable (communément appelée P.S.E., particle strength exchange). Elle est clairement très dépendante de la fonction ζ utilisée et, au cours des tests trois fonctions ζ différentes ont été utilisées, le noyau dit de Fejer, une fraction rationnelle classique en méthode de vortex et la fonction qui nous a servi dans les résultats présentés ici. Le choix de cette dernière fonction est un choix lié à la rapidité et à la précision des calculs. Les démonstrations détaillées et des résultats numériques plus complets feront l’objet d’un prochain article. Références bibliographiques [1] Degond P., Mustieles F.J., Approximation of diffusion equations by deterministic convections of particles, SIAM J. Sci. Statis. Comput. 11 (1990) 293–310. [2] Degond P., Mas-Gallic S., The weighted particle method for convection–diffusion equations, part I: the case of an isotropic viscosity, part II: the anisotropic case, Math. Comput. 53 (1989) 485–526. [3] Fronteau J., Combis X., Hadronic J. 7 (1984) 911–940. [4] Grmela M., Fronteau J., Tellez-Arenas A., Hadronic J. 3 (1980) 1203–1226. [5] Huberson S., Rivoalen E., Hauville F., Simulation numérique des équations de Navier–Stokes 3d par une méthode particulaire, C. R. Acad. Sci. Paris, Série IIb 324 (1997) 543–549. [6] Huberson S., Rivoalen E., Numerical simulation of axisymmetric viscous flows by means of particle method, J. Comput. Phys. 152 (1) (1999) 1–31. [7] Lacombe G., Analyse d’une équation de vitesse de diffusion, C. R. Acad. Sci. Paris, Série I 329 (5) (1999) 383– 386. [8] Lacombe G, Mas-Gallic S., Presentation and analysis of a diffusion velocity method, in: ESAIM Proceedings, Vortex Flow and Related Numerical Methods III, 1998. [9] Mas-Gallic S., Une présentation de la méthode de vitesse diffusion, Conférence en l’honneur de Marc Feix, Orléans, in: Dynamical Systems, Plasmas and Gravitation, Leach P.G.L., Bouquet S.E., Rouet J.L., Fijalkow E. (Eds.), Springer-Verlag, 1997, pp. 74–81. [10] Mas-Gallic S., Raviart P.-A., A particle method for first order symmetric systems, Numer. Math. 51 (1987) 323– 352. [11] Mustieles F.J., Thèse de doctorat, CMAP, École polytechnique, Palaiseau, 1990. [12] Ogami Y., Akamatsu T., Viscous flow simulation using a discrete vortex method, J. Computers and Fluids 19 (1991) 433–441. [13] Rivoalen E., Thèse de doctorat, Le Havre, 1994. [14] Strickland J.H., Kempka S.N., Wolfe W.P., Viscous diffusion of vorticity using the diffusion velocity concept, in: ESAIM Proceedings, Vortex Flow and Related Numerical Methods II, Vol. 1, 1996, pp. 135–151.
376
Une méthode particulaire déterministe pour des équations diffusives non linéaires Pierre-Louis LIONS a , Sylvie MAS-GALLIC b,c a
Ceremade, UMR CNRS 7534, Université Paris-9–Dauphine, place du Maréchal-de-Lattre-de-Tassigny, 75775 Paris cedex 16, France b CMAP, École polytechnique, 91128 Palaiseau cedex, France c Département de mathématiques, Université Evry-Val-d’Essonne, 91025 Evry cedex, France (Reçu le 27 novembre 2000, accepté le 4 décembre 2000)
Résumé.
Nous étudions dans cette Note une méthode particulaire déterministe pour des équations de type chaleur (ou Fokker–Planck) ou de type milieux poreux. Cette méthode est basée sur une approximation de ces équations par des équations de transport non linéaires et nous montrons la convergence de cette approximation. Enfin, nous présentons quelques expériences numériques pour l’équation de la chaleur et pour un exemple d’équations de milieux poreux. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
A deterministic particle method for nonlinear diffusion equations Abstract.
We study in this Note a deterministic particle method for heat (or Fokker–Planck) equations or for porous media equations. This method is based upon an approximation of these equations by nonlinear transport equations and we prove the convergence of that approximation. Finally, we present some numerical experiments for the heat equation and for an example of porous media equations. 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Abridged English version We study in this Note a numerical approximation method for heat type equations (possibly nonlinear). We present this approach on a porous media equation ∂u 1 − ∆u2 = 0 ∂t 2
(1)
∂u − ∆u = 0. ∂t
(2)
and on the usual heat equation
Note présentée par Pierre-Louis L IONS. S0764-4442(00)01795-X/FLA 2001 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés.
369
P.-L. Lions, S. Mas-Gallic
We detail here the case of (1) and refer to the text below (in French) for the case of (2). In order to simplify the presentation and notation, we only consider the case of nonnegative, periodic (of period Ti in each xi , 1 i N ) solutions of (1). The idea of the approximation is to write 1 ∆u2 = div(V u) with V = ∇u, 2 to consider V as a “transport velocity” and finally to smooth V for exampleby convolution, i.e., to consider Vε = V ∗ ζε , where ζε = ε1N ζ ε· , ε > 0, ζ ∈ S(RN ) (for instance. . .) and RN ζ dx = 1. This leads to ∂u = div(Vε u), ∂t
Vε = ∇u ∗ ζε .
(3)
If we prescribe an initial condition u|t=0 = u0 , where u0 is a nonnegative bounded measure, we explain below that there exists a unique solution of (3) which is smooth in (x, t) when u0 is smooth. In addition, (3) is a conservative (nonlinear) transport equation which can be solved by particle methods and we refer the reader to the French version below for a discussion of the origin of this approach. We present below numerical experiments in order to compare the numerically obtained solutions with the solutions of (1) and we show below a rigorous proof of the convergence of solutions uεof (3) to solutions of (1) as ε goes to 0. More precisely, if uε solves (3) with uε |t=0 = uε0 ∈L1 such that uε0 |Log uε0 | < ∞ and uε0 converges weakly in L1 to u0 and if ζε = ρˇε ∗ ρε with ρε = ε1N ρ ε· , ρ ∈ S(RN ), ρ > 0, RN ρ dx = 1, ρˇ = ρ(−·), then uε converges weakly in L1 to the solution of (1) satisfying: u|t=0 = u0 , u ∈ L2t (H1) ∩ Ct (L1 ), supt0 u |Log u| < ∞. In addition, if uε0 converges strongly to u0 and for uε0 Log uε0 → u0 Log u0 , then uε converges to u in L1 , uniformly in t ∈ [0, T ] (for all T ∈ (0, ∞)).
ε
1. Introduction Nous étudions dans cette Note une méthode d’approximation numérique d’équations de type chaleur (non linéaires ou non). Nous présentons cette méthode sur deux exemples à savoir une équation de type milieux poreux ∂u 1 − ∆u2 = 0 ∂t 2
(1)
et l’équation de la chaleur ∂u − ∆u = 0. (2) ∂t Pour simplifier la présentation et les notations, nous ne considérons ici que le cas d’équations posées avec des conditions périodiques (u est périodique de période Ti > 0 en xi pour 1 i N ) et de solutions positives de (1) ou (2). L’idée est d’écrire 1 ∆u2 = div(V u) avec V = ∇u, 2
∆u = div(V u) avec V =
∇u = ∇ Log u, u
de considérer V comme une « vitesse de transport » et enfin de régulariser V par convolution, i.e. de remplacer V par Vε = V ∗ ζε , où ε > 0, ζε = ε1N ζ ε· , ζ ∈ S(RN ) (par exemple. . .) et RN ζ dx = 1. Ceci conduit donc à l’approximation suivante de (1), respectivement (2) ∂u = div(Vε u), ∂t
370
Vε = ∇u ∗ ζε
(3)
Une méthode particulaire déterministe pour des équations diffusives non linéaires
respectivement ∂u = div(Vε u), Vε = (∇ Log u) ∗ ζε . ∂t On peut également réécrire (4) en introduisant v = Log u et on obtient ∂v = Vε · ∇v + div Vε , ∂t
Vε = (∇v) ∗ ζε .
(4)
(5)
ε Remarquons enfin que, dans le cas de (2), on peut également régulariser V en introduisant Vε = ∇u∗ζ u∗ζε ce qui donne à l’évidence un modèle très proche de (4) (avec essentiellement la même précision que (4) par rapport à ε . . .). Nous montrons dans ce qui suit que les solutions de (3) et (4) (ou (5)) existent et convergent, quand ε tend vers 0, vers les solutions de (1) et (2) respectivement. D’autre part, les équations de transport (3) ou (4) peuvent être résolues numériquement par des méthodes particulaires classiques et nous présentons quelques expériences numériques illustrant notre approche. Signalons pour conclure cette introduction que l’approche introduite ici est très naturelle du point de vue physique et a été mise en évidence au début des années 80 par des physiciens des plasmas (voir [3] et [4] par exemple) pour la résolution d’équations cinétiques collisionnelles telle l’équation de Fokker– Planck. Réécrivant l’équation sous forme d’une équation de transport dont la vitesse, déterministe, prend en compte les effets de collision et de diffusion, ils ont ouvert la voie aux méthodes numériques telles que les méthodes particulaires. Depuis son introduction, cette méthode a été régulièrement utilisée, pour des calculs sur les semiconducteurs ([1] et [9], par exemple) mais le plus souvent pour des applications proches de la mécanique des fluides ([5], [6], [13], [12], [14], . . .). Pourtant, excepté quelques résultats partiels récents ([7] et [8]), aucune analyse mathématique n’en avait encore été fournie. Comme il a déjà été vu, le point de départ de la méthode consiste à remplacer l’équation initiale (1) ou (2) par une équation de transport non linéaire (3) ou (4), dans laquelle la vitesse de diffusion est régularisée. La méthode numérique choisie ensuite ici pour résoudre les équations de transport est une méthode particulaire mais il est tout à fait possible d’utiliser d’autres méthodes numériques.
2. Équation de type milieux poreux Commençons par la résolution de (3) en imposant une condition initiale : u|t=0 = u0 , où u0 0 est une mesure bornée. La classe de solutions de (3) considérée correspond à une mesure ( 0) u telle que u(t) soit pour tout t 0 une mesure bornée et u(t) est continue pour la convergence faible des mesures et telle que (3) ait lieu au sens des distributions. P ROPOSITION 1. – Il existe une unique solution u de (3) qui est régulière en (x, t) si u0 est régulière. Esquisse de la démonstration. – Comme on a d’après (3) la conservation, pour tout t 0, de u(t) dx k (et que u 0), on déduit aisément une estimation sur Vε dans L∞ t (Cx ) (pour tout k 0) et l’existence d’une solution s’obtient en considérant pour deux solutions u et v la sans difficulté. L’unicité se démontre quantité M (t) = sup (u − v)ϕ | ϕ ∈ W1,∞ , ∇ϕ ∞ 1 . ✷ La convergence des solutions uε de (3) vers la solution de (1) est assurée par le résultat suivant pour lequel on suppose que ζε = ρˇε ∗ ρε avec ρε = ε1N ρ ε· , ρˇε = ρε (−·), ρ ∈ S(RN ), ρ 0, ρ dx = 1. T HÉORÈME 2. – Soit uε0 ∈ L1 , uε0 0 telle que supε uε0 |Log uε0 | < ∞ ; on note par uε la solution faiblement dans L1 vers la solution de (1) vérifiant : de (3) telle que uε |t=0 ≡ uε0 . Alors, uε converge 2 1 1 plus, la convergence de uε vers u est forte u|t=0 = u0 , u ∈ Lt (H ) ∩ Ct (L ), supt0 u |Log u| < ∞. De 1 ε 1 ε dans Ct (L ) si u0 converge fortement dans L vers u0 et u0 Log uε0 dx converge vers u0 Log u0 dx.
371
P.-L. Lions, S. Mas-Gallic
Remarque 1. – Le résultat qui précède et la démonstration s’adaptent facilement au cas de solutions de signe quelconque : il faut alors remplacer u2 par |u|u et Vε par (∇|u|) ∗ ζε . Remarque 2. – La démonstration du théorème 2 montre que u2 ∗ ρε converge vers u dans L2 (par exemple. . .). Démonstration du théorème 2. – On observe tout d’abord que l’on a, pour tout t 0, 2 d ε ε u Log u dx + ∇(uε ∗ ρε ) dx = 0, dt
(6)
1 d’où l’on déduit des bornes sur uε ∗ ρε dans L2t (H1 ) et sur uε (1 + |Log uε |) dans L∞ t (L ). On remarque ensuite que, pour toute fonction ϕ régulière, on a ε ∇(uε ∗ ζε )u · ∇ϕ dx = ∇(uε ∗ ρε )(uε · ∇ϕ) ∗ ρε dx 1 ε 2 = ∇(u ∗ ρε ) · ∇ϕ dx + ∇(uε ∗ ρε ) (uε · ∇ϕ) ∗ ρε − ∇ϕ(uε ∗ ρε ) . 2 De plus, on vérifie aisément que l’on a (uε · ∇ϕ) ∗ ρε − ∇ϕ(uε ∗ ρε ) C D2 ϕ ∞ ε(uε ∗ ρε ). Ceci permet alors de montrer d’une part que (uε ∗ ρε ) est relativement compacte dans L2 et que l’équation (1) a lieu en faisant tendre ε vers 0. Enfin, la convergence forte se déduit aisément de l’identité (6).
3. Équation de la chaleur Tout d’abord, nous signalons que les estimations a priori données ci-dessous impliquent aisément l’existence et l’unicité d’une solution v de (5) telle que v ∈ Ct (L1 ), ev ∈ Ct (L1 ) et v|t=0 ≡ v0 , où v0 ∈ L1 et ev0 ∈ L1 . De plus, cette solution est régulière dès que v0 l’est. Nous considérons alors v0ε et supposons ε que v0ε converge dans L1 vers v0 et que uε0 = ev0 converge dans L1 vers u0 . Enfin, nous notons par v ε la solution de (5) telle que vε |t=0 ≡ v0ε et par uε = evε . Nous aurons besoin de l’hypothèse suivante sur ζε : ζε = ρˇε ∗ ρε , ρε = ε1N ρ ε· , ρˇε = ρε (−·), ρ ∈ S(RN ), ρ dx = 1, ρ 0 et ∇ ρ = 0 sur { ρ = 0} où ϕ désigne la transformée de Fourier de ϕ. T HÉORÈME 3. – Sous les hypothèses précédentes, v ε et uε convergent dans Ct (L1 ), quand ε tend vers 0, vers, respectivement v et u où u = ev , u est la solution de (2) telle que u|t=0 ≡ u0 et v est la solution de ∂v = |∇v|2 + ∆v, ∂t
v|t=0 = v0 ,
v ∈ Ct L1 ∩ L2t H1 .
(7)
De plus, v ε ∗ ρε converge vers v dans L2t (H1 ). Esquisse de la démonstration du théorème 3. – Dans une première étape, nous établissons des bornes a priori (indépendantes de ε). Comme uε > 0 résout (3), nous avons bien sûr ε ε u (t) dx = u0 dx −→ u0 dx. (8) ε
D’autre part, on déduit de (5) d dt
−v dx + ε
|∇v ε ∗ ρε |2 dx = 0.
(9)
1 ε 2 1 On déduit alors aisément de (8)–(9) des bornes sur uε et v ε dans L∞ t (L ) et sur v ∗ ρε dans Lt (H )) ε + ∞ q (remarquons également que (v ) est bornée dans Lt (L ) pour tout q < ∞). Quitte à extraire des
372
Une méthode particulaire déterministe pour des équations diffusives non linéaires
sous-suites, on peut donc supposer que uε et v ε convergent au sens des mesures vers u faiblement et v qui vérifient d’après (8) et les bornes précédentes : u dx = u0 pour tout t 0, u ev et 1 2 1 v ∈ L∞ t (L ) ∩ Lt (H ). La deuxième étape consiste à montrer que v vérifie l’inéquation suivante ∂v |∇v|2 + ∆v, ∂t
v|t=0 = v0 .
(10)
Pour ce faire, on remarque tout d’abord que les bornes précédentes sur (v ε )+ et v ε ∗ ρε impliquent que v ε ∗ ϕε est bornée dans L2t (Lq ) pour 2 q N2N −2 si N = 3, 2 q < ∞ si N = 2, q = +∞ si N = 1, où · 1 ∞ ϕε = εN ϕ ε , ϕ ∈ L et il existe α > N tel que supz (|z|α |ϕ(z)|) < ∞. À l’aide de ces bornes et des hypothèses sur ρε , on peut établir le fait suivant (∇v ε ) ∗ ζε · ∇v − |∇v ε ∗ ρε |2 −→ 0 ε
dans L1x,t .
L’inégalité (10) s’en déduit alors facilement. Enfin, la conclusion de la démonstration s’obtient en multipliant (10) par ev /(1 + δev ) pour δ > 0, d’où
ev ev δe2v 1 d v 2 2 Log 1 + δe dx = |∇v| − ∇v ∇ dx 0. dx = |∇v| dt δ 1 + δev 1 + δev (1 + δev )2 En faisant tendre δ vers 0, on en déduit pour tout t > 0 u0 dx = u(t) dx ev(t) dx ev0 dx = u0 dx. Ceci implique donc que u ≡ ev et on peut alors aisément conclure. ✷ 4. Expériences numériques Nous avons choisi ici d’utiliser comme méthode numérique, une méthode particulaire. Ce choix a été motivé par le fait que les mesures de Dirac dont les centres évoluent à la vitesse de diffusion, sont solutions exactes au sens des distributions de l’équation de transport non linéaire obtenue. En principe, il apparaît dans l’analyse de convergence de ce genre de méthode (voir [10] par exemple) que l’erreur se décompose en deux parties, dont la première, due à l’approximation du modèle de départ équation (1) (resp. (2)) par l’équation (3) (resp. (4)) est de type asymptotique et la seconde, plus spécifiquement numérique, est due à l’approximation de la condition initiale par une combinaison linéaire de mesures de Dirac. La convergence de l’erreur asymptotique, c’est-à-dire la convergence de la solution de l’équation (3) (resp. (4)) vers la solution de équation (1) (resp. (2)) a été démontrée dans les paragraphes précédents mais l’ordre de convergence en ε n’est pas clair et la forte nonlinéarité du problème rend cette analyse particulièrement difficile. En ce qui concerne l’erreur numérique, elle fait intervenir généralement le nombre, N , de particules et le paramètre ε, la fonction ζε jouant le rôle de fonction de forme de particules. Sans préciser d’avantage les hypothèses, disons que si la fonction ζ et la condition initiale u0 sont assez régulières, on obtient une erreur du type
1 u − uh = O u − uh = O 1 , + O εr , ε m m N (εN )
où u est la solution de l’équation (1) (resp. (2)), uh = p∈N αp δ(x − Xp (t)) son approximation par des
mesures de Dirac, uhε = p∈N αp ζε (x − Xp (t)) sa version régularisée et les positions Xp (t) sont les
373
P.-L. Lions, S. Mas-Gallic
positions des particules. Les normes sont des normes d’espaces de Sobolev habituels et les exposants m et r dépendent, en particulier des fonctions ζ et u0 . La seconde estimation fait apparaître, pour un nombre N de particules donné, l’existence d’une valeur optimale de ε, qui peut être recherchée lors de tests numériques. Les illustrations numériques monodimensionnelles qui suivent confirment, d’une part dans le cas de l’équation des milieux poreux et d’autre part dans le cas de l’équation de la chaleur, les résultats de l’analyse théorique. Les calculs ont été fait avec un nombre fixé de 250 particules, 100 pas de temps et une fonction cut-off P2 par morceaux de classe C1 ; les erreurs sont calculées dans les normes pour lesquelles la convergence a été obtenue. 4.1. Équation de type milieu poreux. – La fonction u0 étant donnée, on résout numériquement l’équation ∂u = div(Vε u), ∂t
Vε = ∇u ∗ ζε
où Vε , vitesse de diffusion définie précédemment a été réécrite en tenant compte des propriétés du produit de convolution. Ceci nous permet de l’appliquer à l’approximation particulaire uh de u et ainsi de définir la vitesse des particules. La position des particules au temps t est donnée par la résolution du système d’équations différentielles dXp (t) = − αq ∇ζε Xp (t) − Xq (t) , dt
Xp (0) = x0p ,
q∈N
auquel un schéma de Runge–Kutta d’ordre deux a été appliqué. Les courbes de la figure 1 montrent que, pour ε = 0.20, la solution autosimilaire dite de Barenblatt est bien calculée ; on constate néanmoins que les points d’intersection de la parabole avec l’axe n’ont pas tout à fait la bonne vitesse sans doute parce que les particules initialement placées en ces deux points n’y sont plus. Pour résoudre cette difficulté, plusieurs solutions sont envisageables. La courbe d’erreur en fonction du paramètre ε, figure 2, montre d’une part sur la courbe de gauche, l’existence d’une valeur optimale de ε recherchée, elle est ici obtenue pour ε = 0.10 et donne d’autre part, sur la courbe de droite, une estimation de l’ordre de convergence en ε qui ici à peu près 1.6. Précisons enfin que le pas de temps ∆t est, comme cela été établi dans [2], lié à la valeur de ε et non pas au nombre N de particules et qu’une étude heuristique nous à fait choisir ici ∆t = ε. 4.2. Équation de la chaleur. – Soit u0 une fonction donnée. Rappelons tout d’abord l’équation à résoudre respectivement ∂u = div Vε u , ∂t
∇u ∗ ζε Vε = , u ∗ ζε
Figure 1. – À gauche, condition initiale ; à droite, solutions approchée (trait plein) et exacte (trait pointillé). Figure 1. – Left, initial condition; right, approximate (plain line) and exact (doted line) solutions.
374
(11)
Une méthode particulaire déterministe pour des équations diffusives non linéaires
Figure 2. – À gauche, erreur en norme L1 en x, L2 en t en fonction du paramètre ε ; à droite, zoom et droite de plus proche pente donnant l’ordre de convergence. Figure 2. – Left, L1 in x, L2 in t error with respect to ε; right, closer view plus nearest line with order of convergence.
Figure 3. – À gauche, condition initiale ; à droite, solutions approchée (trait plein) et exacte (trait pointillé). Figure 3. – Left, initial condition; right, approximate (plain line) and exact (doted line) solutions.
Figure 4. – Erreur en norme L1 ln L1 en x, L2 en t en fonction du paramètre ε. Figure 4. – L1 in x, L2 in t error with respect to ε.
où Vε est la vitesse de diffusion. Cette écriture nous permet de définir la vitesse des particules dont la position au temps t est donnée par la résolution du système d’équations différentielles
dXp q∈N αq ∇ζε (Xp (t) − Xq (t)) (t) = − , Xp (0) = x0p , dt q∈N αq ζε (Xp (t) − Xq (t)) qui est résolu par un schéma de Runge–Kutta d’ordre deux.
375
P.-L. Lions, S. Mas-Gallic
Les courbes, figure 3, montrent que la Gaussienne, solution de l’équation de la chaleur, est bien calculée, la courbe d’erreur en fonction de ε, figure 4, montrant une convergence moins rapide que dans le cas de l’équation précédente. Précisons également, que le pas de temps, ∆t, qui est calculé en fonction de la valeur de ε a été fixé, après une étude heuristique, à ∆t = ε2 /2. Remarque 3. – On peut avoir des problèmes de calcul si la solution est proche de zéro. Pour pallier cette difficulté, un seuil numérique η > 0 est introduit et le système couplé d’équations différentielles ordinaires résolu est en fait le suivant dXp q∈N αq ∇ζε (Xp (t) − Xq (t)) q∈N αq ζε (Xp (t) − Xq (t)) (t) = − , Xp (0) = x0p . 2 2 dt αq ζε (Xp (t) − Xq (t)) + η q∈N
Conclusion. – La méthode particulaire qui a été utilisée ici, est, comme le sont les méthodes de marche aléatoire, une méthode particulaire à poids fixes ce qui la rend beaucoup moins coûteuse en temps de calcul que, par exemple, la méthode à poids variable (communément appelée P.S.E., particle strength exchange). Elle est clairement très dépendante de la fonction ζ utilisée et, au cours des tests trois fonctions ζ différentes ont été utilisées, le noyau dit de Fejer, une fraction rationnelle classique en méthode de vortex et la fonction qui nous a servi dans les résultats présentés ici. Le choix de cette dernière fonction est un choix lié à la rapidité et à la précision des calculs. Les démonstrations détaillées et des résultats numériques plus complets feront l’objet d’un prochain article. Références bibliographiques [1] Degond P., Mustieles F.J., Approximation of diffusion equations by deterministic convections of particles, SIAM J. Sci. Statis. Comput. 11 (1990) 293–310. [2] Degond P., Mas-Gallic S., The weighted particle method for convection–diffusion equations, part I: the case of an isotropic viscosity, part II: the anisotropic case, Math. Comput. 53 (1989) 485–526. [3] Fronteau J., Combis X., Hadronic J. 7 (1984) 911–940. [4] Grmela M., Fronteau J., Tellez-Arenas A., Hadronic J. 3 (1980) 1203–1226. [5] Huberson S., Rivoalen E., Hauville F., Simulation numérique des équations de Navier–Stokes 3d par une méthode particulaire, C. R. Acad. Sci. Paris, Série IIb 324 (1997) 543–549. [6] Huberson S., Rivoalen E., Numerical simulation of axisymmetric viscous flows by means of particle method, J. Comput. Phys. 152 (1) (1999) 1–31. [7] Lacombe G., Analyse d’une équation de vitesse de diffusion, C. R. Acad. Sci. Paris, Série I 329 (5) (1999) 383– 386. [8] Lacombe G, Mas-Gallic S., Presentation and analysis of a diffusion velocity method, in: ESAIM Proceedings, Vortex Flow and Related Numerical Methods III, 1998. [9] Mas-Gallic S., Une présentation de la méthode de vitesse diffusion, Conférence en l’honneur de Marc Feix, Orléans, in: Dynamical Systems, Plasmas and Gravitation, Leach P.G.L., Bouquet S.E., Rouet J.L., Fijalkow E. (Eds.), Springer-Verlag, 1997, pp. 74–81. [10] Mas-Gallic S., Raviart P.-A., A particle method for first order symmetric systems, Numer. Math. 51 (1987) 323– 352. [11] Mustieles F.J., Thèse de doctorat, CMAP, École polytechnique, Palaiseau, 1990. [12] Ogami Y., Akamatsu T., Viscous flow simulation using a discrete vortex method, J. Computers and Fluids 19 (1991) 433–441. [13] Rivoalen E., Thèse de doctorat, Le Havre, 1994. [14] Strickland J.H., Kempka S.N., Wolfe W.P., Viscous diffusion of vorticity using the diffusion velocity concept, in: ESAIM Proceedings, Vortex Flow and Related Numerical Methods II, Vol. 1, 1996, pp. 135–151.
376