A Study of Hyperbolicity of Kinetic Stochastic Galerkin System for the Isentropic Euler Equations with Uncertainty*

Shi JIN Ruiwen SHU

(Dedicated to Professor Andrew J.Majda on the occasion of his 70th birthday)

Abstract The authors study the fluid dynamic behavior of the stochastic Galerkin (SG for short)approximation to the kinetic Fokker-Planck equation with random uncertainty.While the SG system at the kinetic level is hyperbolic,its fluid dynamic limit,as the Knudsen number goes to zero and the underlying kinetic equation approaches to the uncertain isentropic Euler equations,is not necessarily hyperbolic,as will be shown in the case study fashion for various orders of the SG approximations.

Keywords Hyperbolic equations,Uncertainty quantification,Stochastic Galerkin methods

1 Introduction

Hyperbolic conservation laws are one of the classical nonlinear partial differential equations with important applications from gas dynamics,shallow water,combustion,to magnetohydrodynamics.A general system of conservation laws reads

where t ∈R≥0is time,and x ∈Ω ?R is space,assumed to be one-dimensional for simplicity.U(t,x)= (U1(t,x),··· ,Un(t,x))Tis the vector of conserved quantities,and F(U)=(F1(U),··· ,Fn(U))is the flux function,assumed to be at least C1.The system is called hyperbolic if the Jacobian matrix ?UF is real-diagonalizable.For most models from physics in conservation forms,in particular,the compressible Euler equations in gas dynamics,the hyperbolicity condition is satisfied.

Consider the isentropic Euler equations

Here ρ ＞0 is the density of gas,and m is the momentum.The Jacobian of ?UF is given by

Hydrodynamic equations have two components:Conservations and equations of state (or constitutive relations).While the conservation properties are basic physical properties typically satisfied for an ensemble of large number of particles,the equation of state or constitutive relations are usually empirical,thus inevitably contains uncertainty.Uncertainties also arise from initial and/or boundary data,and/or forcing or source terms,due to measurement or modeling errors.By quantifying how these uncertainties affect the behavior of the solution,one can make reliable predictions of the physical quantities modeled by the conservation laws and further calibrate or improve the models.To model the uncertainty,we introduce a random variable z lying in a random space Iz?Rdwith a probability distribution π(z)dz,and any uncertain quantity is described by its dependence on z.For example,uncertain initial data is modeled by Uin= Uin(x,z).As a result,the solution U also depends on z:U = U(t,x,z).The solution U(t,x,z)still satisfies the same equations (1.1),but depends on an extra random parameter z.

In the last few decades,many methods have been developed for uncertainty quantification(UQ for short)(see[8-9,17,28-29]),including the Monte Carlo methods,stochastic collocation(SC for short)methods and stochastic Galerkin (SG for short)methods.The Monte Carlo methods (see [20])solve the deterministic problem on random sample points,and then obtain statistical information by averaging.Stochastic collocation methods (see [1,3,21,30])solve the deterministic problem on some well-chosen sample points (for example,quadrature points,or points chosen by some optimization procedure),and then compute statistical quantities by interpolation.Stochastic Galerkin methods (see [2-3,31])use a finite (K-)dimensional orthonormal basis in the random space Iz(with respect to π(z)dz),expand any function in z in the L2space as a linear combination of this basis.After conducting the Galerkin projection one obtains a deterministic system of the K coefficients in the expansion.A popular choice of basis is the generalized polynomial chaos (gPC for short)basis (see [31]),which are the orthonormal polynomials with respect to π(z)dz.

The main advantage of the Monte Carlo methods is that it maintains a convergence orderfor any dimensional random spaces,thus it is efficient even if the dimension d is large.SC and SG methods can achieve high convergence order if the solution is sufficiently smooth in the random space,but they suffer from the ‘curse of dimensionality’ if d becomes large.

The Monte Carlo and SC methods are non-intrusive methods:One only needs to run the deterministic solver many times to obtain statistical information of the random solution.However,for the intrusive gPC-SG methods,the system of the gPC coefficients is a coupled system,which may be significantly different from the deterministic equations.One usually needs to design new numerical methods to solve it.

When applied to system of hyperbolic conservation laws (1.1),the gPC-SG method gives a system of conservation laws for the gPC coefficients which may lose the hyperbolicity,unless for special classes of systems which include:

· Symmetric conservation laws,i.e.,?UF is symmetric.This includes scalar conservation laws.See [11] for an example.

· Linear conservation laws,i.e.,F is linear.

One example of loss of hyperbolicity is given in [5] for the Euler equations.When hyperbolicity of the gPC-SG method is lost,the initial value problem of the gPC-SG system becomes ill-posed.

One possible fix of the loss of hyperbolicity of the gPC-SG method is proposed by [24].By using entropy variables,they write (1.1)into a symmetric form,and this symmetry can be maintained by the gPC-SG method.However,this method requires the transformation between the conservative variables U and the entropy variables in each time step and every grid point,which is very costly in general.For other efforts in this direction see [7,23,26].

It is well known that many hyperbolic conservation laws can be obtained as the hydrodynamic limit of kinetic equations.For example,the Euler equations are the hydrodynamic limit of the Boltzmann equation (with the hyperbolic scaling)(see [4]),and the isentropic Euler equations (1.2)is the hydrodynamic limit of the kinetic Fokker-Planck equation (see [27])

where v ∈R is the velocity variable,f =f(t,x,v)is the particle distribution function,and Lfis the kinetic Fokker-Planck operator

Here Kn is the Knudsen number,the dimensionless particle mean free path.As Kn →0,the hydrodynamic limit (1.2)is (formally)obtained.

Recently there has been a rapid progress on the gPC-SG methods for multiscale kinetic equations with uncertainty (see [10,13-15,32]).Many gPC-SG methods have the stochastic asymptotic-preserving(s-AP for short)property(see[15]),which means that the method works uniformly from the kinetic regime (Kn = O(1))to the fluid regime (Kn ＜＜1),with a fixed number K of basis functions.Also,the random space regularity of the solutions for various kinetic equations have been proved,in the linear cases (see [12-13,18]),and nonlinear ones(see [16,19,25]).Consequently one can even establish the spectral accuracy,uniform in Kn,of the gPC-SG methods.In the case of nonlinear kinetic equations,such results were obtained for solutions near the global Maxwellian.Note that at the kinetic level,the equation is a scalar equation,and moreover,with a linear convection,thus the gPC-SG approximations naturally preserve the hyperbolicity of the kinetic equations.Therefore,one may be curious whether the hyperbolicity of the kinetic gPC-SG system for a kinetic equation will survive in the fluid dynamic limit.If it does,then this would provide a vehicle to derive stable gPC-SG approximations for the compressible Euler equations,as in their deterministic counterparts,which is known as the kinetic schemes for the Euler equations (see [22]).

In this work,we give a counter-argument to this.For the kinetic Fokker-Planck equation(1.4),we show that the gPC-SG system for (1.4),after taking a formal hydrodynamic limit by sending Kn →0,does not necessarily give a hyperbolic gPC-SG approximation for (1.2).To be precise,our numerical results suggest that the gPC-SG method obtained from a kinetic approach may fail to be hyperbolic even for arbitrarily small randomness,if K = 4.However,for K =2,3,we prove that hyperbolicity holds under reasonable assumptions.

The paper is organized as follows:In Section 2 we introduce the gPC-SG method for the kinetic Fokker-Planck equation (1.4)with uncertainty,and formally derive its hydrodynamic limit as the gPC-SG method for the isentropic Euler equations (1.2); in Section 3 we show that the Jacobian matrix of the gPC-SG system of(1.2),if contains imaginary eigenvalues,the imaginary part will be small if the uncertainty of ρ and m is small; in Section 4 we prove the hyperbolicity of this gPC-SG system for K = 2,3 under reasonable assumptions,and provide numerical evidence suggesting that hyperbolicity may fail for arbitrarily small randomness when K =4.The paper is concluded in Section 5.

2 The gPC-SG Method

In this section we aim to derive a gPC-SG method for the isentropic Euler equations (1.2)with uncertainty as the hydrodynamic limit of a gPC-SG method for the kinetic Fokker-Planck equation (1.4)with uncertainty.

We start by deriving formally the hydrodynamic limit of the deterministic kinetic Fokker-Planck equation (1.4).The Fokker-Planck operator (1.5)conserves mass and momentum:

The local equilibrium of Lf(f),i.e.,the f with Lf(f)=0,is given by the local Maxwellian

with ρ(t,x)and u(t,x)being the local density and bulk velocity.As Kn →0 in (1.4),one can see (formally)that f is at this local equilibrium.Then by taking moments of (1.4)against 1,v and using (2.2),one obtains the hydrodynamic limit (1.2)with m=ρu.

2.1 The gPC-SG method for (1.4)

For simplicity,we assume the random space is one-dimensional.The gPC basis associated to π(z)dz is denoted byis a polynomial of degree k - 1,satisfying the orthonormal property

We expand the solution f(t,x,v,z)to (1.4)into gPC series and truncate:

Here fk(t,x,v)are the gPC coefficients,which no longer depend on the random variable z.

If one directly substitutes (2.4)into (1.4)and conduct Galerkin projection,then the termin the Fokker-Planck operator will becomedz,with ρK,mKdefined similarly to (2.4).This term is not easy to compute numerically from the gPC coefficients f1,··· ,fK.Instead,by conducting the gPC-SG method on the identity

we get an approximation

where

is the vector of the gPC coefficients,and A(ρ)is the multiplication matrix given by

Notice that A(ρ)is symmetric,and positive-definite iffor all z.In this case,

Similarly,we obtain

Using these two approximations,we obtain the gPC system

Remark 2.1There is more than one way to approximate the gPC coefficients of the termf,for example,by A(ρ)-1A(m)However,(2.12)seems to be the only consistent approximation which is easy to compute and conserves momentum (see next section for the conservation properties).

2.2 Properties of the gPC system

(2.12)has the same conservation properties.

Lemma 2.1(2.12)conserves mass and momentum:

ProofIt is clear thatTo see the second equality,

Notice that if we use A(ρ)-1A(m)instead of A(m)A(ρ)-1then the momentum conservation does not hold.

Then,similar to [14],we have the local equilibrium ofgiven by the following.

Lemma 2.2Assume ρK(z)＞0 for all z.Thenimplies

where exp(-A?v)is defined as the solution at s=1 to the system

ProofWe first show that exp(-A?v)is well-defined,i.e.,the linear system of conservation laws (2.17)is hyperbolic.We need to show that the matrix A is real-diagonalizable.To see this,notice that

Then,to see the local equilibrium (2.16),we use the conjugate structure of →L:

Remark 2.2A rigorous proof of (2.19)can be obtained in a similar way to [14].We omit the details.

Using (2.16),one obtains

where we first integrate by parts,and use the fact that only the terms with j = 0,2 do not vanish.Taking moments of(2.11)and using this,one obtains the limiting fluid dynamic system

(2.22)is a system of conservation laws for the 2K functionsandThe Jacobian of the flux function of (2.22)is given by the block matrix

where I stands for the identity matrix,

and

3 Analysis of Hyperbolicity of (2.22)for Small Randomness

In this section we study the eigenvalues of J,in the case when the randomness is small.

Theorem 3.1The imaginary part of any eigenvalue of J is at most O(∈2).

To prove the theorem,we start with the following lemma.

Lemma 3.1For givenwe denote the corresponding J asThen one has

where eig means the eigenvalues of matrix,and the last +α means adding α to each eigenvalue.

ProofExplicit calculation shows that

where B(1),B(2)are those inThus

Without loss of generality,we may assume ρ1=1.To analyze the imaginary parts of eig(J),one can replaceto make m1=0,in view of Lemma 3.1.Since m1=O(1),this replacement does not affect any smallness condition.From now on we will assume

Then we define the following matrix:

Also define

and its inverse is given by

Explicit computation shows

Lemma 3.2Assume (3.4).Then

ProofIt suffices to check B(1)=-A(m)2+O(∈2)and B(2)=2A(m)+O(∈2).

The equation for B(1)is clear since both B(1)and -A(m)2are of order O(∈2).

To check the expansion for B(2),notice that

To analyze the size of imaginary part of eigenvalues of J,we need the following lemma,which is a consequence of the Gershgorin Circle Theorem (see [6]).

Lemma 3.3Suppose A,B are real matrices.Suppose A can be real-diagonalized by PAP-1= D = diag{λ1,··· ,λK} with max(‖P‖,‖P-1‖)≤C1and ‖B - A‖ ≤C2.Here ‖A‖ =Then

(i)each eigenvalue η of B satisfiesfor some i.In particular,the imaginary part of η is at most C;

(ii)if |λi-λj|＞2C for each i ＜j,then the eigenvalues {ηk} of B can be arranged so that|ηk-λk|≤C,k =1,··· ,K,and B is real-diagonalizable.

ProofFirst,

Since‖P(B-A)P-1‖≤‖P‖‖B-A‖‖P-1‖≤C,it follows from the Gershgorin Circle Theorem that each eigenvalue η of B satisfies |η-λi|≤C for some i.This proves (i).

To prove(ii),let ηk(t)be the k-th eigenvalue of(1-t)A+tB,for 0 ≤t ≤1.Then each ηk(t)is a complex-valued continuous function,with ηk(0)=λk.Also,for each fixed k,t,one has the estimate |ηk(t)-λi|≤C for some i.Since |λi-λj|＞2C for each i ＜j,the only possibility is that |ηk(t)-λk| ≤C for 0 ≤t ≤1.This shows that |Re(ηk(1))-λk| ≤C,and thus {ηk(1)},the eigenvalues of B,have distinct real parts.Since B is a real matrix,this implies that B is real-diagonalizable.

Finally we prove Theorem 3.1.

Proof of Theorem 3.1The matrixcan be diagonalized by

where P1is an orthogonal matrix,P2is defined in (3.6),and λ1,··· ,λKstand for the real eigenvalues of the symmetric matrix A(m).This is because

is symmetric(see(3.8)),and has eigenvalues{λk±1,k=1,··· ,K}.It is clear that P1,P2,,are of order O(1).Therefore the conclusion follows from Lemma 3.2 and Lemma 3.3(i).

4 Analysis of Hyperbolicity of (2.22)for Specific Values of K

In this section we analyze the hyperbolicity of(2.22)for specific values of K.We will show:

· For K =2,(2.22)is hyperbolic iffor all z ∈Iz(Theorem 4.1).

· For K = 3,(2.22)may fail to be hyperbolic if one only assumes ρK(z)＞0 (Subsection 4.2.1).However,(2.22)is hyperbolic if ρ1= O(1),m1= O(1)and ρ2,ρ3,m2,m3are small enough (Theorem 4.2),for some special cases of π(z).

· For K ≥4,we provide numerical evidence showing that (2.22)may fail to be hyperbolic even if ρ1=O(1),m1=O(1)and ρk,mk,k =2,··· ,K are small (Subsection 4.3).

4.1 The case K = 2

Lemma 4.1Assume K =2.Define the 2*2 matrix Φij=φj(zi),where z1,z2∈Izare the two-point Gauss quadrature points (with respect to the probability measure π(z)dz).Then for any=(f1,f2)T,we have

ProofWe first show that the matrix Ψ defined by Ψij= φi(zj)wj,where w1,w2are the quadrature weights corresponding to z1,z2,is the inverse of Φ.In fact,

where the second equality is because φiφkis a polynomial of degree at most 2,and the quadrature is exact for polynomials of degree at most 3.Similarly,

since φiφjφkis a polynomial of degree at most 3.

Then,denoting Sias the matrix given by (Si)jk=Sijk,we have

where in the third equality we used ΦΨ=I.Thus

and the conclusion follows.

Then it follows easily.

Theorem 4.1Assume K =2.If ρK(z)＞0 for all z ∈Iz,then (2.22)is hyperbolic.

ProofWe use the same notation as Lemma 4.1.Multiplying by Φ on each equation of(2.22)and using Lemma 4.1 gives

Remark 4.1In fact,this proof shows that for K =2 the gPC system (2.22)is equivalent to a stochastic collocation method at z1,z2.

4.2 The case K = 3

In this subsection we assume K =3 and prove the hyperbolicity of(2.22)for small randomness.Same as the previous subsection,we will assume (3.4)without loss of generality.

Lemma 4.2Assume (3.4)and |m2|+|m3| = 1.Then there exists c1＞0 such that any two eigenvalues λi,λjof the symmetric matrix A(m)satisfy |λi-λj|≥c1.

Proof for Special DistributionsIt is clear that all λiare at most O(1).Therefore,it suffices to show that the discriminant of det(λI -A(m))(as a cubic polynomial of λ)is away from zero.Let us write x=m2,y =m3,a=S223,b=S333for clarity,

Its discriminant is given by the degree-6 homogeneous polynomialHere G(x)is a degree three polynomial,given by the following in the special cases:

· When Iz=[-1,1]with the uniform distribution,φkare the normalized Legendre polynomials,andThen

·When Iz=R with the Gaussian distribution,φkare the normalized Hermite polynomials,andThen

One can easily verify that for both cases G(x)is positive for x ＞0.Therefore disc is always positive,ifSince the set |x|+|y| = 1 is compact and disc is continuous with respect to x and y,disc has a positive lower bound.

Theorem 4.2Assume (3.4)and |ρ2|+|ρ3|+|m2|+|m3| ≤c with c ＞0.If c is small enough,then J defined in (2.23)is real-diagonalizable,i.e.,(2.22)is hyperbolic.

ProofWe denote δ =|m2|+|m3|.Then we claim

where B(1),B(2)are the block matrices appeared in (2.23).The equation for B(1)is clear since both B(1)and -A(m)2are of order O(δ2).To check the equation for B(2),

where ρ′= (0,ρ2,ρ3)T,and we used A(ρ′)= O(c),= O(δ).We also used the fact that c is small enough so that

Now we claim that if c is small enough,thenJ and J satisfy the assumptions of Lemma 3.3(ii).In fact,we already showed in Theorem 3.1 that the transformation matrix P = P1P2of(defined in the proof of Theorem 3.1)is C1= O(1)(following the notation in Lemma 3.3).Then we have ‖J -‖ = C2= O(cδ).Therefore C = C21C2= O(cδ).We also have|μi-μj| ≥c1δ.Therefore,if we choose c small enough such that C = O(cδ)＜,then one has |μi-μj|≥2C.Then the conclusion that J is real-diagonalizable follows from Lemma 3.3(ii).

4.2.1 Loss of hyperbolicity forK = 3 with large uncertainty

We show by counterexample that for K = 3,even if ρK(z)＞0 for all z ∈Iz,(2.22)may still fail to be hyperbolic if the uncertainty of m is large.

We take Iz=[-1,1]with the uniform distribution(so that φk(z)are the normalized Legendre polynomials).We take →ρ=(1,-0.1076151078,-0.01061304986)T,→m=(0,4.678859141,3.49 8325096)T.In this case one clearly has ρK(z)＞0.Numerical result shows that J has the complex eigenvalues -2.172329053±0.009698318358i.

4.3 The case K ≥4

For K ≥ 4,we try to find →ρ,→m with (3.4)and small uncertainty (in the sense thatis small),such that J(→ρ,→m)has non-real eigenvalues,and we want the imaginary part of non-real eigenvalues as large as possible.We take Iz= [-1,1] with the uniform distribution (so that φk(z)are the normalized Legendre polynomials).The algorithm for finding such →ρ,→m and the numerical parameters are described in Appendix.

We take K = 4,5,6 and various values of size of randomness ∈.For each K,we plot the maximal imaginary part of eigenvalues of J found by the algorithm (the variable maximag in the algorithm)versus ∈,see Figure 1.One can clearly see that for each K,maximag is proportional to ∈3.In particular,this result strongly suggests that there are examples for which J has non-real eigenvalues for arbitrarily small size of randomness.Also,this result suggests that Theorem 3.1 may not be sharp,the sharp estimate could be O(∈3).

Figure1 Maximal imaginary part of eigenvalues of J (found by Algorithm 1)versus ∈.From top to bottom:K = 4,5,6.Horizontal axis:The size of randomness ∈.Blue asterisks:Maximal imaginary part of eigenvalues of J (found by Algorithm 1).Red line:slope=3.

5 Conclusion

In this paper,we study the fluid dynamic behavior of the stochastic Galerkin system to the kinetic Fokker-Planck (FP for short)equation with random uncertainty.As the Knudsen number goes to zero,the FP system approaches to the isentropic Euler system with uncertainty.For various orders of the Stochastic Galerkin (SG for short)system for the FP,we show that their fluid dynamic limits,which give rise to SG systems to the Euler equations,are not necessarily hyperbolic.Thus one cannot expect to derive kinetic SG system to the isentropic Euler equations with uncertainty from the kinetic approximation,in contrast to the deterministic case where kinetic approximations provide a robust mechanism to derive kinetic schemes for the Euler equations.

Our analysis was carried out in the case study fashion,for various orders of the SG approximation.So far there has been no general theory for general order of the SG approximation.It will be interesting to develop a more general theory for the problem.It will also be interesting to study the problem for the full Euler equations as the fluid dynamic limit of the Boltzmann equation.

6 Appendix:The Algorithm for Finding J with Non-real Eigenvalues

We use the following algorithm to find J with non-real eigenvalues:For a randomly chosenwe try to adjustto decrease the minimum distance of two distinct real eigenvalues of J as much as possible.This procedure drives J to have multiple eigenvalues,and one expects a small perturbation may make J have non-real eigenvalues,in case there are no obstruction for the eigenvalues of being non-real.Then we try to adjustto increase the imaginary parts of non-real eigenvalues.

Here we provide the code for the algorithm used in Subsection 4.3 to find J with non-real eigenvalues.It looks forsuch that the size of randomness,measured by

is equal to a given ∈＞0.We need three numerical parameters ntrial,niter1and niter2,whose meanings are explained as follows:For each randomly chosenwe try niter1times to adjust it in order to decrease the minimum distance of two distinct real eigenvalues of J.In case we fnidsuch that J has non-real eigenvalues,we then try niter2times to adjust it in order to increase the imaginary parts of non-real eigenvalues.The whole procedure is repeated ntrialtimes and the largest imaginary part of non-real eigenvalues (the variable maximag)is reported.In all the numerical results,we take the numerical parameters ntrial= 100 and niter1=niter2=1000.

AcknowledgementThe authors thank Prof.Pierre Degond for stimulating discussions and encouragement for this work.

Algorithm 1 Look for →ρ,→m such that J has large non-real eigenvalues Require:∈＞0,ntrial,niter1,niter2 ∈Z＞0 Ensure:Print →ρmax,→mmax,maximag such that ρ1 = 1,m1 = 0,K∑(|ρk| + |mk|)= ∈,and eig(J(→ρmax,→mmax))has a complex eigenvalue with imaginary part maximag ＞0,or print ‘Did not found’.1:maximag ←0 2:for itrial =1:ntrial do 3:ρ1 ←1,m1 ←0 4:Randomly generate ρk,mk,k =2,··· ,K such that k=2 K∑(|ρk|+|mk|)=∈5:eigJ ←eig(J(→ρ,→m))6:ind ←min{|λ1-λ2|:λ1,λ2 ∈eigJ,λ1 /=λ2}7:for iiter1 =1:niter1 do 8:if not isreal(eigJ)then 9:break 10:end if 11:Perturb ρk,mk,k = 2,··· ,K into ρ1k,m1k randomly,with size of perturbation ind,such that they do not change sign,and k=2 K∑(|ρk|+|mk|)=∈is kept the same.12:eigJ1 ←eig(J(→ρ1,→m1))13:ind1 ←min{|λ1-λ2|:λ1,λ2 ∈eigJ1,λ1 /=λ2}14:if not isreal(eigJ1)or ind1 ＜ind then 15:→ρ ←→ρ1,→m ←→m1,eigJ ←eigJ1,ind ←ind1 16:end if 17:end for 18:if isreal(eigJ)then 19:continue 20:end if 21:ind ←max{Im(λ):λ ∈eigJ}22:for iiter2 =1:niter2 do 23:Perturb ρk,mk,k = 2,··· ,K into ρ1k,m1k randomly,with size of perturbation ind,such that they do not change sign,and k=2 K∑(|ρk|+|mk|)=∈is kept the same.24:eigJ1 ←eig(J(→ρ1,→m1))25:if isreal(eigJ1)then 26:ind1 ←0 27:else 28:ind1 ←max{Im(λ):λ ∈eigJ1}29:end if 30:if ind1 ＞ind then 31:→ρ ←→ρ1,→m ←→m1,eigJ ←eigJ1,ind ←ind1 32:end if 33:end for 34:if ind ＞maximag then 35:→ρmax ←→ρ,→mmax ←→m,maximag ←ind 36:end if 37:end for 38:if maximag ＞0 then 39:print →ρmax,→mmax,maximag 40:else 41:print ‘Did not found’42:end if k=2