THE LAX-FRIEDRICHS SCHEME FOR INTERACTION BETWEEN THE INVISCID BURGERS EQUATION AND MULTIPLE PARTICLES

. We propose a ﬁnite diﬀerence method based on the Lax-Friedrichs scheme for a model of interaction between multiple solid particles and an in- viscid ﬂuid. The single-particle version has been studied extensively during the past decade. The model studied here consists of the inviscid Burgers equa- tion with multiple nonconservative moving source terms that are singular and account for drag force interaction between the ﬂuid and the particles. Each particle trajectory satisﬁes a diﬀerential equation that ensures conservation of momentum of the entire system. To deal with the singular source terms we discretize a model that associates with each particle an advection PDE whose solution is a shifted Heaviside function. This alternative model is well known but has not previously been used in numerical methods. We propose a def-inition of entropy solution which directly generalizes the previously deﬁned single-particle notion of entropy solution. We prove convergence (along a subsequence) of the Lax-Friedrichs approximations, and also prove that if the set of times where the particle paths intersect has Lebesgue measure zero, then the limit is an entropy solution. We also propose a higher resolution version of the scheme, based on MUSCL processing, and present the results of numerical experiments.


(Communicated by Kenneth Karlsen)
Abstract. We propose a finite difference method based on the Lax-Friedrichs scheme for a model of interaction between multiple solid particles and an inviscid fluid. The single-particle version has been studied extensively during the past decade. The model studied here consists of the inviscid Burgers equation with multiple nonconservative moving source terms that are singular and account for drag force interaction between the fluid and the particles. Each particle trajectory satisfies a differential equation that ensures conservation of momentum of the entire system. To deal with the singular source terms we discretize a model that associates with each particle an advection PDE whose solution is a shifted Heaviside function. This alternative model is well known but has not previously been used in numerical methods. We propose a definition of entropy solution which directly generalizes the previously defined single-particle notion of entropy solution. We prove convergence (along a subsequence) of the Lax-Friedrichs approximations, and also prove that if the set of times where the particle paths intersect has Lebesgue measure zero, then the limit is an entropy solution. We also propose a higher resolution version of the scheme, based on MUSCL processing, and present the results of numerical experiments.
(1.1) Here f (u) = u 2 /2, and δ(x) denotes the Dirac delta measure concentrated at x = 0. The function u = u(x, t) models the velocity of the fluid, h k (t) models the location of the kth solid particle at time t, λ k > 0 is a drag coefficient associated with the 144 JOHN D. TOWERS kth particle, and m k > 0 is the mass of the kth particle. Study of the singleparticle version of (1.1) was initiated in [11], and has been the subject of a number of additional papers.
The fluid velocity is governed by the inviscid Burgers equation u t + f (u) x = 0, and the particle-fluid coupling is due to friction, more specifically the drag terms λ k (u − h k ) which appear in both the PDE and the ODEs in (1.1). Since there is no viscosity, the velocity u(x, t) admits entropy weak solutions, meaning that shock waves occur. This leads to complex interactions between the resulting shock wave and the particles. When multiple particles are present there are interesting features of the solutions that include particles drafting and passing by one another; see Figure 4 or Figure 5.
There are some difficulties associated with (1.1), in addition to the well-known ones associated with a nonlinear conservation law. The source terms on the right side of the first equation are nonconservative products of distributions; their meaning is not immediately clear. The differential equations appearing in the second line are coupled to the conservation law. Due to discontinuities in u the meaning of the right side of the DE's is also not readily apparent. There are related difficulties in designing practical numerical algorithms.
Notwithstanding these difficulties there has been much progress on the singleparticle version of (1.1). A notion of solution has been developed, well-posedness has been proven, and numerical algorithms have been designed whose approximations are known to converge to the unique solution. In this paper we focus on the multiple-particle problem, which has not been studied as thoroughly. We propose a notion of entropy solution suitable for multiple particles, present a Lax-Friedrichs difference scheme for the multiple-particle problem, and prove that the resulting approximations converge to an entropy solution. This is accomplished under the assumption that the particle paths do not intersect except possibly at a set of times whose Lebesgue measure is zero.
Reference [4] developed a unifying framework for the jump conditions that hold across a spatial flux discontinuity for a conservation law with discontinuous flux, using the theory of L 1 -dissipative (L 1 D) admissibility germs. The relevant L 1 D admissible germ for the problem discussed here is G(λ, c), which was identified in [7]. Definition 1.1 (the germ G(λ, c), [7]). The germ G(λ, c) is the subset of R 2 defined by G(λ, c) = (c, c)+{(a, b) ∈ R 2 |b = a−λ}∪{(a, b) ∈ R 2 |a ≥ 0, b ≤ 0, −λ ≤ a+b ≤ λ}. (1.2) Reference [6] gives a definition of entropy solution for the single-particle version of (1.1). The following is a direct generalization of that definition to the multipleparticle problem.
A function u is a solution of the first equation of (1.1) with initial data u 0 if u ∈ L ∞ (Π T ) ∩ C([0, T ]); L 1 loc (R)), if u is a Kružkov entropy solution in Π T \ Γ of the Burgers equation with initial data u 0 , and if for a.e. t ∈ (0, T ) the one-sided traces of u at each particle position satisfy (1.4) (iii) With the notation h = (h 1 , . . . , h K ), a pair (u, h) satisfying (i) and (ii) above is an entropy solution of the system (1.1).
Remark 1. Definition 1.2 requires strong one-sided traces u(h k (t) ± , t) along each path x = h k (t). Assuming that the particle trajectories do not intersect except possibly on a subset of (0, T ) having Lebesgue measure zero, the results of [13] guarantee existence of the required traces. This is due to the regularity of the paths x = h k (t) and the fact that u is a Kružkov entropy solution of the Burgers equation in Π T \ Γ.
Above we have used the notation BV(R) to denote the set of functions of bounded variation on R, i.e., those functions ρ : R → R for which where the sup extends over all M ≥ 1 and all partitions {ξ 0 < ξ 1 < . . . < ξ M } of R. Theorem 1.3 (Main theorem). The Lax-Friedrichs scheme described in Section 2 produces approximations that converge as the mesh size approaches zero, along a subsequence, to a pair (u, h) where u ∈ L ∞ (Π T ) ∩ C([0, T ]; L 1 loc (R)) and h k ∈ W 2,∞ ([0, T ]), k = 1, . . . , K. If the particle trajectories h k (t) do not intersect except possibly on a subset of (0, T ) having Lebesgue measure zero, then (u, h) is an entropy solution in the sense of Definition 1.2.
As mentioned above, there has been significant progress on the single-particle version of (1.1) [1,5,6,7,11]. The study of (1.1) started with reference [11]. Among other things the authors completely solved the Riemann problem for K = 1, and described the asymptotic behavior of solutions.
In reference [5], the authors introduce two finite volume methods for computing approximate solutions. One is a Glimm-like scheme, and the other is a well-balanced scheme that uses nonrectangular space-time cells near the interface. These methods employ random sampling for placing the particle at a mesh interface at each time step. The nonconservative source term is handled by using a certain well-balanced scheme that was analyzed in [7]. They avoid the use of a moving mesh, and also avoid the use of a Riemann solver for the full model. The case of multiple particles is addressed, and is handled via a splitting method.
Reference [14] presents a finite volume scheme that is based on the well-balanced scheme of [5,7], but uses an adaptive stencil as an alternative to using a moving grid. The multiple-particle case is handled by splitting. Reference [7] proves well-posedness for the problem This is a simplification of (1.1), but its analysis provides an important step in analyzing the full problem. As mentioned above the germ G(λ, c), which is required for the correct defintion of entropy solution, was identified in [7]. Reference [6] proves well-posedness of the model (1.1) for K = 1, assuming that the initial data is of bounded variation. Approximate solutions are generated via a wave-front tracking algorithm. Definition 1.2 is a direct generalization of the definition for K = 1 appearing in [6].
Reference [1] presents a class of finite volume schemes for (1.1) when K = 1. The schemes are similar to those in [5], but a moving grid is used, which keeps the particle located at a fixed cell boundary. The approximations are shown to converge to the unique entropy solution.
References [2] and [3] concern a generalized version of (1.1) (again, for K = 1), where the fluid is governed by the inviscid compressible Euler equations.
Reference [10] specifically deals with a multiple-particle problem. The authors prove well-posedness for a version of (1.1) where the particle paths h k (t) are given, i.e., the second equation of (1.1) does not appear.
Let H(·) denote the Heaviside function, i.e., the characteristic function of [0, ∞). The system (1.1) has the following equivalent formulation [5,11]: (1.6) Although the splitting approach for multiple particles used in [5] and [14] gives good numerical results, extending the convergence analysis from the single-particle to the multiple-particle problem seems difficult. Various bounds required for convergence are not preserved by the splitting steps. The numerical schemes in those papers are based on the model (1.1). In this paper we instead discretize (1.6), using Lax-Friedrichs differencing for each of the PDEs. The advantage of this approach is that the case of multiple particles is accommodated without splitting. This makes it possible to obtain a number of estimates which taken together give a convergence proof for the multiple-particle model. On the other hand, while the schemes of [1], [5], and [14] give very sharply resolved shocks at the particle locations, our Lax-Friedrichs method results in a substantial amount of smearing. With this in mind, we additionally propose a higher resolution version of the scheme, based on MUSCL processing. The rest of the paper is organized as follows. In Section 2 we describe the Lax-Friedrichs scheme mentioned above. In Section 3 we prove convergence, modulo a subsequence, of the approximations for u, as well as the approximations for h k . In Section 4 we prove convergence of the approximations for w k . In Section 5 we verify that the subsequential limit u is a Kružkov entropy solution in Π T \ Γ and satisfies the jump condition (1.3). In Section 6 we prove that the limit h k satisfies the differential equation (1.4). Section 6 concludes with the proof of Theorem 1.3. Section 7 describes the MUSCL processing mentioned above. Section 8 presents the results of some numerical experiments.
2. The Lax-Friedrichs scheme applied to (1.6). We use a uniform spatial mesh size ∆x, and temporal step size ∆t. Define . Let χ j (x) denote the characteristic function of I j , and χ n (t) the characteristic function of I n We denote by U n j the finite difference approximation of u(x j , t n ), U n j ≈ u(x j , t n ). Similarly W n k,j ≈ w k (x j , t n ). Let {Q n j } be a griddefined function such as {U n j } or {W n k,j }. We will use the following notational abbreviations: . We extend {U n j } and {W n k,j } from grid-defined functions to functions defined on all of Π T via where c n k ≈ h k (t n ) and h n k ≈ h k (t n ), with the initialization (h 0 k , c 0 k ) = (h k,0 , v k,0 ) Let µ = ∆t/∆x. The algorithm that we propose discretizes the first two equations of (1.6) via the Lax-Friedrichs scheme, the third equation using Euler's method: where q is a parameter. For our purposes q ∈ (0, 1/2]. The numerical fluxes in (2.7) result by applying the Lax-Friedrichs flux [12] to f (u) = u 2 /2 and g n k (w) = c n k w. Remark 2. The scheme (2.6) preserves solutions where the fluid velocity and particle velocities are equal to the same constant: Remark 3. Some explanation of the third equation of (2.6) is in order. Based on the third equation of (1.6), the third equation of (2.6) should be (approximately) equivalent to , a delta function concentrated at x = h k (t n ). In particular, we expect (1/2) j∈Z W n k,j+1 − W n k,j−1 ≈ 1 (in fact this holds with "≈" replaced by "="; this follows from (3.5) of Lemma 3.1), and so we can write the third equation of (2.6) in the form Thus, by defining Clearly there are other, possibly simpler, methods of discretizating the third equation of (1.6).
The reason for choosing this particular approximation is to ensure the discrete conservation of momentum property discussed below.
From the first two equations of (1.1) it follows that, at least formally, the total momentum of the system is conserved: The scheme (2.6) enforces a discrete version of (2.8).
Proposition 1. Assume that there is a 0 < J ∈ Z such that U n j = 0 for |j| > J, and that U n ∞ < ∞. Define the discrete momentum: The discrete momentum is conserved: Proof. Multiplying by ∆x and summing the first equation of (2.6) over j ∈ Z gives Multiplying the third equation of (2.6) by m k and then summing over k gives The proof is completed by adding (2.10) and (2.11). Define Lemma 2.1. Z n j satisfies the following (equivalent) evolution equations: (2.14) (2.15) Evidently (2.13) is a discretization of (2.15).
For the proof of (2.13), we start from the second equality of (2.18) and substitutê The identity (2.13) now follows directly from (2.19).
3. Convergence of u ∆ and h ∆ k . Let ∆ = (∆x, ∆t). For our convergence analysis we will assume that ∆ → 0 with µ fixed, and satisfying the following CFL condition: Additionally we assume that which will be satisfied automatically for ∆ sufficiently small.
The following properties hold: Proof. The proof is by induction on n. Clearly all of (3.3), (3.4), (3.5), and (3.6) hold at n = 0. Assume that those assertions hold at time step n. From (3.1) and the induction hypothesis it follows that To prove that (3.3) holds at time step n + 1 we rewrite (2.14) using incremental coefficients: Invoking the induction hypothesis then completes the proof of (3.3) for n + 1.
Next we prove that (3.5) holds for n+1. We rewrite the second equation of (2.6): is a convex combination of W n k,j−1 , W n k,j , W n k,j+1 , implying that W n+1 k,j ∈ [0, 1] after invoking the induction hypothesis. By differencing (3.11) we get Invoking the induction hypothesis again yields ∆ + W n+1 k,j ≥ 0. Finally, summing (3.13) over j and then applying the induction hypothesis yields j∈Z ∆ + W n+1 k,j = 1. To prove (3.4) holds at n+1, we employ the result of the previous two paragraphs. Recalling (2.12), the proven bound on Z n+1 j is equivalent to It is readily verified that u 0 min ≤ z 0 min and z 0 max ≤ u 0 max + K k=1 λ k . Replacing z 0 min and z 0 max in (3.14), the result is Recalling that λ k > 0 and W n+1 k,j ∈ [0, 1], it is clear that (3.4) holds.
To verify that (3.6) holds for n + 1, we start with the third formula of (2.6), from which it is evident that (3.16) The induction hypothesis yields j∈Z W n k,j+1 − W n k,j−1 = 2, and so after taking absolute values, and applying (3.2), equation (3.16) becomes from which the desired inequality follows readily.
Lemma 3.2. U n j and Z n j satisfy spatial variation bounds: Proof. We claim that the scheme is a so-called Total Variation Decreasing (TVD) scheme with respect to the variable Z n j , i.e., To prove the claim we use (3.8). We have shown that C n j+1/2 , D n j+1/2 ≥ 0. It suffices by a standard result [12, p. 116] to show that C n j+1/2 + D n j+1/2 ≤ 1. Using (3.9) we find that (3.21) Here we have used (3.1) to get the last inequality. The desired bound then results by recalling that q ≤ 1/2. Then by induction it follows from (3.18) that It is readily verified using (2.12) that where the constant B is independent of ∆.
Proof. Rearranging the first equation of (2.6), and using (2.17) to rewrite After taking absolute values, applying the triangle inequality, then using the bounds on c n k andÛ n j provided by Lemma 3.1, we sum over j ∈ Z. The result is where B 1 and B 2 are ∆-independent constants. The proof is completed by invoking Lemma 3.2, along with the observation that j∈Z W n k,j+1 − W n k,j−1 = 2, which follows from (3.5).
Lemma 3.4. The particle velocity approximations satisfy the following bound: Proof. We start with the third formula of (2.6). Subtracting c n k from both sides, taking absolute values, and then using the triangle inequality, the result is (3.28) The proof of (3.27) is completed using (3.4) and (3.6).
Lemma 3.5. The approximations u ∆ converge boundedly a.e. and in L 1 loc (Π T ) as ∆ → 0, along a subsequence, to some u ∈ L ∞ (Π T ) ∩ C([0, T ]; L 1 loc (R)). For each k ∈ {1, . . . , K} the sequence h ∆ k converges (along the same subsequence) in , and c ∆ k converges (also along the same subsequence) to h k in L 1 loc ((0, T )). Proof. The proof is a standard argument (e.g., the proof of Proposition 2.4 of [1]) using Lemmas 3.1, 3.2, and 3.3 for the u portion, and Lemmas 3.1 and 3.4 for the h k portion.
Remark 6. In Sections 5 and 6 we will assume that the particle trajectories do not intersect except possibly on a subset of (0, T ) having Lebesgue measure zero. The convergence result above holds without any assumptions about particle path intersections.
In what follows (u, h) refers to a fixed subsequential limit of the type whose existence is guaranteed by Lemma 3.5. When taking the limit as ∆ → 0 it is understood to be along this fixed subsequence.
4. Convergence of w ∆ k and z ∆ . Lemma 4.1. W n k,j satisfies a spatial variation bound and a time continuity estimate for each k ∈ {1, . . . , K}: Proof. The first part of (4.1) is evident from (3.5). For the second part of (4.1), we write (3.11) in the form Taking absolute values, and recalling from the proof of Lemma 3.1 that α n k , β n k ∈ [0, 1] yields Then summing over j ∈ Z and using j∈Z ∆ + W n k,j = 1, α n k + β n k ≤ 1/2, gives the second part of (4.1) boundedly a.e. and in L 1 loc (Π T ) for each k ∈ {1, . . . , K}.
Proof. Lemma 4.1 along with W n k,j ∈ [0, 1] (Lemma 3.1) guarantees that w ∆ k converges along a subsequence in L 1 loc (R + × R) and boundedly a.e. to some w k ∈ L ∞ (Π T ) ∩ C([0, T ]; L 1 loc (R)). A standard Lax-Wendroff calculation [9] proves that w k is a weak solution of One such weak solution is w k (x, t) = H(x − h k (t)). We will show that this is the only weak solution and the proof will be complete. Assume that w k andw k are both weak solutions of (4.4). This implies that for every φ ∈ C ∞ 0 (R × [0, T ]), It is readily verified that φ t + h k (t)φ x = ψ, φ(·, T ) = 0. Substituting into (4.5), we have Since (4.7) holds for any ψ ∈ C ∞ 0 (R × [0, T ]), we conclude that w =w a.e.
The following lemma is a direct consequence of (2.12), Lemma 3.5, and Lemma 4.2. t) boundedly a.e. and in L 1 loc (Π T ).

5.
Jump and entropy conditions for u. In this section we verify that the subsequential limit u is a Kružkov entropy solution in Π T \ Γ and satisfies the jump condition (1.3).
Since each of the particle paths t → h k (t) is continuous, each F i,j is closed, and thus F is also a closed subset of (0, T ). There are Clearly ∂G/∂U n j ≥ 0 since q ≤ 1/2. For ∂U n+1 j /∂U n j±1 , λ kŴ n k,j .

(5.5)
It is readily verified that each of these partial derivatives is nonnegative using (3.1) and the fact thatŴ n k,j ∈ [0, 1]. The following lemma is a straightforward consequence of (3.5) and Lemma 4.2.
S n j has the following properties:

7)
and as boundedly a.e. and in L 1 loc (Π T ). Lemma 5.3. The following discrete entropy inequalities hold for all κ ∈ [z 0 min , z 0 max ]: Proof. Writing (2.13) in the form Z n+1 j = P (Z n j+1 , Z n j , Z n j−1 ), it is readily apparent that P (κ, κ, κ) = κ. Using this observation the proof is a standard calculation [8,9], using the fact that P is a nondecreasing function of all three arguments (Lemma 5.1).
Proof. We start with the first inequality in (5.8), and use the identity Since ∆ ± S n j ≥ 0, we have (Z n j+1 ∨ κ)∆ + S n j ≥ κ∆ + S n j , (Z n j−1 ∨ κ)∆ − S n j ≥ κ∆ − S n j , (5.11) and so we can replace (5.10) by Following the proof of the Lax-Wendroff theorem [9], let φ be a nonnegative test function with φ(x, 0) = 0, and φ n j := φ(x j , t n ). We multiply (5.12) by φ n j ∆x, and then sum over j ∈ Z, n ≥ 0. After summation by parts the result is (5.14) After simplifying the last integral the result is A similar calculation starting from the second inequality of (5.8) yields (5.16) Recalling Assumption 5.1 and Remark 7, fix an interval I m := (a m , b m ) ⊆ (0, T ) where there are no path intersections, and fix a particle path, indexed by k. For this calculation we will use the abbreviations z ± (t) = z(h k (t) ± , t) and c k (t) = h k (t). The ordering of the particles does not change in I m , so we can assume that the particles are labeled so that where γ k = l<k λ l , and we have abbreviated z ± = z ± (t), c k = c k (t). Another such test function calculation, this time with (5.16) results in Continuing with the abbreviation z ± = z ± (t), c k = c k (t), for a.e. t ∈ I m we have wherec k = c k + γ k . Repeating this calculation with z + ≤ κ ≤ z − , we find that Plugging κ = z − into the first inequality of (5.22) and then into the first inequality of (5.23), and recalling f (z) = z 2 /2, yields The second inequality of (5.22) (for z − < z + ) or the second inequality of (5.23) (for z − > z + ) implies that in either case Substituting κ = z + into the second inequalities of (5.22) and (5.23) yields Finally, with > 0, we substitute κ = z + − into the first inequality of (5.22), and κ = z + + into the first inequality of (5.23). Sending ↓ 0 results in Thus either z + = z − or all of (5.24), (5.25), (5.26), (5.27) hold. Let u ± = u(h k (t) ± , t). Substituting z − = u − + γ k , z + = u + + γ k + λ k into these relationships we have shown that either Recalling Definition 1.1, and that c k = h k (t), it is evident from (5.28), (5.29) that Lemma 5.5. The following discrete entropy inequality holds for each κ ∈ R: We write the first equation of (2.6) in the form λ k c n k −Û n j W n k,j+1 − W n k,j−1 .

(5.33)
Invoking the monotonicity of G (Lemma 5.1), a standard calculation [8,9] yields for κ ∈ U. Substituting V n+1 j = U n+1 j − Q n j , and using the triangle inequality yields (5.31), assuming κ ∈ U. Now take the case where κ / ∈ U, say κ < u 0 min − K k=1 λ k . In that case (5.31) reduces to U n+1 which, recalling the first equation of (2.6), is clearly satisfied. The case where κ > u 0 max + K k=1 λ k is handled similarly. Lemma 5.6. The limit u is a Kružkov entropy solution in Π T \ Γ of the Burgers equation with initial data u 0 .
for every κ ∈ R and every nonnegative test The proof is based on the discrete entropy inequality (5.31). Due to the bounds on U n j and c n k (Lemma 3.1), we have for some B > 0 which independent of ∆, (5.38) Multiplying by φ n j = φ(x j , t n ) and then summing by parts we find that The proof is finished by observing that R H(x−h k (t))φ x dx = 0, since φ(h k (t), t) = 0.
6. Differential equation for h k and proof of the main theorem. In this section we prove that the limit h k satisfies the differential equation (1.4). This section also contains the proof of Theorem 1.3. Assumption 5.1 (restriction on particle intersections) remains in effect in this section.
Lemma 6.1. The limit h k (t) satisfies the differential equation (1.4) for each k ∈ 1, . . . , K and a.e. t ∈ (0, T ). Also, (h k (0), h k (0)) = (h k,0 , v k,0 ). Proof. Fix a particle with index k, 1 ≤ k ≤ K. Let a n k = (c n+1 k − c n k )/∆t. The third equation of (2.6) yields where ψ δ is defined by (5.1). Let ξ(t) ∈ C ∞ 0 ((0, T )) and define ξ n = ξ(t n ). We re-write (6.1) in the form m k a n k = − Next we multiply by ξ n ∆t and sum over n : We solve for c n k −Û n j W n k,j+1 − W n k,j−1 in the first equation of (2.6), and substitute into the first sum on the right side of (6.3). The result is Summing the left side of (6.5) by parts, we find that Letting ∆ ↓ 0 in (6.6), and using c ∆ k → h k , the result is and for S 1 , summation by parts followed by sending ∆ → 0 yields We next estimate S 2 . Fix l = k. It suffices to estimate S 2,l , where Since c n l andÛ n j are bounded (Lemma 3.1), and W n l,j+1 − W n l,j−1 ≥ 0, ψ n j ≥ 0, where B is some positive number independent of δ and ∆. Summation by parts yields j∈Z W n l,j+1 − W n l,j−1 ψ n j = j∈Z W n l,j+1 ψ n j+1 − W n l,j−1 ψ n j−1 − j∈Z W n l,j+1 + W n l,j ∆ + ψ n j . (6.11) The first sum on the right is telescoping and is equal to zero. Thus, referring back to (6.10) we have (6.14) Substituting into (6.13) yields the desired estimate of S 2,l : We claim that S 3 → 0. Since c n k andÛ n j are bounded (Lemma 3.1), and ψ n j ≤ 1, W n k,j+1 − W n k,j−1 ≥ 0, where B is some positive number independent of the mesh size ∆. Using the formula (6.11) with 1 − ψ n j replacing ψ n j , (6.17) In the second term on the right side we have used ∆ + (1 − ψ n j ) = −∆ + ψ n j . The first sum on the right is telescoping and is equal to 2. Thus, referring back to (6.16) we have (W k,j+1 + W k,j ) ∆ + ψ n j /∆x.

Proof of the main theorem.
Proof. Lemma 3.5 provides the convergence portion of Theorem 1.3. That the limit (u, h) is an entropy solution results from Lemmas 5.4, 5.6, and 6.1.
Remark 8. For the single-particle case, Theorem 8 of [6] states that Definition 1.2 is sufficient for uniqueness. Thus if K = 1, the Lax-Friedrichs approximations (u ∆ , h ∆ 1 ) converge to the unique entropy solution, and convergence is along the entire computed sequence, not just a subsequence. 7. Improved resolution via MUSCL processing. It is possible to somewhat reduce the excessively diffusive nature of Lax-Friedrichs differencing without adding too much complexity by using the MUSCL approach. Our incorporation of MUSCL processing is standard [12]. Let M(·, ·) denote the minmod function: We replace the numerical fluxesf n j+1/2 ,ḡ n k,j+1/2 in (2.7) bỹ W n,± k,j = W n k,j ± 1 2 M(∆ + W n k,j , ∆ − W n k,j ).

(7.3)
We do not presently have any convergence results or even stability estimates for the resulting scheme with MUSCL processing incorporated. A moderate amount of numerical experience indicates that the algorithm produces approximations that converge to the same solution as the basic algorithm of Section 2. The exact solution is available for comparison, using the results of [11]. See Figure 1.
The approximations appear to improve when the mesh size is halved, as expected.
It is also apparent that the MUSCL scheme is more accurate than the basic one.
The sharp transition at x ≈ 0.8 is a shock that is collocated with the particle. With our Lax-Friedrichs scheme there is some smearing of the shock. We must rely on a very small mesh size, even with the MUSCL version, to obtain a very sharp transition. The schemes of [1], [5], and [14] resolve this type of shock ( i.e., the shock is collocated with the particle) with no smearing.  As in the previous example the exact solution is available via [11]. This example displays a spurious kink, see Figure 2, that appears in some cases where a particle's velocity h k (t) lies between u(h − k (t), t) and u(h + k (t), t). The kink is probably due to the large numerical viscosity of the Lax-Friedrichs scheme. The size of the kink diminishes, as expected, when the mesh shrinks. Also the MUSCL approximation has a smaller kink than the basic approximation.
Example 8.3. This is a two-particle example with z(x, t) = constant =ẑ. It is possible to explicitly solve this type of problem. With z(x, t) =ẑ, we have u(x, t) =ẑ − λ 1 H(x − h 1 (t)) − λ 2 H(x − h 2 (t)). Thus the problem reduces to determining the particle paths h 1 (t) and h 2 (t). This can be accomplished using the differential equations (1.4), which become Here where (8.5) Assume that the particle trajectories do not intersect except for a finite set of times τ ν with 0 < τ 1 < . . . < τ M < T . Define τ 0 = 0, τ M +1 = T , and let r k = λ k /m k . The solution of (8.3), (8.4), (8.5) can be expressed piecewise. For t ∈ (τ ν , τ ν+1 ) the solution is (1−exp(−r k (t−τ ν )))+ σ k r k (t−τ ν ).  Figure 3 it appears that the MUSCL scheme is more accurate than the basic scheme, as expected. We also see that the discrete L 1 error in u decreases as we decrease the mesh size. Figure 4 shows the approximate and exact particle trajectories. At the level of discretization shown, the particle trajectories produced by the basic scheme do not quite agree with the exact trajectories. This discrepancy diminishes when the mesh size is decreased (not shown), but convergence is slow. For the MUSCL scheme the resolution is better; the exact and computed trajectories are not visually distinguishable at this level of grid refinement.
Example 8.4. This is another two-particle example. This time the particles are initially heading toward each other, and the fluid is initially at rest. Unlike the previous example the true solution is not known. In Figure 5 we show the particle trajectories at three levels of grid refinement. It appears that the particle trajectories are converging as the mesh size is refined. The MUSCL scheme is better able to resolve the fine details of the trajectory, especially after the first crossing of trajectories.
The initial fluid velocity is zero, u 0 (x) = 0. The other parameters of the problem are