the l1-error estimates for a hamiltonian-preserving scheme for the ...
THE L1 -ERROR ESTIMATES FOR A HAMILTONIAN-PRESERVING SCHEME FOR THE LIOUVILLE EQUATION WITH PIECEWISE CONSTANT POTENTIALS ∗ XIN WEN
† AND
SHI JIN
‡
Abstract. We study the l1 -error of a Hamiltonian-preserving scheme, developed in [11], for the Liouville equation with a piecewise constant potential in one space dimension. This problem has important applications in computations of the semiclassical limit of the linear Schr¨ odinger equation through barriers, and of the high frequency waves through interfaces. We use the l1 -error estimates established in [30, 28] for the immersed interface upwind scheme to the linear advection equations with piecewise constant coefficients. We prove that the scheme with the Dirichlet incoming boundary conditions is l1 -convergent for a class of bounded initial data, and derive the one-halfth order l1 -error bounds with explicit coefficients. The initial conditions can be satisfied by applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data, which arises in the semiclassical limit of the linear Schr¨ odinger equation. Key words. Liouville equations, Hamiltonian preserving schemes, piecewise constant potentials, error estimate, half order error bound, semiclassical limit AMS subject classifications. 65M06, 65M12, 65M25, 35L45, 70H99
1. Introduction. In [11], we constructed a class of numerical schemes for the d-dimensional Liouville equation in classical mechanics: ft + v · ∇x f − ∇x V · ∇v f = 0 ,
t > 0,
x, v ∈ Rd ,
(1.1)
where f (t, x, v) is the density distribution of a classical particle at position x, time t and traveling with velocity v. V (x) is the potential. The main interest is in the case of a discontinuous potential V (x), corresponding to a potential barrier. When V is discontinuous, the Liouville equation (1.1) is a linear hyperbolic equation with a measure-valued coefficient. Such a problem cannot be understood mathematically using the renormalized solution by DiPerna and Lions for linear advection equations with discontinuous coefficients [5] (see also [2]). Our approach in [11, 12] to such problems was to provide an interface condition to couple the Liouville equation (1.1) on both sides of the barrier or interface. The interface condition accounts for particle or wave transmission and reflection. Based on this notion of the solution, the socalled Hamiltonian-preserving schemes were constructed in [11, 12], which build this interface condition into the numerical flux. Schemes so constructed provide solutions that are physically relevant for particle or wave reflection and transmission through the barriers or interfaces. The Liouville equation is the phase space representation of Newton’s second law: dx = v, dt
dv = −∇x V , dt
∗ Research supported in part by the Knowledge Innovation Project of the Chinese Academy of Sciences grants K5501312S1, K5502212F1, K7290312G7 and K7502712F7, NSFC grant 10601062, NSF grant DMS-0608720 and NSAF grant 10676017. † Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, P.O. Box 2719, Beijing 100190, China. Email address: [email protected]. ‡ Department of Mathematics, University of Wisconsin, Madison, WI 53706, USA. Email address: [email protected].
1
2
X. WEN AND S. JIN
which is a Hamiltonian system with the Hamiltonian H=
1 2 |v| + V (x). 2
It is known from classical mechanics that the Hamiltonian remains constant across a potential barrier. This is one of the main ingredients in the Hamiltonian-preserving schemes developed in [11, 12]. The two schemes developed in [11]–one based on a finite difference formulation (called Scheme I) and the other on a finite volume formulation (called Scheme II) were proved, in one space dimension with a piecewise constant potential, to be positive, and l1 and l∞ -stable for suitable initial value problems and under a hyperbolic CFL condition (see also [29]). The more sophisticated issue of the l1 -stability of Scheme I was studied in [29]. We proved in [29] that, in the case of a step function potential, Scheme I with the homogeneous Dirichlet incoming boundary conditions is l1 -stable under a certain condition on the initial data. In particular, we showed that such an initial condition is satisfied when choosing the numerical initial data as the cell averages of a bounded analytical initial data. Since the exact solution, in the aforementioned notion using the interface condition, to the same problem is l1 -contracting, an l1 -convergent scheme should be l1 -stable. Thus the results established in [29] implies that Scheme I meets the necessary condition for the l1 -convergence in the case of the homogeneous Dirichlet boundary conditions and bounded initial data. In this paper we will study the l1 -convergence of Scheme I in the case of a step function potential and the Dirichlet incoming boundary conditions. We will prove that in this case Scheme I is l1 -convergent for a class of bounded initial data. The initial conditions can be satisfied when applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data arising in the semiclassical limit of the linear Schr¨odinger equation. The Liouville equation with a step function (or more generally piecewise constant) potential belongs to hyperbolic equations with singular (discontinuous or measurevalued) coefficients. For the discontinuous coefficient case, the convergence of numerical schemes were widely studied [6, 9, 14, 15, 4, 17, 25, 19, 21, 7, 26, 27, 16, 13, 1, 20]. However, in most cases the convergence rate estimates for numerical schemes were not studied. For the measure-valued coefficient case, recently we have obtained the halfth order l1 -error estimates for the immersed interface upwind scheme which builds the interface condition into the numerical flux, see [30, 28]. In this paper, we extend this result to the more interesting case of the Liouville equation with a step function potential. Since the Liouville equation is defined in the phase space, extra efforts are needed to deal with the discretization in the velocity space. To derive the desired error estimates, we split the computational domain into several parts. In each subdomain, we will introduce a number of linear advection equations with step function coefficients. The immersed interface upwind schemes for these linear advection equations yield the same numerical solution as Scheme I for the Liouville equation with the step function potential. Then the l1 -error estimates for Scheme I can be achieved by applying the l1 -error estimates for the immersed interface upwind schemes and estimating the l1 -errors between the exact solutions of the linear advection equations with step function coefficients and that of the Liouville equation with the step function potential. The first part of the error estimates are obtained by applying the results established in [30, 28] with the aid of the condition on the initial data. The second part of the errors are analyzed also utilizing the condition on the
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
3
initial data. With this approach, we can derive the halfth order l1 -error bound with explicit coefficients for Scheme I. Our proof is based on the step function potential. This result clearly generates to piecewise constant potentials in a straightforward but tedious way. The same approach should also be applicable to analyze the l1 -error estimates for the finite volume Scheme II developed in [11]. But we will not pursue this in this paper. This paper is organized as follows. In Section 2 we review the Hamiltonianpreserving scheme proposed in [11] for the Liouville equation with a discontinuous potential in one space dimension. In Section 3 we recall the l1 -error estimates for the immersed interface upwind scheme established in [30, 28] to the linear advection equation with a step function wave speed. In Section 4 we setup the problem and state the main Theorem of this paper which is proved in Section 5 using the results presented in Section 3. We conclude the paper in Section 6. 2. A Hamiltonian-preserving scheme. In this Section we review the Hamiltonianpreserving scheme proposed in [11] to the Liouville equation in one space dimension ft + ξfx − Vx fξ = 0
(2.1)
with a discontinuous potential V (x). Consider a uniform mesh with grid points at xi+ 12 , i ∈ Z in the x-direction and ξj+ 12 , j ∈ Z in the ξ-direction. The cells are centered at (xi , ξj ), (i, j) ∈ Z2 with xi = 21 (xi+ 12 + xi− 12 ) and ξj = 12 (ξj+ 12 + ξj− 12 ). The mesh size are denoted by ∆x = xi+ 12 − xi− 12 , ∆ξ = ξj+ 21 − ξj− 12 . We also assume a uniform time step ∆t and the discrete times are given by tn = n∆t, n ∈ N ∪ {0}. We assume that the computation is performed in a bounded rectangular domain {(x, y)|x 21 ≤ x ≤ xN + 12 , ξ 21 ≤ ξ ≤ ξM + 12 }.
(2.2)
Let the cell averages of f be 1 fij = ∆x∆ξ
Z
xi+ 1
Z
2
ξj+ 1 2
xi− 1
f (x, ξ, t)dξdx.
ξj− 1
2
2
The 1-d average quantity fi+1/2,j is defined as fi+1/2,j =
1 ∆ξ
Z
ξj+ 1 2
ξj− 1
f (xi+1/2 , ξ, t)dξ .
2
fi,j+1/2 is defined similarly. In classical mechanics, a particle will either cross a potential barrier with a changing momentum, or be reflected, depending on its momentum and the strength of the potential barrier. Figure 2.1 shows the typical situations when a particle moves from left to right at a potential barrier. If the potential is reduced, the particle will cross it with an accelerated speed. If the potential is increased, the particle is either reflected if its momentum is not big enough to overcome the barrier or transmitted otherwise with a reduced velocity. The Hamiltonian H = 12 ξ 2 + V is preserved across the potential barrier: 1 + 2 1 (ξ ) + V + = (ξ − )2 + V − , 2 2
(2.3)
4
X. WEN AND S. JIN
2
+
−
[ξ −2(V
ξ
1/2
− V )]
1)
2)
ξ
[ξ2+2(V− − V+)]1/2
−ξ
3)
ξ
V−
V−
V+
V+
Figure 2.1 Transmission and reflection of a particle at a potential barrier. where the superscripts ± indicate the right and left limits of the quantity at the potential barrier. This property was used in [11] to provide the interface condition for (2.1) at the barrier: f (t, x+ , ξ + ) = f (t, x− , ξ − ) ±
±
±
for transmission ±
f (t, x , ξ ) = f (t, x , −ξ )
(2.4)
for reflection
(2.5)
where ξ + and ξ − are related by the constant Hamiltonian condition (2.3) in the case of transmission. With such an interface condition, we established the well-posedness of the initial value problem to (1.1) with a piecewise constant potential in [12]. The main ingredient in the Hamiltonian-preserving schemes developed in [11], like the earlier work for shallow-water equations [22], was to build into the numerical flux the interface conditions (2.4), (2.5) at the barrier. We now present the first Hamiltonian-preserving scheme, called Scheme I in [11]. Assume that the discontinuous points of the potential V are located at the grid − + points. Let the left and right limits of V at point xi+1/2 be Vi+ respectively. 1 and V i+ 1 2
2
− + Note that if V is continuous at xj+1/2 , then Vi+ . We approximate V by a 1 = V i+ 12 2 piecewise linear function
V (x) ≈
+ Vi−1/2
+
− + Vi+1/2 − Vi−1/2
∆x
(x − xi−1/2 ) .
The flux-splitting, semidiscrete scheme (with time continuous) reads ∂t fij + ξj
− + fi+ − fi− 1 1 ,j ,j 2
2
∆x
−
− + Vi+ 1 − V i− 12 fi,j+ 12 − fi,j− 21 2
∆x
∆ξ
= 0,
where the numerical fluxes fi,j+ 12 are defined using the upwind discretization. Since the characteristics of the Liouville equation may be different on the two sides of a barrier, the corresponding numerical fluxes should also be different. The essential + − , fi− at each cell part of the algorithm is to define the split numerical fluxes fi+ 1 1 2 ,j 2 ,j interface. (2.4) will be used to define these fluxes. Assume V is discontinuous at xi+1/2 . Consider the case ξj > 0. Using upwind − scheme, fi+ = fij . However, 1 ,j 2
+ − + − fi+1/2,j ≡ f (x+ i+1/2 , ξ ) = f (xi+1/2 , ξ )
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
5
where ξ − is obtained from ξ + = ξj from (2.3). Since ξ − may not be a grid point, we have to define it approximately. The first approach is to locate the two cell centers that bound this velocity, then use a linear interpolation to evaluate the needed numerical flux at ξ − . The case of ξj < 0 is treated similarly. The algorithm to generate the numerical flux is given in [11]. Here we present the simplified algorithm for the case + − being discussed in this paper. Vi+ 1 > V i+ 12 2 Algorithm I • ξj > 0 − = fij , fi+ 1 2 ,j r ³ ´ − + o if ξj > 2 Vi+ , 1 − V 1 i+ 2 2 r ³ ´ + − ξ − = ξj2 + 2 Vi+ 1 − V 1 i+ 2
2
if ξk ≤ ξ − < ξk+1 for some k − − −ξk + then fi+ = ξk+1∆ξ−ξ fik + ξ ∆ξ fi,k+1 1 ,j o else
2
+ fi+ = fi+1,k where ξk = −ξj 1 ,j 2
o end • ξj < 0 + fi+ = fi+1,j , 1 2 ,j r ³ ´ − + ξ + = − ξj2 + 2 Vi+ 1 − V 1 i+ 2
2
if ξk ≤ ξ + < ξk+1 for some k + −ξ + −ξk − fi+1,k + ξ ∆ξ fi+1,k+1 then fi+ = ξk+1 1 ∆ξ 2 ,j After the spatial discretization is specified, one can use any time discretization for the time derivative. In [11] we proved that, when the first order upwind scheme is used spatially, and the forward Euler method is used in time, and the potential V has a single jump, Scheme I is positive and l∞ -contracting under the CFL condition: ¯ − ¯ ¯ V 1 −V + 1 ¯ ¯ i+ 2 i− 2 ¯ maxi ¯ ¯ ∆x ¯ ¯ maxj |ξj | < 1. ∆t + (2.6) ∆ξ ∆x In [29] we proved that when the potential is a step function, the same scheme is l1 stable under the CFL condition (2.6) and a suitable condition on the numerical initial data. ¯ − ¯ ¯ V 1 −V + 1 ¯ ¯ i+ 2 i− 2 ¯ Note that the quantity ¯ ¯ represents the gradient of the potential at its ∆x ¯ ¯ smooth point, which has a finite upper bound. Thus the scheme satisfies a hyperbolic CFL condition. 3. The l1 -error estimates for the immersed interface upwind scheme. In this Section we present the l1 -error estimates for the immersed interface upwind scheme established in [30, 28] to the linear advection equation ∂u ∂u + c(x) = 0, ∂t ∂x
t > 0, x ∈ R,
(3.1)
6
X. WEN AND S. JIN
u|t=0 = u0 (x),
(3.2)
with a step function wave speed ( c(x) =
c−
x0
c
.
(3.3)
We assume that c(x) has a definite sign. This avoids the mathematical difficulties of dealing with non-uniqueness or measure-valued solutions, see for examples [3, 23]. Without loss of generality we assume c(x) > 0. We consider the following interface condition for (3.1)-(3.3) requiring u being continuous across the interface u(0+ , t) = u(0− , t).
(3.4)
In [30] we have considered a more general class of interface conditions including (3.4). In this paper we only present the results corresponding to the interface condition (3.4), which will be used to derive the l1 -error estimates for Scheme I. The exact solution of (3.1)-(3.3) with the interface condition (3.4) can be constructed using the method of characteristics: u0 (x − c− t) x c+ t Consider the uniform mesh introduced in Section 2 in x and t-directions. We assume that the interface x = 0 is located at a grid point. We introduce quantities ∆t ∆t λ− = c− ∆x , λ+ = c+ ∆x . The condition 0 < λ− , λ+ < 1 is the CFL condition. The immersed interface upwind scheme which builds the interface condition (3.4) into the upwind difference scheme (with the forward Euler time discretization) for the equation (3.1)-(3.3) reads n vin+1 = (1 − λ− )vin + λ− vi−1 ,
if xi < 0,
(3.6)
vin+1
if xi > 0,
(3.7)
= (1 − λ
+
)vin
where vi0 =
1 ∆x
+
Z
n λ+ vi−1 ,
xi+ 1 2
u0 (x)dx.
(3.8)
xi− 1 2
To compare the numerical solution computed from (3.6)-(3.8) with the exact solution (3.5), we introduce v(x, t) = vin ,
for (x, t) ∈ [xi− 12 , xi+ 12 ) × [tn , tn+1 ).
The following Theorem was proven in [30] Theorem 3.1. The upwind difference scheme (3.6)-(3.8), under the CFL condition 0 < λ− , λ+ < 1, has the following l1 -error to the exact solution (3.5): · µ ½ + ¾¶ ¸ cM √ cM c ∆x , kv(·, tn ) − u(·, tn )kl1 (R) ≤ ku0 kBV (R) γm m ∆x + 1 + m + max − , 1 c c c (3.9)
7
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
where s m
c
−
+
= min{c , c },
c
M
−
+
γm =
= max{c , c },
µ ¶ 2 m ∆t c 1 − cm tn+1 , e ∆x
and the definition of the BV norm on R is given by 1 ku0 (· + h) − u0 (·)kl1 (R) . |h| |h|6=0
ku0 kBV (R) = sup
(3.10)
Remark 3.1. In the case of c− = c+ , namely a constant wave speed, the coefficient before ∆x in (3.9) can be sharpened. In this paper we will also use (3.9) for the constant wave speed case since this will not affect the leading order coefficient in the l1 -error estimates for Scheme I. 4. Setup of the problem and the main result. In this Section we establish the l1 -error estimates for Scheme I (with the first order numerical flux and the forward Euler method in time) under suitable conditions on the initial data. We consider the case when V (x) is a step function, with a jump −D, D > 0 at x = 0. Namely V (0− ) − V (0+ ) = D. Let the computational domain be confined in the rectangular domain (2.2). We ∆t employ the uniform mesh introduced in Section 2. Define mesh ratios λtx = ∆x , λξx = ∆ξ 1 ∆x , assumed to be fixed. Let the potential barrier x = 0 be at a grid point xm+ 2 . Then the point-wise values of V satisfy − + Vm+ = D, 1 − V m+ 1 2
2
± − Vi+ , i < m, 1 = V m+ 1 2
2
± + Vi+ , i > m, 1 = V m+ 1 2
2
where the superscripts −, + represent the left and √ right limits at √ x = 0. We consider the typical situation that ξ1 < − 2D, ξM > 2D, so that all possible particle behaviors, including both transmission and reflection, are included. To √ simplify the discussion, we choose the mesh such that 0 and ± 2D are grid points in the ξ-direction. Define the indices I0 , I+ satisfying √ ξI0 + 21 = 0, ξI+ + 12 = 2D. Define the domain
(
Db =
r³
(x, y)|x < 0, ξ < −
ξ 12
´2
) − 2D
.
Due to the velocity change across the barrier at x = 0, Db represents the area where particles come from outside of the domain [x 12 , xN + 21 ] × [ξ 12 , ξM + 12 ]. In order to implement Scheme I conveniently, we need to exclude this domain from the computational domain. Define the index Ib satisfying r³ ´ ξIb − 32 < − and the domain
ξ 12
2
− 2D ≤ ξIb − 12
o n b b = (x, ξ)|x < 0, ξ < ξI − 1 . D b 2
8
X. WEN AND S. JIN
Then we choose the computational domain as b b. DC = {(x, ξ)|x 21 ≤ x ≤ xN + 12 , ξ 12 ≤ ξ ≤ ξM + 12 } \ D b b. Figure 4.1 depicts DC and D We consider the Dirichlet boundary conditions at the incoming boundaries and assume that the initial data satisfy these boundary conditions: f (x 12 , ξ, t) = f (x 12 , ξ, 0),
0 < ξ < ξM + 12 ,
f (xN + 12 , ξ, t) = f (xN + 12 , ξ, 0),
ξ 12 < ξ < 0.
(4.1) (4.2)
Remark 4.1. When V is a step function, (2.1) becomes ft + ξfx = 0 where ξ only serves as a parameter. However, when implementing the interface conditions (2.4), (2.5) numerically, the mesh in ξ is needed. This is the main new difficulty when comparing with the problem presented in Section 3. 4.1. The initial data assumption. We now impose assumptions on the initial data. We assume the initial data are given on the rectangular domain (2.2). We have the following assumption: Assumption 4.1. The initial data f (x, ξ, 0) have bounded variation in the x-direction and is Lipschitz continuous in the ξ-direction. Namely ∀ξ ∈ [ξ 21 , ξM + 12 ], h i h i |f (x, ξ 0 , 0) − f (x, ξ 00 , 0)| ≤ B|ξ 0 − ξ 00 |, ∀x ∈ x 21 , xN + 12 , ξ 0 , ξ 00 ∈ ξ 21 , ξM + 12 .
kf (., ξ, 0)kBV ([x 1 ,xN + 1 ]) ≤ A, 2
(4.3)
2
(4.4)
Remark 4.2. When arising in the semiclassical limit of the linear Schr¨ odinger equation, the Liouville equation is supplied with the measure-valued initial data [8, 18], which does not satisfy Assumption 4.1. However, in [10], a decomposition of the initial data was introduced, which allows one to solve the semiclassical limit problem with bounded initial data which do satisfy Assumption 4.1. Specifically, if one needs to solve the Liouville equation (2.1) with the measure-valued initial data f (x, ξ, 0) = ρ0 (x)δ(ξ − v0 (x)) ,
(4.5)
the decomposition technique proposed in [10] suggests that one just solves two functions satisfying the same Liouville equation with initial data φ(x, ξ, 0) = ρ0 (x) ,
ψ(x, ξ, 0) = ξ − v0 (x) .
(4.6)
Then the measure-valued solution to the Liouville equation with the initial data (4.5) is simply f (x, ξ, t) = φ(x, ξ, t)δ(ψ(x, ξ, t)). The evaluations of the moments of f then resort to the numerical approximations to the delta function integrals.
9
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
In solving for φ, ψ, the initial data (4.6) satisfy Assumption 4.1 if the initial density and velocity ρ0 (x), v0 (x) have bounded variations. For this situation we have established the l1 -error estimates for Scheme I. Remark 4.3. It can be checked that the initial data satisfying Assumption 4.1 is bounded in both l∞ and l1 -norms on DC . When the numerical initial data are chosen to be cell averages of such initial data, the l1 -stability of Scheme I with homogeneous Dirichlet boundary conditions was established in [29]. 4.2. The exact solution. For our proof, we need to introduce the partition of DC as o n o n ¯ ¯ Dl+ = (x, ξ) ¯ x 12 < x < 0, 0 < ξ < ξM + 12 , Dl− = (x, ξ) ¯ x 12 < x < 0, ξIb − 12 < ξ < 0 , n o n √ √ o ¯ ¯ Dr+ = (x, ξ) ¯ 0 < x < xN + 21 , 2D < ξ < ξM + 12 , Dr− = (x, ξ) ¯ 0 < x < xN + 12 , ξ 21 < ξ < − 2D , n √ o √ ¯ Drr = (x, ξ) ¯ 0 < x < xN + 12 , − 2D < ξ < 2D and the extension of the initial data f (x, ξ, 0) x 12 ≤ x ≤ xN + 12 f (x 12 , ξ, 0) x < x 12 fb0 (x, ξ) = x > xN + 12 f (xN + 21 , ξ, 0)
for x ∈ R, ξ 21 ≤ ξ ≤ ξM + 12 . (4.7)
Figure 4.1 shows a sketch of the partition of DC . ξM + 1 2
Dr+ Dl+
ξI+ + 1 =
√
2
DC
0 = ξI0 + 1 2
Drr √ − 2D
Dl− Dr−
ξI b − 1 2
bb D
ξ1 2
x1
xm+ 1 = 0
2
xN + 1
2
Figure 4.1
2
b b. Sketch of partition of DC and D
The exact solution of (2.1) with the step function potential V (x), when using the interface conditions (2.4) (2.5), can be obtained from the initial data f (x, ξ, 0) and boundary conditions (4.1), (4.2) by the method of characteristics: 1) for (x, ξ) ∈ Dl+ ∪ Dr− , f (x, ξ, tn ) = fb0 (x − ξtn , ξ) ; 2) for (x, ξ) ∈ Drr , f (x, ξ, tn ) =
(
fb0 (x − ξtn , ξ) fb0 (−x + ξtn , −ξ)
(4.8)
x > ξtn ; x < ξtn
(4.9)
2D
10
X. WEN AND S. JIN
3) for (x, ξ) ∈ Dr+ , µ√ ¶ p p ξ 2 −2D 2 2 fb0 x − ξ − 2D tn , ξ − 2D ξ f (x, ξ, tn ) = b f0 (x − ξtn , ξ)
0 < x < ξtn
;
x > ξtn
4) for (x, ξ) ∈ Dl− , ¶ µ √ p p 2 fb0 − ξ +2D x + ξ 2 + 2D tn , − ξ 2 + 2D ξ f (x, ξ, tn ) = b f0 (x − ξtn , ξ)
(4.10)
ξtn < x < 0
.
x < ξtn (4.11)
4.3. The numerical solution. Denote µj = λtx |ξj | ,
1 ≤ j ≤ M.
(4.12)
Under the CFL condition (2.6), µj < 1 for 1 ≤ j ≤ M . Since Vx (x) = 0 except at x = xm+1/2 , with the boundary conditions (4.1), (4.2), Scheme I on DC is given by: 1) if 0 < ξj < ξM + 12 , i 6= m + 1, n+1 n n gij = (1 − µj )gij + µj gi−1,j ;
(4.13)
2) if ξIb − 12 < ξj < 0, i < m or ξ 12 < ξj < 0, i > m, n+1 n n gij = (1 − µj )gij + µj gi+1,j ;
3) if
√
(4.14)
2D < ξj < ξM + 12 ,
³ ´ n+1 n n n gm+1,j = (1 − µj )gm+1,j + µj cj,1 gm,d + c g ; j,2 m,d +1 j j
4) if 0 < ξj
0, x ∈ R, fb0 (x, ξ)dξ
(5.3) x 0, x ∈ R, ∂t ∂x ( ξJj,p x0 Z ξ 1 j+ 1 2 ujr,+,p |t=0 = f j (x, ξ)dξ, ∆ξ ξj− 1 r,+,p 2 ( b0 (x, ξ − ξj + ξJ ) f x 0, x ∈ R, ∂t ( ¯ ∂x¯ ¯ξ J ¯ x0 Z ξ 1 j+ 1 2 j f j (x, ξ)dξ, ul,−,p |t=0 = ∆ξ ξj− 1 l,−,p 2 ( b f0 (−x, ξ − ξj + ξJj,p ) j fl,−,p (x, ξ) = fb0 (−x, ξ)
(5.74) (5.75) (5.76) x ξ j tn
.
x ≤ ξj t n
Then ¯ ¯ 7 ¯ −,j,p ¯ (x) − x¯ ≤ ∆ξtn ¯ Hξ 2
(5.82)
for x ≤ 0, ξj− 12 ≤ ξ ≤ ξj+ 21 , Ib ≤ j ≤ I0 , p = 1, 2. This Lemma can be proved using similar technique for proving Lemma 5.3 and we omit the detailed proof. Lemma 5.7. Define functions Z 0 ¯ ¯ ¯ ej ¯ j Gp (ξ) = (5.83) ¯fl,−,p (x, ξ, tn ) − fbl,− (x, ξ, tn )¯ dx, x1
2
for ξj− 12 ≤ ξ ≤ ξj+ 12 , Ib ≤ j ≤ I0 , p = 1, 2, where j fel,−,p (x, ξ, t) =
fbl,− (x, ξ, t) =
³ ´ p fb0 ξJj,p x − ξJj,p t, − ξ 2 + 2D ξj
x > ξj t
b f0 (x − ξj t, ξ) x ≤ ξj t µ √ ¶ p p 2 fb0 − ξ +2D x + ξ 2 + 2D t, − ξ 2 + 2D ξ
,
b f0 (x − ξt, ξ)
(5.84) x > ξt
.
(5.85)
x ≤ ξt
Then ³ ´ √ Gjp (ξ) ≤ 2A + 2DB 7∆ξtn . Proof. From definition (5.85) and conditions (4.3), (4.4) one obtains √ kfbl,− (·, ξ, t)kBV (R) ≤ 2A + 2DB.
(5.86)
(5.87)
We use the fact that ³ ´ j fel,−,p (x, ξ, tn ) = fbl,− Hξ−,j,p (x), ξ, tn ,
x ∈ R, p = 1, 2.
(5.88)
From (5.83), (5.88) and (5.87), applying Lemma 5.6 and Lemma A.1 in the Appendix gives (5.86).
22
X. WEN AND S. JIN
We now give estimates for El− . In comparison with schemes (4.14), (4.17), (4.20), (4.21), (4.23) and (5.79), (5.80), (5.76), (5.77), (4.7), one can check that for Ib ≤ j ≤ I0 g(x, ξ, tn ) =
2 X
j cj,p vl,−,p (−x, tn ),
p=1
a.e. for ξj− 12 < ξ < ξj+ 21 , x 12 < x < 0.
(5.89)
Together with (5.73), (5.89) and the expression (4.11) of f (x, ξ, tn ) in Dl− , one obtains El− ≤ El−,1 + El−,2 + El−,3 ,
(5.90)
where El−,1
=
I0 Z X
=
ξj+ 1
Z
2
ξj− 1 2
I0 Z X j=Ib
by
2
I0 Z X j=Ib
El−,3 =
Z
2
ξj− 1
j=Ib
El−,2
ξj+ 1
ξj+ 1
2
ξj− 1 2
Z
2 X
0
x 1 p=1 2 2 X
0
x 1 p=1 2 2 X
0
x 1 p=1 2
¯ ¯ ¯ j ¯ cj,p ¯vl,−,p (−x, tn ) − ujl,−,p (−x, tn )¯ dxdξ,
(5.91)
¯ ¯ ¯ ¯ j cj,p ¯ujl,−,p (−x, tn ) − fel,−,p (x, ξ, tn )¯ dxdξ,
(5.92)
¯ ¯ ¯ j ¯ cj,p ¯fel,−,p (x, ξ, tn ) − fbl,− (x, ξ, tn )¯ dxdξ.
(5.93)
j Applying Theorem 3.1, the l1 -error between vl,−,p and ujl,−,p , p = 1, 2 are given
"r j kvl,−,p (·, tn )−ujl,−,p (·, tn )kl1 (R)
≤
kujl,−,p (·, 0)kBV (R)
# ¯ ¯ ¯ ¯¶ µ ¯ ξJj,p ¯ 2tn+1 ¯ξJj,p ¯ √ ¯ ∆x . p ∆x + 2 + ¯¯ e ξj ¯ |ξj | (5.94)
Similar to the derivation of (5.64), one has kujl,−,p (·, 0)kBV (R) ≤ 2A +
√
2DB.
(5.95)
Applying (5.94) and (5.95) to (5.91), by summation of j and using the inequalities (5.65) and (5.70) one obtains El−,1
´ r 2t
¸ 14 ¯ ¯ · ³ ´2 √ ¯ ¯ 1 1 ξ ξ − 2D ∆x ¯ 2¯ 2 e ³ ´¯ ¯ µ 1 ¶ √ ¯ ¯ + 2A + 2DB ¯ξ 12 ¯ ln ∆x + O(∆x). ∆x
³
√
≤ 4A + 2 2DB
n+1
(5.96)
Using Lemma 5.5, (5.67) and condition (4.4), from the definitions (5.76), (5.77), (5.78) and (5.84), one has ¯ ¯ ¯ j ¯ j ¯ul,−,p (−x, tn ) − fel,−,p (x, ξ, tn )¯ ≤ 2B∆ξ, for x ≤ 0, ξj− 21 ≤ ξ ≤ ξj+ 12 . Thus one has ¯ ¯¯ ¯ ¯ ¯¯ ¯ El−,2 ≤ ¯ξIb − 12 ¯ ¯x 12 ¯ 2B∆ξ.
(5.97)
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
23
Applying Lemma 5.7 to (5.93) yields ¯ ³ ´¯ √ ¯ ¯ El−,3 ≤ 7 2A + 2DB ¯ξIb − 12 ¯ tn ∆ξ.
(5.98)
Combining (5.90), (5.96), (5.97), (5.98) leads to El−
´ r 2t
¸ 14 ¯ ¯ ·³ ´2 √ ¯ ¯ ≤ 4A + 2 2DB ∆x ¯ξ 12 ¯ ξ 12 − 2D e µ ¶ ³ ´¯ ¯ √ 1 ¯ ¯ + 2A + 2DB ¯ξ 12 ¯ ln ∆x + O(∆x). ∆x ³
√
n+1
(5.99)
Finally combining (5.1), (5.20), (5.21), (5.33), (5.72), (5.99), (5.12) completes the proof for Theorem 4.1. Remark The error terms derived in (4.25) include the halfth order terms ¡ 1 5.1. ¢ and the ln ∆x ∆x terms. Thus the leading error terms are of halfth order and Scheme I has a halfth order convergence rate. Note that this is sharp, since even for the discontinuous solution to linear hyperbolic equation with a smooth coefficient, one cannot expect a better convergence order in the l1 -norm [24]. 6. Conclusion. In this paper we derived the l1 -error estimates for a Hamiltonianpreserving scheme, developed in [11], for the Liouville equation with a piecewise constant potential in one space dimension. The Hamiltonian-preserving scheme is designed by incorporating into the numerical fluxes the particle behavior–transmission and reflection– at the potential barrier. We proved that, with the Dirichlet incoming boundary conditions and for a class of bounded initial data, the numerical solution by the Hamiltonian-preserving scheme converges in l1 -norm to the solution of the Liouville equation–defined by using the interface condition that accounts for particle transmission and reflection. The initial data conditions can be satisfied by applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data arising in the semiclassical limit of the linear Schr¨odinger equation. The strategy for the error analysis in this paper is to apply the l1 -error estimates established in [30, 28] for the immersed interface upwind scheme to the linear advection equations with piecewise constant coefficients. This is a natural approach since the solution of the Liouville equation with a step function potential satisfies linear advection equations with piecewise constant coefficients on the bicharacteristics. To apply this strategy, we split the computational domain into several parts. In each subdomain, we introduced a number of linear advection equations with step function coefficients to which the immersed interface upwind schemes yield the same numerical solution as the Hamiltonian-preserving scheme for the Liouville equation with the step function potential. Then the l1 -error estimates for the Hamiltonian-preserving scheme were derived by applying the l1 -error estimates for the immersed interface upwind schemes and comparing the exact solutions of the linear advection equations and that of the Liouville equation in each of these subdomains. The condition on the initial data were used in deriving these estimates. As a result, we obtained the half order l1 -error bound with explicit coefficients for the Hamiltonian-preserving scheme. The problem under study has important applications in the computation of the semiclassical limit of the linear Schr¨odinger equation through barriers, and more generally, in the computation of high frequency waves through interfaces.
24
X. WEN AND S. JIN
Appendix In this Appendix we prove a property of BV functions on R. Lemma A.1. Let f (x) be a BV function on R, H(x) be a function on [a, b] satisfying |H(x) − x| ≤ HC ,
for x ∈ [a, b],
where HC is a positive constant. Then kf (·) − f (H(·))kL1 ([a,b]) ≤ 2HC kf kBV (R) .
(A.1)
Proof. ∀x ∈ [a, b], one has |f (x) − f (H(x))| ≤ kf kBV ([x−HC ,x+HC ]) . Since f (x) has a bounded variation on R, f 0 (x) exists a.e., and Z kf kBV ([x−HC ,x+HC ]) =
x+HC
|f 0 (y)| dy.
x−HC
Thus Z kf (·) − f (H(·))kL1 ([a,b]) ≤
b
Z
a
x+HC
|f 0 (y)| dydx.
(A.2)
x−HC
Introduce the function ( g(x, y) =
|f 0 (y)|
if |x − y| ≤ HC
0
else
.
Then g(x, y) ∈ L1 ([a, b] × [a − HC , b + HC ]). Applying Fubini Theorem gives Z Z
b a
Z
x+HC
Z |f 0 (y)| dydx =
x−HC b+HC Z b
=
b
Z
b+HC
g(x, y)dydx a a−HC Z b+HC Z y+HC
g(x, y)dxdy ≤ a−HC
Z
a b+HC
= 2HC a−HC
a−HC
|f 0 (y)| dxdy
y−HC
|f 0 (y)| dy ≤ 2HC kf kBV (R) .
(A.3)
Combining (A.2) and (A.3) gives (A.1). REFERENCES e, and G.D.V. Gowda, Godunov-type methods for conservation laws with [1] Adimurthi, J. Jaffr´ a flux function discontinuous in space, SIAM J. Numer. Anal. 42(1), 179-208, 2004. [2] L. Ambrosio, Transport equation and Cauchy problem for BV vector fields, Invent. Math., 158, 227-260, 2004. [3] F. Bouchut and F. James, One-dimensional transport equations with discontinuous coefficients. Nonlinear Anal. 32 (1998), no. 7, 891–933. [4] G.M. Coclite and N.H. Risebro, Conservation laws with time dependent discontinuous coefficients, SIAM J. Math. Anal. 36(4), 1293-1309, 2005.
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
25
[5] R.J. DiPerna and P.-L. Lions, Ordinary differential equations, transport theory, and Sobolev spaces, Invent. Math. 98, 511-547, 1989. [6] T. Gimse, Conservation laws with discontinuous flux functions, SIAM J. Math. Anal. 24(2), 279-289, 1993. [7] L. Gosse and F. James, Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients, Math. Comp. 69(231), 987-1015, 2000. [8] P. G´ erard, P.A. Markowich, N.J. Mauser and F. Poupaud, Homogenization limits and Wigner transforms, Comm. Pure Appl. Math., 50(1997), pp. 321-377. [9] T. Gimse and N.H. Risebro, Solution of the Cauchy problem for a conservation law with a discontinuous flux function, SIAM J. Math. Anal. 23, 635-648, 1992. [10] S. Jin, H.L. Liu, S. Osher and R. Tsai, Computing multivalued physical observables for the semiclassical limit of the Schrodinger equation, J. Comp. Phys. 205, 222-241, 2005. [11] S. Jin and X. Wen, The Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials, Comm. Math. Sci. 3, 285-315, 2005. [12] S. Jin and X. Wen, A Hamiltonian-preserving scheme for the Liouville equation of geometrical optics with partial transmissions and reflections, SIAM J. Num. Anal. 44, 1801-1828, 2006. [13] K.H. Karlsen, C. Klingenberg and N.H. Risebro, A relaxation scheme for conservation laws with a discontinuous coefficient, Math. Comp. 73(247), 1235-1259,2004. [14] C. Klingenberg and N.H. Risebro, Convex conservation laws with discontinuous coefficients. Existence, uniqueness and asymptotic behavior, Comm. Partial Differential Equations 20(11-12), 1959-1990,1995. [15] C. Klingenberg and N.H. Risebro, Stability of a resonant system of conservation laws modeling polymer flow with gravitation, J. Differential Equations 170(2), 344-380,2001. [16] K.H. Karlsen, N.H. Risebro and J.D. Towers, Upwind difference approximations for degenerate parabolic convection-diffusion equations with a discontinuous coefficient, IMA J. Numer. Anal. 22(4), 623-664,2002. [17] K.H. Karlsen and J.D. Towers, Convergence of the Lax-Friedrichs scheme and stability for conservation laws with a discontinuous space-time dependent flux, Chinese Ann. Math. Ser. B. 25(3), 287-318, 2004. [18] P.L. Lions and T. Paul, Sur les measures de Wigner, Revista. Mat. Iberoamericana 9, 1993, 553-618. [19] L.W. Lin, B.J. Temple, and J.H. Wang, A comparison of convergence rates for Godunov’s method and Glimm’s method in resonant nonlinear systems of conservation laws, SIAM J. Numer. Anal. 32(3), 824-840,1995. [20] S. Mishra, Convergence of upwind finite difference schemes for a scalar conservation law with indefinite discontinuities in the flux function, SIAM J. Numer. Anal. 43(2), 559-577, 2005. [21] D.N. Ostrov, Viscosity solutions and convergence of monotone schemes for synthetic aperture radar shape-from-shading equations with discontinuous intensities, SIAM J. Appl. Math. 59(6), 2060-2085, 1999. [22] B. Perthame and C. Simeoni, A kinetic scheme for the Saint-Venant system with a source term, CALCOLO 38(4), 201-231 (2001). [23] F. Poupaud and M. Rascle, Measure solutions to the linear multi-dimensional transport equation with non-smooth coefficients, Comm. Partial Differential Equations 22 (1997), no. 1-2, 337–358. √ [24] T. Tang and Z.H. Teng, The sharpness of Kuznetsov’s O( ∆x) L1-error estimate for monotone difference schemes, Math. Comp. 64, 581-589, (1995). [25] B. Temple, Global solution of the Cauchy problem for a class of 2 × 2 nonstrictly hyperbolic conservation laws, Adv. in Appl. Math. 3(3), 335-375, 1982. [26] J.D. Towers, Convergence of a difference scheme for conservation laws with a discontinuous flux, SIAM J. Numer. Anal. 38(2), 681-698, 2000. [27] J.D. Towers, A difference scheme for conservation laws with a discontinuous flux - the nonconvex case, SIAM J. Numer. Anal. 39(4), 1197-1218, 2001. [28] X. Wen, Convergence of an immersed interface upwind scheme for linear advection equations with piecewise constant coefficients II: Some related binomial coefficients inequalities, preprint. [29] X. Wen and S. Jin, The L1 -stability of a Hamiltonian-preserving scheme for the Liouville equation with discontinuous potentials, preprint. [30] X. Wen and S. Jin, Convergence of an immersed interface upwind scheme for linear advection equations with piecewise constant coefficients I: L1 -error estimates, J. Comput. Math. 26(1), 1-22, 2008.
† AND
SHI JIN
‡
Abstract. We study the l1 -error of a Hamiltonian-preserving scheme, developed in [11], for the Liouville equation with a piecewise constant potential in one space dimension. This problem has important applications in computations of the semiclassical limit of the linear Schr¨ odinger equation through barriers, and of the high frequency waves through interfaces. We use the l1 -error estimates established in [30, 28] for the immersed interface upwind scheme to the linear advection equations with piecewise constant coefficients. We prove that the scheme with the Dirichlet incoming boundary conditions is l1 -convergent for a class of bounded initial data, and derive the one-halfth order l1 -error bounds with explicit coefficients. The initial conditions can be satisfied by applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data, which arises in the semiclassical limit of the linear Schr¨ odinger equation. Key words. Liouville equations, Hamiltonian preserving schemes, piecewise constant potentials, error estimate, half order error bound, semiclassical limit AMS subject classifications. 65M06, 65M12, 65M25, 35L45, 70H99
1. Introduction. In [11], we constructed a class of numerical schemes for the d-dimensional Liouville equation in classical mechanics: ft + v · ∇x f − ∇x V · ∇v f = 0 ,
t > 0,
x, v ∈ Rd ,
(1.1)
where f (t, x, v) is the density distribution of a classical particle at position x, time t and traveling with velocity v. V (x) is the potential. The main interest is in the case of a discontinuous potential V (x), corresponding to a potential barrier. When V is discontinuous, the Liouville equation (1.1) is a linear hyperbolic equation with a measure-valued coefficient. Such a problem cannot be understood mathematically using the renormalized solution by DiPerna and Lions for linear advection equations with discontinuous coefficients [5] (see also [2]). Our approach in [11, 12] to such problems was to provide an interface condition to couple the Liouville equation (1.1) on both sides of the barrier or interface. The interface condition accounts for particle or wave transmission and reflection. Based on this notion of the solution, the socalled Hamiltonian-preserving schemes were constructed in [11, 12], which build this interface condition into the numerical flux. Schemes so constructed provide solutions that are physically relevant for particle or wave reflection and transmission through the barriers or interfaces. The Liouville equation is the phase space representation of Newton’s second law: dx = v, dt
dv = −∇x V , dt
∗ Research supported in part by the Knowledge Innovation Project of the Chinese Academy of Sciences grants K5501312S1, K5502212F1, K7290312G7 and K7502712F7, NSFC grant 10601062, NSF grant DMS-0608720 and NSAF grant 10676017. † Institute of Computational Mathematics, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, P.O. Box 2719, Beijing 100190, China. Email address: [email protected]. ‡ Department of Mathematics, University of Wisconsin, Madison, WI 53706, USA. Email address: [email protected].
1
2
X. WEN AND S. JIN
which is a Hamiltonian system with the Hamiltonian H=
1 2 |v| + V (x). 2
It is known from classical mechanics that the Hamiltonian remains constant across a potential barrier. This is one of the main ingredients in the Hamiltonian-preserving schemes developed in [11, 12]. The two schemes developed in [11]–one based on a finite difference formulation (called Scheme I) and the other on a finite volume formulation (called Scheme II) were proved, in one space dimension with a piecewise constant potential, to be positive, and l1 and l∞ -stable for suitable initial value problems and under a hyperbolic CFL condition (see also [29]). The more sophisticated issue of the l1 -stability of Scheme I was studied in [29]. We proved in [29] that, in the case of a step function potential, Scheme I with the homogeneous Dirichlet incoming boundary conditions is l1 -stable under a certain condition on the initial data. In particular, we showed that such an initial condition is satisfied when choosing the numerical initial data as the cell averages of a bounded analytical initial data. Since the exact solution, in the aforementioned notion using the interface condition, to the same problem is l1 -contracting, an l1 -convergent scheme should be l1 -stable. Thus the results established in [29] implies that Scheme I meets the necessary condition for the l1 -convergence in the case of the homogeneous Dirichlet boundary conditions and bounded initial data. In this paper we will study the l1 -convergence of Scheme I in the case of a step function potential and the Dirichlet incoming boundary conditions. We will prove that in this case Scheme I is l1 -convergent for a class of bounded initial data. The initial conditions can be satisfied when applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data arising in the semiclassical limit of the linear Schr¨odinger equation. The Liouville equation with a step function (or more generally piecewise constant) potential belongs to hyperbolic equations with singular (discontinuous or measurevalued) coefficients. For the discontinuous coefficient case, the convergence of numerical schemes were widely studied [6, 9, 14, 15, 4, 17, 25, 19, 21, 7, 26, 27, 16, 13, 1, 20]. However, in most cases the convergence rate estimates for numerical schemes were not studied. For the measure-valued coefficient case, recently we have obtained the halfth order l1 -error estimates for the immersed interface upwind scheme which builds the interface condition into the numerical flux, see [30, 28]. In this paper, we extend this result to the more interesting case of the Liouville equation with a step function potential. Since the Liouville equation is defined in the phase space, extra efforts are needed to deal with the discretization in the velocity space. To derive the desired error estimates, we split the computational domain into several parts. In each subdomain, we will introduce a number of linear advection equations with step function coefficients. The immersed interface upwind schemes for these linear advection equations yield the same numerical solution as Scheme I for the Liouville equation with the step function potential. Then the l1 -error estimates for Scheme I can be achieved by applying the l1 -error estimates for the immersed interface upwind schemes and estimating the l1 -errors between the exact solutions of the linear advection equations with step function coefficients and that of the Liouville equation with the step function potential. The first part of the error estimates are obtained by applying the results established in [30, 28] with the aid of the condition on the initial data. The second part of the errors are analyzed also utilizing the condition on the
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
3
initial data. With this approach, we can derive the halfth order l1 -error bound with explicit coefficients for Scheme I. Our proof is based on the step function potential. This result clearly generates to piecewise constant potentials in a straightforward but tedious way. The same approach should also be applicable to analyze the l1 -error estimates for the finite volume Scheme II developed in [11]. But we will not pursue this in this paper. This paper is organized as follows. In Section 2 we review the Hamiltonianpreserving scheme proposed in [11] for the Liouville equation with a discontinuous potential in one space dimension. In Section 3 we recall the l1 -error estimates for the immersed interface upwind scheme established in [30, 28] to the linear advection equation with a step function wave speed. In Section 4 we setup the problem and state the main Theorem of this paper which is proved in Section 5 using the results presented in Section 3. We conclude the paper in Section 6. 2. A Hamiltonian-preserving scheme. In this Section we review the Hamiltonianpreserving scheme proposed in [11] to the Liouville equation in one space dimension ft + ξfx − Vx fξ = 0
(2.1)
with a discontinuous potential V (x). Consider a uniform mesh with grid points at xi+ 12 , i ∈ Z in the x-direction and ξj+ 12 , j ∈ Z in the ξ-direction. The cells are centered at (xi , ξj ), (i, j) ∈ Z2 with xi = 21 (xi+ 12 + xi− 12 ) and ξj = 12 (ξj+ 12 + ξj− 12 ). The mesh size are denoted by ∆x = xi+ 12 − xi− 12 , ∆ξ = ξj+ 21 − ξj− 12 . We also assume a uniform time step ∆t and the discrete times are given by tn = n∆t, n ∈ N ∪ {0}. We assume that the computation is performed in a bounded rectangular domain {(x, y)|x 21 ≤ x ≤ xN + 12 , ξ 21 ≤ ξ ≤ ξM + 12 }.
(2.2)
Let the cell averages of f be 1 fij = ∆x∆ξ
Z
xi+ 1
Z
2
ξj+ 1 2
xi− 1
f (x, ξ, t)dξdx.
ξj− 1
2
2
The 1-d average quantity fi+1/2,j is defined as fi+1/2,j =
1 ∆ξ
Z
ξj+ 1 2
ξj− 1
f (xi+1/2 , ξ, t)dξ .
2
fi,j+1/2 is defined similarly. In classical mechanics, a particle will either cross a potential barrier with a changing momentum, or be reflected, depending on its momentum and the strength of the potential barrier. Figure 2.1 shows the typical situations when a particle moves from left to right at a potential barrier. If the potential is reduced, the particle will cross it with an accelerated speed. If the potential is increased, the particle is either reflected if its momentum is not big enough to overcome the barrier or transmitted otherwise with a reduced velocity. The Hamiltonian H = 12 ξ 2 + V is preserved across the potential barrier: 1 + 2 1 (ξ ) + V + = (ξ − )2 + V − , 2 2
(2.3)
4
X. WEN AND S. JIN
2
+
−
[ξ −2(V
ξ
1/2
− V )]
1)
2)
ξ
[ξ2+2(V− − V+)]1/2
−ξ
3)
ξ
V−
V−
V+
V+
Figure 2.1 Transmission and reflection of a particle at a potential barrier. where the superscripts ± indicate the right and left limits of the quantity at the potential barrier. This property was used in [11] to provide the interface condition for (2.1) at the barrier: f (t, x+ , ξ + ) = f (t, x− , ξ − ) ±
±
±
for transmission ±
f (t, x , ξ ) = f (t, x , −ξ )
(2.4)
for reflection
(2.5)
where ξ + and ξ − are related by the constant Hamiltonian condition (2.3) in the case of transmission. With such an interface condition, we established the well-posedness of the initial value problem to (1.1) with a piecewise constant potential in [12]. The main ingredient in the Hamiltonian-preserving schemes developed in [11], like the earlier work for shallow-water equations [22], was to build into the numerical flux the interface conditions (2.4), (2.5) at the barrier. We now present the first Hamiltonian-preserving scheme, called Scheme I in [11]. Assume that the discontinuous points of the potential V are located at the grid − + points. Let the left and right limits of V at point xi+1/2 be Vi+ respectively. 1 and V i+ 1 2
2
− + Note that if V is continuous at xj+1/2 , then Vi+ . We approximate V by a 1 = V i+ 12 2 piecewise linear function
V (x) ≈
+ Vi−1/2
+
− + Vi+1/2 − Vi−1/2
∆x
(x − xi−1/2 ) .
The flux-splitting, semidiscrete scheme (with time continuous) reads ∂t fij + ξj
− + fi+ − fi− 1 1 ,j ,j 2
2
∆x
−
− + Vi+ 1 − V i− 12 fi,j+ 12 − fi,j− 21 2
∆x
∆ξ
= 0,
where the numerical fluxes fi,j+ 12 are defined using the upwind discretization. Since the characteristics of the Liouville equation may be different on the two sides of a barrier, the corresponding numerical fluxes should also be different. The essential + − , fi− at each cell part of the algorithm is to define the split numerical fluxes fi+ 1 1 2 ,j 2 ,j interface. (2.4) will be used to define these fluxes. Assume V is discontinuous at xi+1/2 . Consider the case ξj > 0. Using upwind − scheme, fi+ = fij . However, 1 ,j 2
+ − + − fi+1/2,j ≡ f (x+ i+1/2 , ξ ) = f (xi+1/2 , ξ )
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
5
where ξ − is obtained from ξ + = ξj from (2.3). Since ξ − may not be a grid point, we have to define it approximately. The first approach is to locate the two cell centers that bound this velocity, then use a linear interpolation to evaluate the needed numerical flux at ξ − . The case of ξj < 0 is treated similarly. The algorithm to generate the numerical flux is given in [11]. Here we present the simplified algorithm for the case + − being discussed in this paper. Vi+ 1 > V i+ 12 2 Algorithm I • ξj > 0 − = fij , fi+ 1 2 ,j r ³ ´ − + o if ξj > 2 Vi+ , 1 − V 1 i+ 2 2 r ³ ´ + − ξ − = ξj2 + 2 Vi+ 1 − V 1 i+ 2
2
if ξk ≤ ξ − < ξk+1 for some k − − −ξk + then fi+ = ξk+1∆ξ−ξ fik + ξ ∆ξ fi,k+1 1 ,j o else
2
+ fi+ = fi+1,k where ξk = −ξj 1 ,j 2
o end • ξj < 0 + fi+ = fi+1,j , 1 2 ,j r ³ ´ − + ξ + = − ξj2 + 2 Vi+ 1 − V 1 i+ 2
2
if ξk ≤ ξ + < ξk+1 for some k + −ξ + −ξk − fi+1,k + ξ ∆ξ fi+1,k+1 then fi+ = ξk+1 1 ∆ξ 2 ,j After the spatial discretization is specified, one can use any time discretization for the time derivative. In [11] we proved that, when the first order upwind scheme is used spatially, and the forward Euler method is used in time, and the potential V has a single jump, Scheme I is positive and l∞ -contracting under the CFL condition: ¯ − ¯ ¯ V 1 −V + 1 ¯ ¯ i+ 2 i− 2 ¯ maxi ¯ ¯ ∆x ¯ ¯ maxj |ξj | < 1. ∆t + (2.6) ∆ξ ∆x In [29] we proved that when the potential is a step function, the same scheme is l1 stable under the CFL condition (2.6) and a suitable condition on the numerical initial data. ¯ − ¯ ¯ V 1 −V + 1 ¯ ¯ i+ 2 i− 2 ¯ Note that the quantity ¯ ¯ represents the gradient of the potential at its ∆x ¯ ¯ smooth point, which has a finite upper bound. Thus the scheme satisfies a hyperbolic CFL condition. 3. The l1 -error estimates for the immersed interface upwind scheme. In this Section we present the l1 -error estimates for the immersed interface upwind scheme established in [30, 28] to the linear advection equation ∂u ∂u + c(x) = 0, ∂t ∂x
t > 0, x ∈ R,
(3.1)
6
X. WEN AND S. JIN
u|t=0 = u0 (x),
(3.2)
with a step function wave speed ( c(x) =
c−
x0
c
.
(3.3)
We assume that c(x) has a definite sign. This avoids the mathematical difficulties of dealing with non-uniqueness or measure-valued solutions, see for examples [3, 23]. Without loss of generality we assume c(x) > 0. We consider the following interface condition for (3.1)-(3.3) requiring u being continuous across the interface u(0+ , t) = u(0− , t).
(3.4)
In [30] we have considered a more general class of interface conditions including (3.4). In this paper we only present the results corresponding to the interface condition (3.4), which will be used to derive the l1 -error estimates for Scheme I. The exact solution of (3.1)-(3.3) with the interface condition (3.4) can be constructed using the method of characteristics: u0 (x − c− t) x c+ t Consider the uniform mesh introduced in Section 2 in x and t-directions. We assume that the interface x = 0 is located at a grid point. We introduce quantities ∆t ∆t λ− = c− ∆x , λ+ = c+ ∆x . The condition 0 < λ− , λ+ < 1 is the CFL condition. The immersed interface upwind scheme which builds the interface condition (3.4) into the upwind difference scheme (with the forward Euler time discretization) for the equation (3.1)-(3.3) reads n vin+1 = (1 − λ− )vin + λ− vi−1 ,
if xi < 0,
(3.6)
vin+1
if xi > 0,
(3.7)
= (1 − λ
+
)vin
where vi0 =
1 ∆x
+
Z
n λ+ vi−1 ,
xi+ 1 2
u0 (x)dx.
(3.8)
xi− 1 2
To compare the numerical solution computed from (3.6)-(3.8) with the exact solution (3.5), we introduce v(x, t) = vin ,
for (x, t) ∈ [xi− 12 , xi+ 12 ) × [tn , tn+1 ).
The following Theorem was proven in [30] Theorem 3.1. The upwind difference scheme (3.6)-(3.8), under the CFL condition 0 < λ− , λ+ < 1, has the following l1 -error to the exact solution (3.5): · µ ½ + ¾¶ ¸ cM √ cM c ∆x , kv(·, tn ) − u(·, tn )kl1 (R) ≤ ku0 kBV (R) γm m ∆x + 1 + m + max − , 1 c c c (3.9)
7
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
where s m
c
−
+
= min{c , c },
c
M
−
+
γm =
= max{c , c },
µ ¶ 2 m ∆t c 1 − cm tn+1 , e ∆x
and the definition of the BV norm on R is given by 1 ku0 (· + h) − u0 (·)kl1 (R) . |h| |h|6=0
ku0 kBV (R) = sup
(3.10)
Remark 3.1. In the case of c− = c+ , namely a constant wave speed, the coefficient before ∆x in (3.9) can be sharpened. In this paper we will also use (3.9) for the constant wave speed case since this will not affect the leading order coefficient in the l1 -error estimates for Scheme I. 4. Setup of the problem and the main result. In this Section we establish the l1 -error estimates for Scheme I (with the first order numerical flux and the forward Euler method in time) under suitable conditions on the initial data. We consider the case when V (x) is a step function, with a jump −D, D > 0 at x = 0. Namely V (0− ) − V (0+ ) = D. Let the computational domain be confined in the rectangular domain (2.2). We ∆t employ the uniform mesh introduced in Section 2. Define mesh ratios λtx = ∆x , λξx = ∆ξ 1 ∆x , assumed to be fixed. Let the potential barrier x = 0 be at a grid point xm+ 2 . Then the point-wise values of V satisfy − + Vm+ = D, 1 − V m+ 1 2
2
± − Vi+ , i < m, 1 = V m+ 1 2
2
± + Vi+ , i > m, 1 = V m+ 1 2
2
where the superscripts −, + represent the left and √ right limits at √ x = 0. We consider the typical situation that ξ1 < − 2D, ξM > 2D, so that all possible particle behaviors, including both transmission and reflection, are included. To √ simplify the discussion, we choose the mesh such that 0 and ± 2D are grid points in the ξ-direction. Define the indices I0 , I+ satisfying √ ξI0 + 21 = 0, ξI+ + 12 = 2D. Define the domain
(
Db =
r³
(x, y)|x < 0, ξ < −
ξ 12
´2
) − 2D
.
Due to the velocity change across the barrier at x = 0, Db represents the area where particles come from outside of the domain [x 12 , xN + 21 ] × [ξ 12 , ξM + 12 ]. In order to implement Scheme I conveniently, we need to exclude this domain from the computational domain. Define the index Ib satisfying r³ ´ ξIb − 32 < − and the domain
ξ 12
2
− 2D ≤ ξIb − 12
o n b b = (x, ξ)|x < 0, ξ < ξI − 1 . D b 2
8
X. WEN AND S. JIN
Then we choose the computational domain as b b. DC = {(x, ξ)|x 21 ≤ x ≤ xN + 12 , ξ 12 ≤ ξ ≤ ξM + 12 } \ D b b. Figure 4.1 depicts DC and D We consider the Dirichlet boundary conditions at the incoming boundaries and assume that the initial data satisfy these boundary conditions: f (x 12 , ξ, t) = f (x 12 , ξ, 0),
0 < ξ < ξM + 12 ,
f (xN + 12 , ξ, t) = f (xN + 12 , ξ, 0),
ξ 12 < ξ < 0.
(4.1) (4.2)
Remark 4.1. When V is a step function, (2.1) becomes ft + ξfx = 0 where ξ only serves as a parameter. However, when implementing the interface conditions (2.4), (2.5) numerically, the mesh in ξ is needed. This is the main new difficulty when comparing with the problem presented in Section 3. 4.1. The initial data assumption. We now impose assumptions on the initial data. We assume the initial data are given on the rectangular domain (2.2). We have the following assumption: Assumption 4.1. The initial data f (x, ξ, 0) have bounded variation in the x-direction and is Lipschitz continuous in the ξ-direction. Namely ∀ξ ∈ [ξ 21 , ξM + 12 ], h i h i |f (x, ξ 0 , 0) − f (x, ξ 00 , 0)| ≤ B|ξ 0 − ξ 00 |, ∀x ∈ x 21 , xN + 12 , ξ 0 , ξ 00 ∈ ξ 21 , ξM + 12 .
kf (., ξ, 0)kBV ([x 1 ,xN + 1 ]) ≤ A, 2
(4.3)
2
(4.4)
Remark 4.2. When arising in the semiclassical limit of the linear Schr¨ odinger equation, the Liouville equation is supplied with the measure-valued initial data [8, 18], which does not satisfy Assumption 4.1. However, in [10], a decomposition of the initial data was introduced, which allows one to solve the semiclassical limit problem with bounded initial data which do satisfy Assumption 4.1. Specifically, if one needs to solve the Liouville equation (2.1) with the measure-valued initial data f (x, ξ, 0) = ρ0 (x)δ(ξ − v0 (x)) ,
(4.5)
the decomposition technique proposed in [10] suggests that one just solves two functions satisfying the same Liouville equation with initial data φ(x, ξ, 0) = ρ0 (x) ,
ψ(x, ξ, 0) = ξ − v0 (x) .
(4.6)
Then the measure-valued solution to the Liouville equation with the initial data (4.5) is simply f (x, ξ, t) = φ(x, ξ, t)δ(ψ(x, ξ, t)). The evaluations of the moments of f then resort to the numerical approximations to the delta function integrals.
9
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
In solving for φ, ψ, the initial data (4.6) satisfy Assumption 4.1 if the initial density and velocity ρ0 (x), v0 (x) have bounded variations. For this situation we have established the l1 -error estimates for Scheme I. Remark 4.3. It can be checked that the initial data satisfying Assumption 4.1 is bounded in both l∞ and l1 -norms on DC . When the numerical initial data are chosen to be cell averages of such initial data, the l1 -stability of Scheme I with homogeneous Dirichlet boundary conditions was established in [29]. 4.2. The exact solution. For our proof, we need to introduce the partition of DC as o n o n ¯ ¯ Dl+ = (x, ξ) ¯ x 12 < x < 0, 0 < ξ < ξM + 12 , Dl− = (x, ξ) ¯ x 12 < x < 0, ξIb − 12 < ξ < 0 , n o n √ √ o ¯ ¯ Dr+ = (x, ξ) ¯ 0 < x < xN + 21 , 2D < ξ < ξM + 12 , Dr− = (x, ξ) ¯ 0 < x < xN + 12 , ξ 21 < ξ < − 2D , n √ o √ ¯ Drr = (x, ξ) ¯ 0 < x < xN + 12 , − 2D < ξ < 2D and the extension of the initial data f (x, ξ, 0) x 12 ≤ x ≤ xN + 12 f (x 12 , ξ, 0) x < x 12 fb0 (x, ξ) = x > xN + 12 f (xN + 21 , ξ, 0)
for x ∈ R, ξ 21 ≤ ξ ≤ ξM + 12 . (4.7)
Figure 4.1 shows a sketch of the partition of DC . ξM + 1 2
Dr+ Dl+
ξI+ + 1 =
√
2
DC
0 = ξI0 + 1 2
Drr √ − 2D
Dl− Dr−
ξI b − 1 2
bb D
ξ1 2
x1
xm+ 1 = 0
2
xN + 1
2
Figure 4.1
2
b b. Sketch of partition of DC and D
The exact solution of (2.1) with the step function potential V (x), when using the interface conditions (2.4) (2.5), can be obtained from the initial data f (x, ξ, 0) and boundary conditions (4.1), (4.2) by the method of characteristics: 1) for (x, ξ) ∈ Dl+ ∪ Dr− , f (x, ξ, tn ) = fb0 (x − ξtn , ξ) ; 2) for (x, ξ) ∈ Drr , f (x, ξ, tn ) =
(
fb0 (x − ξtn , ξ) fb0 (−x + ξtn , −ξ)
(4.8)
x > ξtn ; x < ξtn
(4.9)
2D
10
X. WEN AND S. JIN
3) for (x, ξ) ∈ Dr+ , µ√ ¶ p p ξ 2 −2D 2 2 fb0 x − ξ − 2D tn , ξ − 2D ξ f (x, ξ, tn ) = b f0 (x − ξtn , ξ)
0 < x < ξtn
;
x > ξtn
4) for (x, ξ) ∈ Dl− , ¶ µ √ p p 2 fb0 − ξ +2D x + ξ 2 + 2D tn , − ξ 2 + 2D ξ f (x, ξ, tn ) = b f0 (x − ξtn , ξ)
(4.10)
ξtn < x < 0
.
x < ξtn (4.11)
4.3. The numerical solution. Denote µj = λtx |ξj | ,
1 ≤ j ≤ M.
(4.12)
Under the CFL condition (2.6), µj < 1 for 1 ≤ j ≤ M . Since Vx (x) = 0 except at x = xm+1/2 , with the boundary conditions (4.1), (4.2), Scheme I on DC is given by: 1) if 0 < ξj < ξM + 12 , i 6= m + 1, n+1 n n gij = (1 − µj )gij + µj gi−1,j ;
(4.13)
2) if ξIb − 12 < ξj < 0, i < m or ξ 12 < ξj < 0, i > m, n+1 n n gij = (1 − µj )gij + µj gi+1,j ;
3) if
√
(4.14)
2D < ξj < ξM + 12 ,
³ ´ n+1 n n n gm+1,j = (1 − µj )gm+1,j + µj cj,1 gm,d + c g ; j,2 m,d +1 j j
4) if 0 < ξj
0, x ∈ R, fb0 (x, ξ)dξ
(5.3) x 0, x ∈ R, ∂t ∂x ( ξJj,p x0 Z ξ 1 j+ 1 2 ujr,+,p |t=0 = f j (x, ξ)dξ, ∆ξ ξj− 1 r,+,p 2 ( b0 (x, ξ − ξj + ξJ ) f x 0, x ∈ R, ∂t ( ¯ ∂x¯ ¯ξ J ¯ x0 Z ξ 1 j+ 1 2 j f j (x, ξ)dξ, ul,−,p |t=0 = ∆ξ ξj− 1 l,−,p 2 ( b f0 (−x, ξ − ξj + ξJj,p ) j fl,−,p (x, ξ) = fb0 (−x, ξ)
(5.74) (5.75) (5.76) x ξ j tn
.
x ≤ ξj t n
Then ¯ ¯ 7 ¯ −,j,p ¯ (x) − x¯ ≤ ∆ξtn ¯ Hξ 2
(5.82)
for x ≤ 0, ξj− 12 ≤ ξ ≤ ξj+ 21 , Ib ≤ j ≤ I0 , p = 1, 2. This Lemma can be proved using similar technique for proving Lemma 5.3 and we omit the detailed proof. Lemma 5.7. Define functions Z 0 ¯ ¯ ¯ ej ¯ j Gp (ξ) = (5.83) ¯fl,−,p (x, ξ, tn ) − fbl,− (x, ξ, tn )¯ dx, x1
2
for ξj− 12 ≤ ξ ≤ ξj+ 12 , Ib ≤ j ≤ I0 , p = 1, 2, where j fel,−,p (x, ξ, t) =
fbl,− (x, ξ, t) =
³ ´ p fb0 ξJj,p x − ξJj,p t, − ξ 2 + 2D ξj
x > ξj t
b f0 (x − ξj t, ξ) x ≤ ξj t µ √ ¶ p p 2 fb0 − ξ +2D x + ξ 2 + 2D t, − ξ 2 + 2D ξ
,
b f0 (x − ξt, ξ)
(5.84) x > ξt
.
(5.85)
x ≤ ξt
Then ³ ´ √ Gjp (ξ) ≤ 2A + 2DB 7∆ξtn . Proof. From definition (5.85) and conditions (4.3), (4.4) one obtains √ kfbl,− (·, ξ, t)kBV (R) ≤ 2A + 2DB.
(5.86)
(5.87)
We use the fact that ³ ´ j fel,−,p (x, ξ, tn ) = fbl,− Hξ−,j,p (x), ξ, tn ,
x ∈ R, p = 1, 2.
(5.88)
From (5.83), (5.88) and (5.87), applying Lemma 5.6 and Lemma A.1 in the Appendix gives (5.86).
22
X. WEN AND S. JIN
We now give estimates for El− . In comparison with schemes (4.14), (4.17), (4.20), (4.21), (4.23) and (5.79), (5.80), (5.76), (5.77), (4.7), one can check that for Ib ≤ j ≤ I0 g(x, ξ, tn ) =
2 X
j cj,p vl,−,p (−x, tn ),
p=1
a.e. for ξj− 12 < ξ < ξj+ 21 , x 12 < x < 0.
(5.89)
Together with (5.73), (5.89) and the expression (4.11) of f (x, ξ, tn ) in Dl− , one obtains El− ≤ El−,1 + El−,2 + El−,3 ,
(5.90)
where El−,1
=
I0 Z X
=
ξj+ 1
Z
2
ξj− 1 2
I0 Z X j=Ib
by
2
I0 Z X j=Ib
El−,3 =
Z
2
ξj− 1
j=Ib
El−,2
ξj+ 1
ξj+ 1
2
ξj− 1 2
Z
2 X
0
x 1 p=1 2 2 X
0
x 1 p=1 2 2 X
0
x 1 p=1 2
¯ ¯ ¯ j ¯ cj,p ¯vl,−,p (−x, tn ) − ujl,−,p (−x, tn )¯ dxdξ,
(5.91)
¯ ¯ ¯ ¯ j cj,p ¯ujl,−,p (−x, tn ) − fel,−,p (x, ξ, tn )¯ dxdξ,
(5.92)
¯ ¯ ¯ j ¯ cj,p ¯fel,−,p (x, ξ, tn ) − fbl,− (x, ξ, tn )¯ dxdξ.
(5.93)
j Applying Theorem 3.1, the l1 -error between vl,−,p and ujl,−,p , p = 1, 2 are given
"r j kvl,−,p (·, tn )−ujl,−,p (·, tn )kl1 (R)
≤
kujl,−,p (·, 0)kBV (R)
# ¯ ¯ ¯ ¯¶ µ ¯ ξJj,p ¯ 2tn+1 ¯ξJj,p ¯ √ ¯ ∆x . p ∆x + 2 + ¯¯ e ξj ¯ |ξj | (5.94)
Similar to the derivation of (5.64), one has kujl,−,p (·, 0)kBV (R) ≤ 2A +
√
2DB.
(5.95)
Applying (5.94) and (5.95) to (5.91), by summation of j and using the inequalities (5.65) and (5.70) one obtains El−,1
´ r 2t
¸ 14 ¯ ¯ · ³ ´2 √ ¯ ¯ 1 1 ξ ξ − 2D ∆x ¯ 2¯ 2 e ³ ´¯ ¯ µ 1 ¶ √ ¯ ¯ + 2A + 2DB ¯ξ 12 ¯ ln ∆x + O(∆x). ∆x
³
√
≤ 4A + 2 2DB
n+1
(5.96)
Using Lemma 5.5, (5.67) and condition (4.4), from the definitions (5.76), (5.77), (5.78) and (5.84), one has ¯ ¯ ¯ j ¯ j ¯ul,−,p (−x, tn ) − fel,−,p (x, ξ, tn )¯ ≤ 2B∆ξ, for x ≤ 0, ξj− 21 ≤ ξ ≤ ξj+ 12 . Thus one has ¯ ¯¯ ¯ ¯ ¯¯ ¯ El−,2 ≤ ¯ξIb − 12 ¯ ¯x 12 ¯ 2B∆ξ.
(5.97)
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
23
Applying Lemma 5.7 to (5.93) yields ¯ ³ ´¯ √ ¯ ¯ El−,3 ≤ 7 2A + 2DB ¯ξIb − 12 ¯ tn ∆ξ.
(5.98)
Combining (5.90), (5.96), (5.97), (5.98) leads to El−
´ r 2t
¸ 14 ¯ ¯ ·³ ´2 √ ¯ ¯ ≤ 4A + 2 2DB ∆x ¯ξ 12 ¯ ξ 12 − 2D e µ ¶ ³ ´¯ ¯ √ 1 ¯ ¯ + 2A + 2DB ¯ξ 12 ¯ ln ∆x + O(∆x). ∆x ³
√
n+1
(5.99)
Finally combining (5.1), (5.20), (5.21), (5.33), (5.72), (5.99), (5.12) completes the proof for Theorem 4.1. Remark The error terms derived in (4.25) include the halfth order terms ¡ 1 5.1. ¢ and the ln ∆x ∆x terms. Thus the leading error terms are of halfth order and Scheme I has a halfth order convergence rate. Note that this is sharp, since even for the discontinuous solution to linear hyperbolic equation with a smooth coefficient, one cannot expect a better convergence order in the l1 -norm [24]. 6. Conclusion. In this paper we derived the l1 -error estimates for a Hamiltonianpreserving scheme, developed in [11], for the Liouville equation with a piecewise constant potential in one space dimension. The Hamiltonian-preserving scheme is designed by incorporating into the numerical fluxes the particle behavior–transmission and reflection– at the potential barrier. We proved that, with the Dirichlet incoming boundary conditions and for a class of bounded initial data, the numerical solution by the Hamiltonian-preserving scheme converges in l1 -norm to the solution of the Liouville equation–defined by using the interface condition that accounts for particle transmission and reflection. The initial data conditions can be satisfied by applying the decomposition technique proposed in [10] for solving the Liouville equation with measure-valued initial data arising in the semiclassical limit of the linear Schr¨odinger equation. The strategy for the error analysis in this paper is to apply the l1 -error estimates established in [30, 28] for the immersed interface upwind scheme to the linear advection equations with piecewise constant coefficients. This is a natural approach since the solution of the Liouville equation with a step function potential satisfies linear advection equations with piecewise constant coefficients on the bicharacteristics. To apply this strategy, we split the computational domain into several parts. In each subdomain, we introduced a number of linear advection equations with step function coefficients to which the immersed interface upwind schemes yield the same numerical solution as the Hamiltonian-preserving scheme for the Liouville equation with the step function potential. Then the l1 -error estimates for the Hamiltonian-preserving scheme were derived by applying the l1 -error estimates for the immersed interface upwind schemes and comparing the exact solutions of the linear advection equations and that of the Liouville equation in each of these subdomains. The condition on the initial data were used in deriving these estimates. As a result, we obtained the half order l1 -error bound with explicit coefficients for the Hamiltonian-preserving scheme. The problem under study has important applications in the computation of the semiclassical limit of the linear Schr¨odinger equation through barriers, and more generally, in the computation of high frequency waves through interfaces.
24
X. WEN AND S. JIN
Appendix In this Appendix we prove a property of BV functions on R. Lemma A.1. Let f (x) be a BV function on R, H(x) be a function on [a, b] satisfying |H(x) − x| ≤ HC ,
for x ∈ [a, b],
where HC is a positive constant. Then kf (·) − f (H(·))kL1 ([a,b]) ≤ 2HC kf kBV (R) .
(A.1)
Proof. ∀x ∈ [a, b], one has |f (x) − f (H(x))| ≤ kf kBV ([x−HC ,x+HC ]) . Since f (x) has a bounded variation on R, f 0 (x) exists a.e., and Z kf kBV ([x−HC ,x+HC ]) =
x+HC
|f 0 (y)| dy.
x−HC
Thus Z kf (·) − f (H(·))kL1 ([a,b]) ≤
b
Z
a
x+HC
|f 0 (y)| dydx.
(A.2)
x−HC
Introduce the function ( g(x, y) =
|f 0 (y)|
if |x − y| ≤ HC
0
else
.
Then g(x, y) ∈ L1 ([a, b] × [a − HC , b + HC ]). Applying Fubini Theorem gives Z Z
b a
Z
x+HC
Z |f 0 (y)| dydx =
x−HC b+HC Z b
=
b
Z
b+HC
g(x, y)dydx a a−HC Z b+HC Z y+HC
g(x, y)dxdy ≤ a−HC
Z
a b+HC
= 2HC a−HC
a−HC
|f 0 (y)| dxdy
y−HC
|f 0 (y)| dy ≤ 2HC kf kBV (R) .
(A.3)
Combining (A.2) and (A.3) gives (A.1). REFERENCES e, and G.D.V. Gowda, Godunov-type methods for conservation laws with [1] Adimurthi, J. Jaffr´ a flux function discontinuous in space, SIAM J. Numer. Anal. 42(1), 179-208, 2004. [2] L. Ambrosio, Transport equation and Cauchy problem for BV vector fields, Invent. Math., 158, 227-260, 2004. [3] F. Bouchut and F. James, One-dimensional transport equations with discontinuous coefficients. Nonlinear Anal. 32 (1998), no. 7, 891–933. [4] G.M. Coclite and N.H. Risebro, Conservation laws with time dependent discontinuous coefficients, SIAM J. Math. Anal. 36(4), 1293-1309, 2005.
L1 ERROR ESTIMATES FOR LIOUVILLE EQUATION
25
[5] R.J. DiPerna and P.-L. Lions, Ordinary differential equations, transport theory, and Sobolev spaces, Invent. Math. 98, 511-547, 1989. [6] T. Gimse, Conservation laws with discontinuous flux functions, SIAM J. Math. Anal. 24(2), 279-289, 1993. [7] L. Gosse and F. James, Numerical approximations of one-dimensional linear conservation equations with discontinuous coefficients, Math. Comp. 69(231), 987-1015, 2000. [8] P. G´ erard, P.A. Markowich, N.J. Mauser and F. Poupaud, Homogenization limits and Wigner transforms, Comm. Pure Appl. Math., 50(1997), pp. 321-377. [9] T. Gimse and N.H. Risebro, Solution of the Cauchy problem for a conservation law with a discontinuous flux function, SIAM J. Math. Anal. 23, 635-648, 1992. [10] S. Jin, H.L. Liu, S. Osher and R. Tsai, Computing multivalued physical observables for the semiclassical limit of the Schrodinger equation, J. Comp. Phys. 205, 222-241, 2005. [11] S. Jin and X. Wen, The Hamiltonian-preserving schemes for the Liouville equation with discontinuous potentials, Comm. Math. Sci. 3, 285-315, 2005. [12] S. Jin and X. Wen, A Hamiltonian-preserving scheme for the Liouville equation of geometrical optics with partial transmissions and reflections, SIAM J. Num. Anal. 44, 1801-1828, 2006. [13] K.H. Karlsen, C. Klingenberg and N.H. Risebro, A relaxation scheme for conservation laws with a discontinuous coefficient, Math. Comp. 73(247), 1235-1259,2004. [14] C. Klingenberg and N.H. Risebro, Convex conservation laws with discontinuous coefficients. Existence, uniqueness and asymptotic behavior, Comm. Partial Differential Equations 20(11-12), 1959-1990,1995. [15] C. Klingenberg and N.H. Risebro, Stability of a resonant system of conservation laws modeling polymer flow with gravitation, J. Differential Equations 170(2), 344-380,2001. [16] K.H. Karlsen, N.H. Risebro and J.D. Towers, Upwind difference approximations for degenerate parabolic convection-diffusion equations with a discontinuous coefficient, IMA J. Numer. Anal. 22(4), 623-664,2002. [17] K.H. Karlsen and J.D. Towers, Convergence of the Lax-Friedrichs scheme and stability for conservation laws with a discontinuous space-time dependent flux, Chinese Ann. Math. Ser. B. 25(3), 287-318, 2004. [18] P.L. Lions and T. Paul, Sur les measures de Wigner, Revista. Mat. Iberoamericana 9, 1993, 553-618. [19] L.W. Lin, B.J. Temple, and J.H. Wang, A comparison of convergence rates for Godunov’s method and Glimm’s method in resonant nonlinear systems of conservation laws, SIAM J. Numer. Anal. 32(3), 824-840,1995. [20] S. Mishra, Convergence of upwind finite difference schemes for a scalar conservation law with indefinite discontinuities in the flux function, SIAM J. Numer. Anal. 43(2), 559-577, 2005. [21] D.N. Ostrov, Viscosity solutions and convergence of monotone schemes for synthetic aperture radar shape-from-shading equations with discontinuous intensities, SIAM J. Appl. Math. 59(6), 2060-2085, 1999. [22] B. Perthame and C. Simeoni, A kinetic scheme for the Saint-Venant system with a source term, CALCOLO 38(4), 201-231 (2001). [23] F. Poupaud and M. Rascle, Measure solutions to the linear multi-dimensional transport equation with non-smooth coefficients, Comm. Partial Differential Equations 22 (1997), no. 1-2, 337–358. √ [24] T. Tang and Z.H. Teng, The sharpness of Kuznetsov’s O( ∆x) L1-error estimate for monotone difference schemes, Math. Comp. 64, 581-589, (1995). [25] B. Temple, Global solution of the Cauchy problem for a class of 2 × 2 nonstrictly hyperbolic conservation laws, Adv. in Appl. Math. 3(3), 335-375, 1982. [26] J.D. Towers, Convergence of a difference scheme for conservation laws with a discontinuous flux, SIAM J. Numer. Anal. 38(2), 681-698, 2000. [27] J.D. Towers, A difference scheme for conservation laws with a discontinuous flux - the nonconvex case, SIAM J. Numer. Anal. 39(4), 1197-1218, 2001. [28] X. Wen, Convergence of an immersed interface upwind scheme for linear advection equations with piecewise constant coefficients II: Some related binomial coefficients inequalities, preprint. [29] X. Wen and S. Jin, The L1 -stability of a Hamiltonian-preserving scheme for the Liouville equation with discontinuous potentials, preprint. [30] X. Wen and S. Jin, Convergence of an immersed interface upwind scheme for linear advection equations with piecewise constant coefficients I: L1 -error estimates, J. Comput. Math. 26(1), 1-22, 2008.
