6
C. R. Acad. Sci. Paris, t. 330, Série I, p. 517–522, 2000 Analyse numérique/Numerical Analysis Méthode spectrale rapide pour l’équation de Fokker–Planck–Landau Lorenzo PARESCHI a , Giovanni RUSSO b , Giuseppe TOSCANI c a Départment de mathématiques, Université de Ferrara, via Machiavelli 35, 44100 Ferrara, Italie Courriel : [email protected] b Départment de mathématiques, Université de L’Aquila, via Vetoio, loc. Coppito, 67010 L’Aquila, Italie Courriel : [email protected] c Départment de mathématiques, Université de Pavia, via Ferrata 1, 27100 Pavia, Italie Courriel : [email protected] (Reçu le 20 juillet 1999, accepté après révision le 31 janvier 2000) Résumé. Dans cette Note, nous introduisons une nouvelle méthode spectrale pour l’équation de Fokker–Planck–Landau (FPL). La méthode permet, contrairement au coût usuel en O(n 2 ) de l’opérateur de collision de FPL, d’obtenir des solutions numériques spectralement précises avec simplement O(n log 2 n) opérations. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS Fast spectral methods for the Fokker–Planck–Landau equation Abstract. In this Note we present a new spectral method for the Fokker–Planck–Landau (FPL) equation. The method allows to obtain spectrally accurate numerical solutions with simply O(n log 2 n) operations in contrast with the usual O(n 2 ) cost of the FPL collision operator. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS Abridged English version We will adapt to the FPL equation a new method for the numerical solution of the space homogeneous Boltzmann equation, based on a Fourier spectral approximation of the collision operator in the velocity space, recently introduced in [9,10]. In the FPL case, the method gives a computational cost of order O(n log 2 n), where n is the number of Fourier modes, in contrast with the usual quadratic cost of deterministic methods, and it can be designed to preserve mass and to approximate with spectral accuracy momentum and energy. Note présentée par Philippe G. CIARLET. S0764-4442(00)00212-3/FLA 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés. 517

Méthode spectrale rapide pour l'équation de Fokker–Planck–Landau

Embed Size (px)

Citation preview

C. R. Acad. Sci. Paris, t. 330, Série I, p. 517–522, 2000Analyse numérique/Numerical Analysis

Méthode spectrale rapide pour l’équationde Fokker–Planck–LandauLorenzo PARESCHI a, Giovanni RUSSOb, Giuseppe TOSCANIc

a Départment de mathématiques, Université de Ferrara, via Machiavelli 35, 44100 Ferrara, ItalieCourriel : [email protected]

b Départment de mathématiques, Université de L’Aquila, via Vetoio, loc. Coppito, 67010 L’Aquila, ItalieCourriel : [email protected]

c Départment de mathématiques, Université de Pavia, via Ferrata 1, 27100 Pavia, ItalieCourriel : [email protected]

(Reçu le 20 juillet 1999, accepté après révision le 31 janvier 2000)

Résumé. Dans cette Note, nous introduisons une nouvelle méthode spectrale pour l’équation deFokker–Planck–Landau (FPL). La méthode permet, contrairement au coût usuel enO(n2)de l’opérateur de collision de FPL, d’obtenir des solutions numériques spectralementprécises avec simplementO(n log2 n) opérations. 2000 Académie des sciences/Éditionsscientifiques et médicales Elsevier SAS

Fast spectral methods for the Fokker–Planck–Landau equation

Abstract. In this Note we present a new spectral method for the Fokker–Planck–Landau(FPL)equation. The method allows to obtain spectrally accurate numerical solutions with simplyO(n log2 n) operations in contrast with the usualO(n2) cost of the FPL collision operator. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS

Abridged English version

We will adapt to the FPL equation a new method for the numerical solution of the space homogeneousBoltzmann equation, based on a Fourier spectral approximation of the collision operator in the velocityspace, recently introduced in [9,10]. In the FPL case, the method gives a computational cost of orderO(n log2 n), wheren is the number of Fourier modes, in contrast with the usual quadratic cost ofdeterministic methods, and it can be designed to preserve mass and to approximate with spectral accuracymomentum and energy.

Note présentée par Philippe G. CIARLET .

S0764-4442(00)00212-3/FLA 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés. 517

L. Pareschi

1. Introduction

Nous considérons l’équation de Landau–Fokker–Planck (FPL) pour la fonction de distributionf , quidépend de la vitessev dansRd, avecd> 2, et du tempst dansR+ :

∂f

∂t=QL(f, f),

QL(f, f)(v) =∇ ·∫

dv∗A(v − v∗)(f(v∗)∇f(v)− f(v)∇v∗f(v∗)

). (1.1)

Dans (1.1), la matriceA est symétrique, semi-définie positive, et dépend de l’interaction entre les particules,A(z) = Ψ(|z|)Π(z) pour une fonction positiveΨ. Pour les interactions décroissant selon une puissance dela distance :

Πij(z) =

(δij −

zizj|z|2

), Ψ(|z|) = Λ|z|γ+2, γ >−3,

où Λ est une constante positive. Les cas les plus intéressants pour les applications sont les casγ = −3,d= 3, i.e. le cas coulombien.

Plusieurs approches numériques de l’équation de FPL ont été développées durant les dix dernières années[1–3,5,7,8]. Ici, nous adaptons à l’équation de FPL une nouvelle méthode pour les solutions numériquesde l’équation de Boltzmann homogène en espace, basée sur une approximation spectrale de l’opérateur decollision dans l’espace des vitesses, récemment introduite dans [9,10]. Dans le cas de l’équation de FPL, laméthode a un coût enO(n log2 n), oùn est le nombre de modes de Fourier, contrairement au coût usuel enO(n2), et elle peut préserver la masse et approcher avec une précision spectrale la quantité de mouvementet l’énergie. Finalement, signalons que l’on peut obtenir une réduction comparable du coût de calcul enutilisant les techniques de [3,7].

2. Extrapolation spectrale de l’équation de FPL

Soitg = v− v∗. L’opérateur de collision de FPL (1.1) peut s’écrire comme :

QL(f, f)(v) =∇ ·∫

dgA(g)(f(v+ g)∇f(v)− f(v)∇gf(v+ g)

).

SoitB(0,R) la boule de rayonR de centre l’origine. Sisupp (f(v))⊂B(0,R), alors

QL(f, f)(v) =∇ ·∫B(0,2R)

A(g)(f(v+ g)∇f(v)− f(v)∇gf(v+ g)

)dg, (2.1)

avecv + g ∈ B(0,3R).Dans le reste de l’article, nous restreignons la fonction de distribution sur le cube[−T,T ]3 avecT > 3R,

supposant quef(v) = 0 sur[−T,T ]3rB(0,R), et en la prolongeant par périodicité sur[−T,T ]3. Commeobservé dans [10], il suffit de prendreT > 2R pour éviter les intersections des regions oùf est différent dezero. Pour simplifier les formules, on suppose queT = π, et doncR= π/2.

La fonction de distributionf est approchée par une somme partielle de séries de Fourier

fN(v) =N∑

k=−Nfk eik·v , fk =

1

(2π)3

∫Ω

f(v) e−ik·v dv, (2.2)

518

Méthode spectrale rapide pour l’équation de Fokker–Planck–Landau

oùΩ = [−π,π]3, k = (k1, . . . , kd), et, pour simplifier les notations,

N∑k=−N

≡N∑

k1=−N· · ·

N∑kd=−N

.

Un système d’EDO’s pour les coefficientsfk, k =−N, . . . ,N , est obtenu à partir de [4,6] :∫Ω

(∂fN∂t−QL(fN , fN )

)e−ik·v dv = 0. (2.3)

En substituant l’expression (2.2) dans (2.3), et en utilisant les propriétés orthogonales des polynômestrigonométriques, nous obtenons le procédé

dfkdt

=

N∑`+m=k`,m=−N

f` fmβL(`,m), k =−N, . . . ,N, (2.4)

où βL(`,m) = BL(`,m)− BL(m,m) et lesmodes du noyaude FPLBL(`,m) sont donnés par :

BL(`,m) =

∫B(0,π)

Ψ(|g|)[`2 −

(` · g|g|

)2]eig·m dg. (2.5)

Remarquons que (2.5) sont des quantités scalaires indépendantes de la fonctionf , avec la propriété

BL(`,m) = BL(−`,m) = BL(`,−m) = BL(`,−m),

BL(m,m) dépend seulement de|m|.(2.6)

La dernière égalité de l’équation (2.6) montre que les coefficients sont réels, et est une conséquence dufait que la fonctionΨ dépend seulement de|z|. De plus, si le noyauΨ(|g|) = Λ|g|2+γ , alors∣∣BL(`,m)

∣∣6 `2BL(0,0),

BL(0,0) = 4ππγ+5

γ + 5Λ.

3. Propriétés principales

Pour chaquet> 0, fN(v, t) est un polynôme trigonométrique de degréN env, i.e.fN (t) ∈ PN , où

PN = span

eik·v | −N 6 kj 6N, j = 1,2,3.

De plus, soitPN : L2(Ω)→ PN la projection orthogonale surPN pour le produit scalaire deL2(Ω)(cf. (2.3)) : ⟨

f −PNf,φ⟩

= 0, ∀φ ∈ PN .

519

L. Pareschi

Avec cette définitionPNf = fN , oùfN est la série de Fourier partielle def (2.2) et la méthode définiepar (2.4) peut s’écrire sous la forme équivalente :

∂fN∂t

=QLN(fN , fN ),

où QLN(fN , fN) = PNQL(fN , fN) et QL(·, ·) est l’opérateur de FPL avec uncut-off sur la vitesserelative (2.1). Il est facile de vérifier d’après la propriété (2.6) que l’opérateur de collision extrapoléQLN(fN , fN) préserve la masse dans le temps, i.e.

d

dt

∫Ω

fN (v) dv = 0.

Maintenant, si nous notonsHrp(Ω), oùr > 0 est un entier, le sous-espace de l’espace de SobolevHr(Ω),

qui est formé des fonctions périodiques [4,6], nous avons [12]

PROPOSITION 3.1. –Soitf ∈H2p(Ω), g ∈ L2(Ω), alors∥∥QL(f, g)

∥∥26C ‖g‖1 ‖f‖H2

p.

Pour simplifier les notations, nous noteronsQL(f) au lieu deQL(f, f). La qualité de l’approximationen normeL2 de l’opérateur de collision de FPLQL(f) parQLN(fN ), est donnée par :

THÉORÈME 3.2. –Soitf ∈H2p(Ω), alors pour toutr > 0,

∥∥QL(f)−QLN(fN )∥∥

26C

(‖f − fN‖H2

p+‖QL(fN )‖Hrp

N r

),

oùC depend de‖f‖2.

Démonstration. –D’abord, nous pouvons diviser l’erreur en deux parties∥∥QL(f)−QLN(fN )∥∥

26∥∥QL(f)−QL(fN )

∥∥2

+∥∥QL(fN )−QLN(fN )

∥∥2.

Maintenant, on a clairementQL(fN ) ∈ P2N et doncQL(fN) est périodique et régulier ainsi que sesdérivées comme dans [4] :∥∥QL(fN )−QLN(fN )

∥∥26 C

N r

∥∥QL(fN )∥∥

Hrp, ∀r > 0.

Par application de la proposition 3.1 et de l’identité

QL(f)−QL(g) =QL(f + g, f − g),

nous avons ∥∥QL(f)−QL(fN )∥∥

2=∥∥QL(f + fN , f − fN )

∥∥2

6C ‖f + fN‖1‖f − fN‖H2p6 2C1‖f‖2‖f − fN‖H2

p. 2

Enfin, il suffit de voir que

‖f − fN‖H2p6 C

N r−2‖f‖Hrp,

pour avoir la précision spectrale de l’approximation.

520

Méthode spectrale rapide pour l’équation de Fokker–Planck–Landau

COROLLAIRE 3.3. –Soitf ∈Hrp(Ω), r > 2, alors

∥∥QL(f)−QLN(fN)∥∥

26 C

N r−2

(‖f‖Hrp +

∥∥QL(fN )∥∥

Hrp

).

D’après le corollaire précédent, il découle que

∣∣⟨QL(f), ϕ⟩−⟨QLN(fN), ϕ

⟩∣∣6 C

N r−2‖ϕ‖2

(‖f‖Hrp +

∥∥QL(fN )∥∥

Hrp

),

et donc, en prenantϕ= v, v2, on obtient la précision spectrale pour la quantité de mouvement et l’énergie.

4. Mise en place de l’algorithme rapide

Réécrivons le schéma (2.4) sous la forme,k =−N, . . . ,N :

dfkdt

=

N∑m=−N

fk−m fmβL(k−m,m). (4.1)

Dans l’expression précédente, nous supposons que les coefficients de Fourier sont nuls pour|kj | > N ,j = 1,2,3. L’évaluation suivante de (4.1) demande exactementO(n2) opérations, oùn= Nd et d> 2 estla dimension de l’espace des vitesses. Cependant, grâce à la structure des modes du noyau de FPL, il estpossible de réduire le calcul à seulementO(n log2n) opérations. En fait,BL(`,m) se sépare en

BL(`,m) = `2F (m)−d∑

p, q=1

`p `q Ipq(m) = `2 F (m)− `I(m) `T , (4.2)

où `T représente le vecteur`= (`1, `2, `3), I = (Ipq) est une matrice symétrique3× 3, et

F (m) =

∫B(0,π)

|g|2+γ eig·m dg, (4.3)

Ipq(m) =

∫B(0,π)

|g|γ gp gq eig·m dg, p, q = 1, . . . , d.

Nous pouvons donc écrire

βL(`,m) = `2 F (m)− `I(m) `T − BL(m,m).

L’algorithme demande donc le calcul ded(d+ 1)/2 + 1 convolutions pour le terme de gain (le nombred’éléments distincts deI plus un) et une seule convolution pour le terme de perte. Il est bien connu qu’ilest possible de calculer une convolution en seulementO(n log2 n) opérations. Pour les détails de la mise enplace de cette technique standard pour des méthodes spectrales, nous invitons le lecteur à se reporter à [4].

Maintenant, soitΨ(|g|) = Λ|g|γ+2. Pour simplifier les notations, nous prenonsΛ = 1. Pour la mise enplace de l’algorithme, nous devons évaluer les quantités (4.2)–(4.3). Par souci de simplicité, nous traitonsici seulement le cas de la dimension2 : d= 2. Un calcul simple montre que

F (m) = F(|m|)

= 2π

∫ π

0

rγ+3J0

(|m|r

)dr,

oùJ0 est la fonction de Bessel d’ordre0. De plus, on peut montrer que :

521

L. Pareschi

I11(m) =1

2

[F(|m|)

+m2

1 −m22

|m|2 G(|m|)],

I22(m) =1

2

[F(|m|)− m2

1 −m22

|m|2 G(|m|)],

I12(m) = I21(m) =m1m2

|m|2 G(|m|),

G(|m|) = 2π

∫ π

0

rγ+3C(|m|r

)dr,

et

C(x) =1

∫ 2π

0

cos(x cosφ) cos(2φ) dφ.

Donc le calcul des modes du noyau de FPL en deux dimensions réduit simplement le calcul des intégralesà une variableF (|m|) etG(|m|). Ces quantités peuvent être calculées très précisément et stockées dans destableaux à deux dimensions.

Références bibliographiques

[1] Berezin Yu.A., KhudickV.N., Pekker M.S., Conservative finite difference schemes for the Fokker–Planck equationnot violating the law of an increasing entropy, J. Comput. Phys. 69 (1987) 163–174.

[2] Buet C., Cordier S., Conservative and entropy decaying numerical scheme for the isotropic Fokker–Planck–Landau equation, J. Comput. Phys. 145 (1) (1998) 228–245.

[3] Buet C., Cordier S., Degond P., Lemou M., Fast algorithms for numerical, conservative, and entropy approxima-tions of the Fokker–Planck equation, J. Comput. Phys. 133 (1997) 310–3229.

[4] Canuto C., Hussaini M.Y., Quarteroni A., Zang T.A., Spectral Methods in Fluid Dynamics, Springer-Verlag, NewYork, 1988.

[5] Degond P., Lucquin-Desreux B., An entropy scheme for the Fokker–Planck collision operator of plasma kinetictheory, Numer. Math. 68 (1994) 239–262.

[6] Gottlieb D., Orszag S.A., Numerical Analysis of Spectral Methods: Theory and Applications, SIAM CBMS-NSFSeries, 1977.

[7] Lemou M., Multipole expansions for the Fokker–Planck–Landau operator, Numer. Math. 78 (1998) 597–618.[8] Lucquin-Desreux B., Discrétization de l’opérateur de Fokker–Planck dans le cas homogène, C. R. Acad. Sci. Paris,

Série I 314 (1992) 407–411.[9] Pareschi L., Perthame B., A Fourier spectral method for homogeneous Boltzmann equations, Transp. Theo. Stat.

Phys. 25 (1996) 369–383.[10] Pareschi L., Russo G., Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the

collision operator, SIAM J. Numer. Anal. (to appear).[11] Pareschi L., G. Russo, Toscani G., Fast spectral methods for the Fokker–Planck–Landau collision operator of

plasma physics, J. Comput. Phys. (2000) (submitted).[12] Pareschi L., Toscani G., Villani C., Spectral methods for the non-cut-off Boltzmann equation and numerical

grazing collision limit, Preprint, 1999.

522