7 Model reduction of a kinetic swarming model by 1 0 2 operator projection n a J Junming Duan∗, Yangyu Kuang∗, Huazhong Tang† 1 1 January 12, 2017 ] A N . Abstract h t This paper derives the arbitrary order globally hyperbolic moment system a m for a non-linear kinetic description of the Vicsek swarming model by using the [ operatorprojection. Itisbuiltonourcarefulstudyofafamilyofthecomplicate Gradtypeorthogonalfunctionsdependingonaparameter(angleofmacroscopic 1 v velocity). We calculate their derivatives with respect to the independent vari- 8 able,andprojectionofthosederivatives,theproductofvelocityandbasis,and 8 collisionterm. Themomentsystemisalsoprovedtobehyperbolic,rotationalin- 8 variant,andmass-conservative. TherelationshipbetweenGradtypeexpansions 2 0 in different parameter is also established. A semi-implicit numerical scheme is 1. presentedtosolveaCauchyproblemofourhyperbolicmomentsysteminorder 0 to verify the convergence behavior of the moment method. It is also compared 7 tothespectralmethodforthekineticequation. Theresultsshowthatthesolu- 1 : tions of our hyperbolic moment system converge to the solutions of the kinetic v equation for the Vicsek model as the order of the moment system increases, i X andthemomentmethodcancapturekeyfeaturessuchasvortexformationand r traveling waves. a Keywords: Moment method, hyperbolicity, kinetic equation, modelreduc- tion, operator projection 1 Introduction Swarm behaviour, or swarming, is a collective behaviour exhibited by entities, par- ticularly animals, of similar size which aggregate together, perhaps milling about the samespotorperhapsmovingenmasseormigratinginsomedirection[1]. Itisahighly ∗LMAM,SchoolofMathematicalSciences,PekingUniversity,Beijing100871,P.R.China. Email: [email protected];[email protected] †HEDPS,CAPT&LMAM,SchoolofMathematicalSciences,PekingUniversity,Beijing100871, P.R. China; School of Mathematics and Computational Science, Xiangtan University, Xiangtan 411105,HunanProvince,P.R.China. Email: [email protected] 1 interdisciplinary topic. Some works studied kinetic models for swarming [7, 10, 17], but few did a numerical investigation. The first numerical method was presented in [13] for a kinetic description of the Vicsek swarming model. The main contribution wastouseaspectralrepresentationlinkedwithadiscreteconstrainedoptimizationto compute those interactions. Unfortunately, only first-order accurate upwind method was used to approximated the transport term. Thekinetictheoryhasbeenwidelystudiedandplayedanimportantroleinmany fields during several decades, see e.g. [8, 9]. The kinetic equation can determine the distribution function hence the transport coefficients, however such task is not so easy. The moment method [15, 16] is a model reduction for the kinetic equation by expanding the distribution function in terms of tensorial Hermite polynomials and introducing the balance equations corresponding to higher order moments of the distribution function. One major disadvantage of the Grad moment method is the lossofhyperbolicity,whichwillcausethesolutionblow-upwhenthedistributionisfar away from the equilibrium state. Increasing the number of moments could not avoid suchblow-up[6]. Uptonow,therehasbeensomelatestprogressontheGradmoment method for the kinetic equation. Numerical regularized moment method of arbitrary orderwasstudiedforBoltzmann-BGKequation[5]andforhighMachnumberflow[6]. BasedontheobservationthatthecharacteristicpolynomialofthefluxJacobianinthe Grad moment system did not depend on the intermediate moments, a regularization waspresentedin[2,3,4]fortheone-andmulti-dimensionalGradmomentsystemsto achieve global hyperbolicity. The quadrature based projection methods were used to derivehyperbolicsystemsforthesolutionoftheBoltzmannequation[18,19]byusing thequadratureruleinsteadoftheexactintegration. Inthe1Dcase,itissimilartothe regularization in [2]. Those contributions led to well understanding the hyperbolicity oftheGradmomentsystems. Basedontheoperatorprojection, ageneralframework of model reduction technique was recently given in [12]. It projected the time and spacederivativesinthekineticequationintoafinite-dimensionalweightedpolynomial space synchronously, and could give most of the existing moment systems in the literature. Recently, such model reduction method was also successfully extended to the 1D special relativistic Boltzmann equation and the globally hyperbolic moment model of arbitrary order was derived in [20]. The aim of this paper is to extend the model reduction method by the operator projection to the two-dimensional kinetic description of the Vicsek swarming model and derive corresponding globally hyperbolic moment system of arbitrary order. The paper is organized as follows. Section 2 introduces the kinetic and macroscopic equa- tionsfortheVicsekmodel. Section3givesafamilyoforthogonalfunctionsdependent onaparameterandtheirpropertiesandderivesthearbitraryordergloballyhyperbolic momentsystemofthekineticdescriptionfortheVicsekmodel. Section4investigates the mathematical properties of moment system, including: hyperbolicity, rotational 2 invariance, mass-conservation and relationship between Grad type expansions in dif- ferent parameter. Section 5 presents a semi-implicit numerical scheme and conducts a numerical experiment to check the convergence of the proposed hyperbolic moment system. Section 6 concludes the paper. 2 Kinetic and macroscopic equations for the Vicsek model This section introduces the non-linear kinetic and macroscopic equations for the Vic- sek model. 2.1 Kinetic equation The kinetic equation for the Vicsek model can be written as follows [13] ∂tf +ω·∇xf +∇ω ·(F[f]f)=σ∆ωf, (2.1) where f(t,x,ω) is the particle distribution function depending on the time t, spatial variable x∈Rd, and the unit velocity vector ω, the parameter σ is a scaled diffusion constant describing the intensity of the noise with the Brownian motion, the vector F[f](t,x,ω) is the mean-field interaction force between the particles and given by J(t,x) F[f](t,x,ω)=(Id−ω⊗ω)Ω(x), Ω(t,x)= , (2.2) |J(t,x)| the notation Id is the identity operator, Ω(t,x) is the mean velocity, and J(t,x) denotes the mean flux at x and is defined by (cid:90) (cid:90) J(t,x)= K(y−x)ωf(t,y,ω)dydω. (2.3) Rd Sd−1 Here K(y −x) is the characteristic function of the ball B(0,R) = {x : |x| ≤ R}, i.e. K(x) = 1|x|<R, and R is the radius of the alignment interactions between the particles. The vector field F[f](t,x,ω) tends to align the particles to the direction Ω which is the director of the particle flux J and becomes spatially local in the large scale limit of space and time so that J can be approximated by (cid:90) J(t,x)= ωf(t,x,ω)dω. (2.4) Sd−1 In the numerical computations, J(t,x) can be considered within the approximation (2.4)becausethespatialmeshstepsizemaybesmallerthantheradiusofballB(0,R) [13]. It has been shown [14] that the non-linear kinetic equation (2.1) with (2.2) and (2.3) or(2.4) has anon-negative global weak solutionin C(0,T;L1(D)∩L∞(cid:0)(0,T)× 3 D(cid:1) with D = Rd×Sd−1, for any time T, given non-negative initial value f(0,x,ω) in L1(D)∩L∞(D) and J(t,x) which is always not equal to 0 for t∈[0,T]. Throughout the paper, we will only consider the case of (2.4) and the direction of mean velocity of f becomes (cid:82) ωfdω Ω= (cid:82)Sd−1 . | ωfdω| Sd−1 Moreover, the kinetic equation (2.1) is rewritten as following ∂ f +ω·∇ f =Q(f), (2.5) t x where the collision term Q(f) is defined by Q(f)=−∇ω ·(F[f]f)+σ∆ωf. (2.6) Lemma 1 ([13]). The collision term Q(f) can be rewritten as follows (cid:18) (cid:19) (cid:0) f (cid:1) Q(f)=σ∇ω · MΩ∇ω , (2.7) M Ω and satisfies (cid:90) dω Q(f)f ≤0, (2.8) M Sd−1 Ω where ω·Ω MΩ(ω)=C0exp( σ ), (2.9) is the equilibrium function, also known as the Von Mises distribution, and C is a 0 constant of normalization. Remark 2.1. The equilibria of operator Q are given by the set {ρMΩ|ρ ∈ R,Ω ∈ Sd−1} and forms a set of dimension d, but the collisional invariants of Q are only of dimension 1. Particularly, Q does preserve the mass but cannot preserve the flux, that is, for a general f, it holds (cid:90) (cid:90) Q(f)dω =0, Q(f)ωdω (cid:54)=0. Sd−1 Sd−1 Remark 2.2. In the case of d=2, if letting Ω=(cosθ¯,sinθ¯)T, then the equilibrium can be expressed as follows (cid:18)cos(θ−θ¯)(cid:19) M(θ−θ¯)=C exp , (2.10) 0 σ and the kinetic equation (2.5) becomes f +cosθf +sinθf =Q(f), t x y Q(f)=σ∂θ(cid:18)Mθ¯∂θ(cid:0)Mfθ¯(cid:1)(cid:19)=∂θ(sin(θ−θ¯)f)+σ∂θ2f, (2.11) where θ and θ¯ denote the angles of microscopic and macroscopic velocities, respec- tively. 4 2.2 Macroscopic equations The kinetic equation (2.5) is written at the microscopic level, i.e. at time and length scales which are characteristic of the dynamics of the individual particles. When investigating the dynamics of the system at large time and length scales compared with the scales of the individuals, a set of new variables t˜=εt and x˜ =εx has to be introduced[11], whereεdenotestheratiobetweenmicroandmacrovariables, ε(cid:28)1. In this new set of variables, the kinetic equation (2.5) is written (after dropping the tildes for clarity) as follows 1 ∂ fε+ω·∇ fε = Q(fε). (2.12) t x ε In the limit ε→0, fε converges locally in space to an equilibrium state in the local space, that is fε →f(t,x,ω)=ρ(t,x)MΩ(ω). (2.13) (cid:82) At this moment, the evolution of macroscopic density ρ = fdω and velocity Ω Sd−1 is described by the following equations ∂ ρ+∇ ·(c ρΩ)=0, (2.14) t x 1 ρ(∂ Ω+c (Ω·∇ )Ω)+λ(Id−Ω⊗Ω)∇ ρ=0, (2.15) t 2 x x |Ω|=1, (2.16) where c ,c , and λ are three constants depending on σ. Such macroscopic system is 1 2 hyperbolic but non-conservative, and the operator Id−Ω⊗Ω ensures the constraint (2.16). Remark 2.3. In the 2D case, the macroscopic density ρ and angle of velocity θ¯are related to the particle distribution function f by (cid:90) 2π ρ(t,x)= f(t,x,θ)dθ, (2.17a) 0 (cid:32) (cid:33) (cid:90) 2π cosθ J(t,x)= f(t,x,θ)dθ, (2.17b) sinθ 0 J(t,x) Ω=(cosθ¯,sinθ¯)T = . (2.17c) |J(t,x)| 3 Derivation of moment system This section derives the moment system for the 2D kinetic description of Vicsek swarming model by using the operator projection [12, 20]. For the sake of simplicity, we assume θ¯= 0 in Sections 3.1 and 3.2. In fact, the general case of θ¯(cid:54)= 0 can be converted into such simple case by operating a translation transformation θ (cid:55)→θ−θ¯. 5 3.1 Orthogonal functions The Hermite polynomials are used to derive the Grad’s moment system [15, 16], where the distribution function is assumed to close to a local Maxwellian. Here we need to find the orthogonal functions with respect to the weight M(θ) defined in [0,2π], because θ ∈[0,2π]. Those functions, denoted by Hc(θ),Hc(θ),··· ,Hc (θ),··· ,Hs(θ),··· ,Hs(θ),··· , 0 1 N 1 N are built on the trigonometric functions {1,cosθ,··· ,cos(Nθ),··· ;sinθ,··· ,sin(Nθ),···} byusingtheSchmitorthogonalprocess,wherethesuperscriptc(resp. s)denotesthe functions consisting of a linear combination of cos(kθ) (resp. sin(kθ)). Because the function M(θ) is even, it holds (cid:90) 2π sin(mθ)cos(nθ)M(θ)dθ =0, ∀m,n∈N, 0 so in the process of Schimidt orthogonalization, the linear combination of cos(kθ) is orthogonal to that of sin(kθ), so that there are two sets of “independent” or- thogonal functions. The first 2N + 1 functions are Hc(θ),Hc(θ),··· ,Hc (θ) and 0 1 N Hs(θ),··· ,Hs(θ), and can be expressed in terms of the trigonometric functions 1 N {1,cosθ,··· ,cos(Nθ), sinθ,··· ,sin(Nθ)} as follows (cid:16) (cid:17)T (cid:16) (cid:17)T Hc(θ),Hc(θ),··· ,Hc (θ) =A 1,cosθ,··· ,cos(Nθ) , (3.1) 0 1 N N (cid:16) (cid:17)T (cid:16) (cid:17)T Hs(θ),··· ,Hs(θ) =B sinθ,··· ,sin(Nθ) , (3.2) 1 N N where A ∈ R(N+1)×(N+1),B ∈ RN×N, and both A and B are invertible lower N N N N triangular matrix. This set of functions satisfy (cid:90) 2π Hl (θ)Hl(θ)M(θ)dθ =δ , m,n∈{0,1,··· ,N},l∈{c,s}. m n m,n 0 When l=s, m,n(cid:54)=0. For the sake of simplicity, define (cid:16) (cid:17)T (cid:16) (cid:17)T 1,cosθ,··· ,cos(Nθ) (cid:44)Ec (θ), sinθ,··· ,sin(Nθ) (cid:44)Es (θ), N N (cid:16) (cid:17)T (cid:16) (cid:17)T Hc(θ),Hc(θ),··· ,Hc (θ) (cid:44)Hc (θ), Hs(θ),Hs(θ),··· ,Hs(θ) (cid:44)Hs (θ), 0 1 N N 1 2 N N (cid:16) (cid:17)T (cid:16) (cid:17)T Hc(θ),Hc(θ),··· (cid:44)Hc(θ), Hs(θ),Hs(θ),··· (cid:44)Hs(θ), 0 1 1 2 then (3.1) and (3.2) can be rewritten as follows Hc (θ)=A Ec (θ), Hs (θ)=B Es (θ). (3.3) N N N N N N Itisworthnotingthatthefirstpolynomialis1,andwillbeappliedtothecalculation of density ρ. 6 Remark 3.1. The coefficients A and B in (3.3) can be calculated by using the N N following regular modified cylindrical Bessel function of order n 1 (cid:90) π I (x)= exp(xcosθ)cos(nθ)dθ. (3.4) n π 0 In fact, because of the identities 1 sin(mθ)sin(nθ)=− [cos((m+n)θ)−cos((m−n)θ)], 2 1 cos(mθ)cos(nθ)= [cos((m+n)θ)+cos((m−n)θ)], 2 the Schmidt process is carried out with the aid of I (x), and then the calculation of n A and B is completed by the existing packages. N N 3.2 Hilbert space and orthonormal basis On the interval [0,2π], define Hilbert space H and inner product with respect to the weight M(θ) by (cid:26) (cid:12)(cid:90) 2π dθ (cid:27) H= f(cid:12) f2(θ) <∞ , (3.5) (cid:12) M(θ) 0 (cid:90) 2π dθ (cid:104)f(θ),g(θ)(cid:105) (cid:44) f(θ)g(θ) . (3.6) M(θ) M(θ) 0 Such inner product is symmetric, i.e. (cid:104)f,g(cid:105) =(cid:104)g,f(cid:105) . M(θ) M(θ) Moreover, due to (2.7), the inner product still satisfies (cid:90) 2π (cid:18) f (cid:19) (cid:16) g (cid:17) (cid:104)Q(f),g(cid:105) =−σ M(θ)∂ ∂ dθ =(cid:104)f,Q(g)(cid:105) , (3.7) M(θ) θ M θ M M(θ) 0 and (cid:90) 2π (cid:18) f (cid:19) (cid:104)Q(f),f(cid:105) =−σ M(θ)∂2 dθ (cid:54)0. (3.8) M(θ) θ M 0 Take a basis of H as Pc(θ),··· ,Pc(θ)···, Ps(θ),··· ,Ps(θ)···, and denote 0 N 1 N (cid:16) (cid:17)T Pc(θ),Pc(θ),··· (cid:44)Pc(θ), (3.9a) 0 1 (cid:16) (cid:17)T Ps(θ),Ps(θ),··· (cid:44)Ps(θ). (3.9b) 1 2 Such basis can be generated by the previous orthogonal functions Hc(θ) and Hs(θ) as follows Pc(θ)=Hc(θ)M(θ), Ps(θ)=Hs(θ)M(θ). (3.10) 7 Lemma 2. The functions Pc(θ) and Ps(θ) form an orthonormal basis of H, and satisfy the following properties (cid:104)cos(kθ)M(θ),Pc(θ)(cid:105) =0, k (cid:54)n−1, n M(θ) (cid:104)sin(kθ)M(θ),Ps(θ)(cid:105) =0, k (cid:54)n−1, n M(θ) (3.11) (cid:104)cos(kθ)M(θ),Ps(θ)(cid:105) =0, ∀k,n, n M(θ) (cid:104)sin(kθ)M(θ),Pc(θ)(cid:105) =0, ∀k,n. n M(θ) Because {Pc,Ps} is an orthonormal basis, one has H=span{Pc,Pc,Pc,··· ;Ps,···}, 0 1 2 1 and any function f ∈H may be expressed as follows ∞ ∞ (cid:88) (cid:88) f = fcPc+ fsPs, (3.12) k k k k k=0 k=1 where fc =(cid:104)f,Pc(cid:105) and fs =(cid:104)f,Ps(cid:105). k k k k Take a subspace of H as H =span{Pc,Pc,··· ,Pc;Ps,··· ,Ps}. N 0 1 N 1 N It is obvious that {Pc (θ),Ps (θ)} forms an orthonormal basis of H , where N N N (cid:16) (cid:17)T Pc(θ),Pc(θ),··· ,Pc(θ) (cid:44)Pc (θ), (3.13a) 0 1 N N (cid:16) (cid:17)T Ps(θ),Ps(θ),··· ,Ps(θ) (cid:44)Ps (θ). (3.13b) 1 2 N N For any f ∈H, it is expanded as follows ∞ ∞ f(t,x,θ)=(cid:88)fc(t,x,θ¯)Pc(θ−θ¯)+(cid:88)fs(t,x,θ¯)Ps(θ−θ¯), k k k k k=0 k=1 so one can define projection operator Π [θ¯]:H(cid:55)→H by N N N N Π [θ¯]f =(cid:88)fc(t,x,θ¯)Pc(θ−θ¯)+(cid:88)fs(t,x,θ¯)Ps(θ−θ¯). (3.14) N k k k k k=0 k=1 (cid:16) (cid:17)T Let f˜ = fc,··· ,fc,fs,··· ,fs , the above equation can be rewritten as follows 0 N 1 N Π f =f˜TP (θ−θ¯). (3.15) N N In the case of no confusion, the symbol θ¯will be ignored in the projection operator, that is, Π =Π [θ¯]. N N 8 In the following let us projecting the product of cosθ or sinθ and orthonormal basis. Let (cid:32) (cid:33) Pc (θ) P (θ)= N , (3.16) N Ps (θ) N andtheithcomponentofP isdenotedby(P ) ,equaltoPc(θ),ifi(cid:54)N,otherwise N N i i Ps (θ). i−N Lemma 3 (Projecting product of velocity and basis). The result of the operator Π N acting on the product of velocity and basis is Π [θ¯]cosθP (θ−θ¯)=Jc(θ¯)P (θ−θ¯), (3.17a) N N N Π [θ¯]sinθP (θ−θ¯)=Js(θ¯)P (θ−θ¯), (3.17b) N N N where Jc(θ¯) ∈ R(2N+1)×(2N+1) and Js(θ¯) ∈ R(2N+1)×(2N+1). Both matrices Jc(θ¯) and Js(θ¯) are symmetric and thus can be real diagonalizable. For any eigenvalue λ, |λ|(cid:54)1, and Jc(θ¯),Js(θ¯) has the following form Jc(θ¯)=cosθ¯J −sinθ¯J , Js(θ¯)=sinθ¯J +cosθ¯J , (3.18) 1 2 1 2 where J and J satisfy 1 2 Π [0]cosθP (θ)=J P (θ), Π [0]sinθP (θ)=J P (θ). N N 1 N N N 2 N Proof: According to (3.12) and (3.14), the (i,j)th component of matrix Jc(θ¯) is calculated as follows Jicj(θ¯)=(cid:104)cosθ(PN)i(θ−θ¯),(PN)j(θ−θ¯)(cid:105)M(θ−θ¯). (3.19) Using the definition of inner product (3.6) gives Jc(θ¯)=Jc(θ¯), ij ji so Jc(θ¯) is symmetric and can be written as follows Jc(θ¯)=(cid:82)2πcosθP PT dθ. 0 N NM Because P is an orthonormal basis, (cid:82)2πP PT dθ = I. For any λ ∈ R, and N 0 N NM non-zero vector x∈R2N+1, one has (cid:18)(cid:90) 2π dθ(cid:19) (cid:90) 2π dθ xT(λI−Jc(θ¯))x=xT (λ−cosθ)P PT x= (λ−cosθ)(xTP )2 . N NM N M 0 0 When λ > 1 or λ < −1, λ−cosθ > 0 or λ−cosθ < 0, so that the above is greater than 0 or less than 0. Thus the matrix λI −Jc(θ¯) is positive definite or negative definite, λI−Jc(θ¯) is non singular, so λ>1 or λ<−1 is not eigenvalue. Therefore, |λ|(cid:54)1. For Js(θ¯), the conclusion can be similar proved. 9 Becausecosθ =cosθ¯cos(θ−θ¯)−sinθ¯sin(θ−θ¯),sinθ =sinθ¯cos(θ−θ¯)+cosθ¯sin(θ− θ¯), using the definition of J ,J gives 1 2 Jicj(θ¯)=(cid:104)cosθ(PN)i(θ−θ¯),(PN)j(θ−θ¯)(cid:105)M(θ−θ¯) =cosθ¯((cid:104)cosθ(P ) (θ),(P ) (θ)(cid:105) )−sinθ¯((cid:104)sinθ(P ) (θ),(P ) (θ)(cid:105) ) N i N j M(θ) N i N j M(θ) =cosθ¯(J ) −sinθ¯(J ) . 1 ij 2 ij Thus one gets Jc(θ¯)=cosθ¯J −sinθ¯J . 1 2 Similarly, one has Js(θ¯)=sinθ¯J +cosθ¯J . 1 2 Lemma 4 (Projecting derivatives of basis). The projection of derivatives of basis function P (θ−θ¯) is N Π d(P (θ−θ¯))=D˜P (θ−θ¯)d(θ−θ¯), (3.20) N N N where D˜ ∈R(2N+1)×(2N+1) is a constant matrix, and its (i,j)th element is (cid:90) 2π (P ) (θ) − (P ) (θ)d N j . N i M(θ) 0 Proof: Similar to the calculation of (3.19), the (i,j)th component of matrix D˜ is calculated as follows D˜ij =(cid:104)d(PNi(θ−θ¯)),PNj(θ−θ¯)(cid:105)M(θ−θ¯) (cid:90) 2π d((P ) (θ)) (cid:90) 2π (P ) (θ) = (P ) (θ) N i =− (P ) (θ)d N j , N j M(θ) N i M(θ) 0 0 where the periodicity of the basis functions and M(θ) has been used in the second equal sign while the integration by parts is used in the third equal sign. Lemma 5 (Projecting collision term). The result of the operator Π acting on the N collision term is Π Q(P (θ−θ¯))=Q P (θ−θ¯), (3.21) N N N N where Q ∈ R(2N+1)×(2N+1) is a symmetric and negative semidefinite matrix, and N its first column and first row elements are zeros. Proof: Using the property of inner product (3.7), one has (QN)ij =(cid:104)Q((PN)i(θ−θ¯)),(PN)j(θ−θ¯)(cid:105)M(θ−θ¯) =(cid:104)(PN)i(θ−θ¯),Q((PN)j(θ−θ¯))(cid:105)M(θ−θ¯) =(cid:104)Q((PN)j(θ−θ¯)),(PN)i(θ−θ¯)(cid:105)M(θ−θ¯) =(QN)ji. 10