A PARTICLE METHOD AND NUMERICAL STUDY OF A QUASILINEAR PARTIAL DIFFERENTIAL EQUATION

. We present a particle method for studying a quasilinear partial differential equation (PDE) in a class proposed for the regularization of the Hopf (inviscid Burger) equation via nonlinear dispersion-like terms. These are obtained in an advection equation by coupling the advecting ﬁeld to the advected one through a Helmholtz operator. Solutions of this PDE are “regularized” in the sense that the additional terms generated by the coupling prevent solution multivaluedness from occurring. We propose a particle algorithm to solve the quasilinear PDE. “Particles” in this algorithm travel along characteristic curves of the equation, and their positions and momenta determine the solution of the PDE. The algorithm follows the particle trajectories by integrating a pair of integro-diﬀerential equations that govern the evolution of particle positions and momenta. We introduce a fast summation algorithm that reduces the compu- tational cost from O ( N 2 ) to O ( N ), where N is the number of particles, and illustrate the relation between dynamics of the momentum-like characteristic variable and the behavior of the solution of the PDE.


1.
Introduction. Consider the quasilinear partial differential equation (PDE) We refer to the above equation as the Helmholtz-regularized Hopf equation. The equation is the b = 0 member of the b-family m t + um x + b u x m = 0, m = u − α 2 u xx , α > 0, (2) described in [7]. The present paper will focus only on equation (1), the b = 0 case.
With the definition of the one dimensional Helmholtz operator as and its Green's function, one can represent the solution of equation (1) u in terms of m as u(x, t) = 1 2α ∞ −∞ e −|x−y|/α m(y, t)dy.
The Leray-type [10] smoothed velocity u has been recently proposed for its role of regularizing the inviscid Burgers (Hopf) equation.
With nonzero viscosity, the Burgers equation is the simplest model combing the nonlinear propagation effect and the diffusive effect, while its inviscid version for ν = 0, also known as Hopf equation, is the simplest model of shock forming hyperbolic equation. In [11,12], attempts were made to compare regularization effects induced by the viscosity and the filtered velocity (4). It is well-known that one can study the Burgers equation to determine the physically relevant solutions of the Hopf equation. The regularized Burgers equation (1), written in the form u t + uu x = α 2 (u xxt + uu xxx ) , (7) is a quasilinear equation that consists of the Hopf equation plus O(α 2 ) nonlinear terms. An approach, similar to that for the Burgers equation, is proposed (see, e.g., [1]) to show that solutions of this equation converge strongly to physically relevant weak solutions of the Hopf equation (5) as α → 0, provided the initial data u(0, x) are in a suitable function space. Thus, equation (7) has been proposed as an alternative to Burgers equation (5) in this respect. Part of the interest in studying this particular class of regularizations stems from the fact that, unlike its Burger counterpart, it can be endowed with interesting mathematical structures [7] shared with the rest of the family (2).
It is worth noting that in the absence of either of the two terms on the right hand site of (7), the equation becomes Both of these are dispersive equations (as linearization around a constant solution readily shows) and have zero-dispersion limits. One can expect that high-frequency oscillation will occur when α is small. However, when the two third derivative terms, uu xxx and u xxt , appear together in equation (7), the delicate balance between the two terms will cancel the high-frequency oscillation expected to arise in passing to the zero-dispersion limits. For most numerical schemes, certain care is required to discretize the two third-order terms in order to achieve a higher-order accuracy, while maintaining the delicate balance of the two terms. A dispersion-relationpreserving scheme is developed in [6] that maintains the balance between the two dispersive terms and is suitable to study the long-time solution behavior of equation (7). The study in this paper, however, takes a different route by developing a particle method that avoids any discretization of the spatial derivative terms in equation (7). "Particles" in this algorithm travel along characteristic curves of the equation, and their positions and momenta determine the solution of the PDE. The algorithm follows the particle trajectories by integrating a pair of integro-differential equations that govern the evolution of particle positions and momenta. The particle method is a self-adaptive algorithm, since particles normally cluster in areas with sharp gradients. The particle method reveals dynamics of particle position and momenta, which sheds light on solution behaviors of the PDE we study. Such a special feature distinguishes the particle method from Eulerian methods.
The rest of the paper is organized as follows. In Section 2, we define the characteristic variables and derive the evolution (integro-differential) equations for the characteristic variables. We introduce a finite dimensional particle system in Section 3 by discretizing the integro-differential equations. We derive an implicit analytic solution for a special case of two-particle dynamics, the particle-antiparticle collision, which, although very simple, is remarkably informative on the trend of PDE solutions with more general initial conditions. We verify the particle method in Section 4 by comparing the numerical solution of the particle-antiparticle collision with the implicit exact solution. Then we present an example with an initial condition that is a Gaussian hump. We introduce a fast summation algorithm that reduces the computational cost from O(N 2 ) to O(N ), where N is the number of particles. Using these more general initial data, we verify the proposed method by comparing our numerical solution with the solution computed by a two-step iterative method. Finally, we illustrate the relation between dynamics of the momentum-like characteristic variable and solution behaviors of the PDE.
2. Characteristic variables. In this section we introduce a particle method for the regularized Burgers equation (1) following the one proposed in [3] (and studied extensively in, e.g., [4,5,6,8,9]) for a shallow water equation that shares a similar structure. By introducing the characteristics variable (curve) equation (1) is equivalent to where the material (total) derivative is defined as Equation (11) suggests that the variable m is a constant along the characteristic curve, or m(q, t) = m 0 (ξ).
Equations (11) and (4) imply thaṫ Defining an auxiliary variable taking a time derivative about the above equation, and using equation (14) yields the evolution equatioṅ Equations (14), (15), and (16) yield a pair of continuous evolution equations for p and q: After solving the above evolution equation, solutions of the PDE (1) u are found by Here the characteristic q(ξ, t) plays the role of positions conjugate to the momentumlike variable p(ξ, t). The initial conditions for p(ξ, t) and q(ξ, t) are The choice of initial condition for the position variable, dictated by the characteristics condition, implies q ξ (ξ, 0) = 1, so that the constraint, obtained from equation (15), is maintained at all times of existence of the solution (q(ξ, t), p(ξ, t)). Thus, the momentum variable p(ξ, t) could be eliminated from the system to obtain an evolution equation containing only the dependent variable q(ξ, t) and its first derivative with respect to the initial label ξ. Vanishing of this derivative generically corresponds to crossing of characteristics curves, with loss of uniqueness of solutions ξ(x, ·) to the equation x = q(ξ, ·) that defines the characteristic map. Constraint (20) can then be used to show that if the initial condition p(ξ, 0) is bounded, then q ξ (·, t) is bounded away from zero, thereby preventing characteristics from crossing, as long as p(·, t) does not have zeros. We argue, by contradiction, that for a finite ξ the derivative q ξ cannot vanish in any finite time, provided the initial data p(ξ, 0) = 0. Suppose that t = T is the first time when crossing of characteristics occurs. Let Ξ be the location from which this curve emanates at t = 0. That is, q ξ (Ξ, T ) = 0, and q ξ (ξ, t) = 0 for all finite ξ and 0 < t < T . For simplicity, suppose that p(ξ, 0) is sign definite and p(ξ, 0) > 0; definition (20) then implies that q ξ (ξ, t) > 0 for t ≤ T , as long as p(ξ, t) > 0 for t ≤ T . If we can show that p(ξ, t) > 0 for t ≤ T , then we have the contradiction. To show this, let the kernel in equation (4) be where α > 0. Since 0 < G(ξ, η) ≤ 1 2α , theq equation satisfies the inequalitẏ where By Gronwall's inequality, it follows that for all characteristic initial conditions ξ ∈ R, and times t ≤ T . We now argue that p(ξ, t) > 0 for 0 ≤ t ≤ T . For the time period 0 ≤ t ≤ T , the sign definiteness of p(·, t) and theṗ equation in (17) implẏ by theq equation in (17). Using Gronwall inequality again, we obtain Contradiction could be resolved by failure of the assumption q ξ < ∞ for t ≤ T : Suppose that there could exist a time T 1 < T and a pointξ such that q ξ (ξ, where p(ξ, 0) is sign definite, p(ξ, t) could change sign. Thus, q ξ could change sign. However, it is easy to show that p(ξ, t) < ∞, and hence q ξ (ξ, t) < ∞ by (20). In fact, similar to (25), for the time interval 0 ≤ t ≤ T , By Gronwall's inequality, we obtain for all ξ and 0 ≤ t ≤ T . This contradicts q ξ (ξ, T 1 ) = ∞.
3. Finite dimensional particle system. The integrals in system (17) can be approximated by their Riemann sums, thereby yielding discrete systems for "particles" with coordinates where ξ i = Ξ + ih for some real Ξ, step-size h > 0 and i = 1, · · · , N . The discretized version of the system results in the finite dimensional dynamical system of N particles, System (34) constitutes our particle method for solving the regularized equation (7). Solution of the PDE are found by which is often referred to as "peakon" superposition in the literature.
3.1. Two-particle dynamics. Substituting the scaled variables q * , p * , x * , t * , and u * into system (34) and equation (35), where yields (drop ' * ', herein and after) and The scaled system of evolution equations (37) corresponds to the case α = 1 in the PDE (7) It is worth noting that making use of spacial and temporal scalings we can find solutions of equation (7) from equation (39). That is, given a solution u(x, t) of (39), the solution of (7)ũ(x, t) can be found bỹ We now study the dynamics of two particles (N = 2) for the evolution equations (37). Suppose that two particles are initially well separated and have speeds c 1 and c 2 , with c 1 > c 2 , and c 1 > 0, so that they collide. Define the new variables which obey the equationṡ Consider the special case P ≡ 0 (p 1 = −p 2 at all times). Equation (42) becomeṡ Equation (43) implies dp dq Solving equation (44), we obtain for some constant C > 0. Substituting equation (45) into the theq equation in (43), we obtain dq Considering the case q < 0 for some t and letting u ≡ e −|q| , we integrate equation for some t 0 . Hence an implicit expression of the solution q is given by Note that equation (48) implies when q → 0 − , or Equation (50) suggests that, in the course of particle and antiparticle collision, q → 0 at the rate of O(t −2 ), when t → ∞. On the other hand, equation (45) implies that p → 0 at the rate of O(t −1 ), when t → ∞.
4. Numerical investigation. The particle algorithm based on system (34) can be readily implemented numerically with an appropriate ODE integrator. A simple test of its performance is provided by the two particle (or peakon) collision case itself.
In Figure 1, we solve the evolution equation (37) with the initial momenta p 1 = −1,  (45) and (48), within fourteen decimal digits (machine precision). We summarize the calculation in Table 1.  Table 1. Computed p i and q i in the example of perfectly antisymmetric particle-antiparticle collision.

4.2.
Gaussian disturbance initial data. Next, we validate the particle method on the case of the initial condition corresponding to a Gaussian disturbance where w is related to the initial width of the Gaussian disturbance. We choose w = 20 and α = 1. Figure 2(a) compares the evolution results by using the proposed particle method and the two-step iterative algorithm developed in [6] for the solutions at the final time t = 300 and t = 600. As demonstrated by the figures, to within graphical accuracy the two independent numerical methods are identical. The simulations are performed using particle number (grid cells) N = 10000 in the domain [−50, 150]. In Figure 2(b), we show particle trajectories of about 50 particles in the (x, t) plane. The formation of the sharp gradient (quasi-shock) in Figure 2(a) is associated with the particle clustering in Figure 2(b). The clustering process acts as a high particle-density front sweeping through the computational domain. This front corresponds to the trailing edge of a quasi-shock wave leaving behind a rarefaction wave. As a consequence of the particle clustering, we implement a redistribution algorithm to prevent two particles from occupying one location at any finite time: when two particles, with positions q i and q i+1 , are too close to be distinguished within machine precision, we replace them with one particle at the same location carrying a momentum equal to the sum of p i and p i+1 . After we carry out this replacement, we relabel the rest of the particles from the original i + 2, . . . , N/2 to i + 1, . . . , N/2 − 1. The end result is to reduce the dimension of the system of ordinary differential equations (ODEs) from 2N to 2N − 2 by combining any two clustering particles. Of course, this method is somewhat crude, as this process of replacement depletes the total number of particles. While it is not too difficult to implement a true redistribution that conserves particles, we leave this out for simplicity. The redistribution algorithm results in the following fast summation algorithm.  Figure 2. (a) Comparison of the evolution results between the particle method and the two-step iterative algorithm developed in [6] for the solutions at the final time t = 300 and t = 600. (b) Particle trajectories of about 50 particles in the (x, t) plane. The formation of the sharp gradient (quasi-shock) in (a) is associated with the particle clustering in (b).

4.3.
Fast summation algorithm. The redistribution algorithm prevents particle collision. Without particle collision, we strip the absolute value notation in the power of the exponential function, which in turn makes a recursion relation for evaluating the sums possible. With the help of this recursion formula, the total operations needed for performing the summation is reduced to O(N ) for the Nparticle system. Because particles do not collide, the particle method (34) has q i > q j if i > j and vice-versa if i < j. Hence equation (34) can be written aṡ Define new variables: Equation (53) then becomeṡ One can see that with a pre-computed f l and f r the number of operations needed for the Riemann sum is O(N ) for the N × N system of equations. Since the operations required for f l or f r are also growing as O(N ), the total number of operations is O(N ).
We now establish the recursion relation for f l and f r : Similarly, Because e −(qi−qi+1)/α leads to an exponential growth, for numerical stability the recursion relation for f r is better solved backward as 4.4. Dynamics of momentum-like variable. In the paper [6], we demonstrate that for smooth nearly square-wave initial data, where w is related to the width of the profile of the hyperbolic tangent function, the regularized Burger equation (7) exhibits markedly different behavior from the corresponding Hopf equation, for an arbitrary small α. The mechanism of such solution behavior is appreciated when looking at the finite dimensional dynamical system of the particle method. Figure 3(a) shows the velocity u of the nearly square-wave initial data with w = 0.1, while Figure 3(b) is the initial condition of the momentumlike variable p. The initial distribution of p exhibits an "antiparticle-particle" structure on the left and a "particle-antiparticle" structure on the right. We show in [6] that when the initial distribution of p is associated with the "antiparticle-particle" structure, the left region of the nearly square wave grows into a hump almost immediately after the initial time. The amplitude of the hump depends on both w and the parameter α. Our numerical experiments indicate that as w decreases or α increases, the amplitudes of the humps also become larger. The mechanism of such a growth is related to the magnitude of the second derivative of u, which increases when w decreases. Hence, for a fixed α, no matter how small α is, m = u − α 2 u xx departs significantly from u.
To investigate further the mechanism of such a growth, we consider the following initial data The initial p distributions of equations (60) and (61) are reminiscent of the left and right region of p in Figure 3(b), respectively. The initial condition of the momentum-like variable p. The initial distribution of p exhibits a "antiparticleparticle" structure on the left and a "particle-antiparticle" structure on the right.   Figure 4(a,b) show the evolutions of the velocity u and the momentum-like variable p, respectively, with the initial data (60). The simulations show that with the antiparticle-particle structure, the initial waves travel away from each other and the magnitudes of the peaks grow and then decay, while the magnitude of the pvariable continues to grow. The parameters α and w in the simulations are both set to be one. Since the initial p is equal to m 0 = u 0 − α 2 u 0xx , the growth of amplitude depends on the magnitude of the second derivative of u initially, which increases when w decreases, and on the parameter α. Hence when α increases, or w decreases, the initial growth of the antiparticle-particle waves becomes more significant. Such behavior distinguishes the regularization effect by the diffusion term νu xx in the Burgers equation (6) from that by the third-order derivative terms α(uu xxx + u xxt ) in the regularized Burgers equation (7). For the diffusion term, when ν increases, the amplitude of waves will be dissipated with time more noticeably, regardless of whether the initial condition is an antiparticle-particle wave or a particle-antiparticle wave.   Figure 4. The negative sign in the initial data (61) turns the anitparticle-particle wave structure upside down. Such a particleantiparticle structure behaves similarly to the particle-antiparticle collision described in Section 3.1. One can expect that both velocity and momentum approach zero as time approaches infinity. Figure 5's panels (a) and (b) show that both amplitudes of the velocity u and the momentum-like variable p decay with time.
The implication of the above examples is that for the nearly square-wave initial condition (59), in the limiting case when w → 0, the left edge of the wave (resembling to the Riemann problem u l < u r ) might grow indefinitely for any fixed α. Based on a rescaling of the independent spatial variable with w,x = x/w, so that the Helmholtz operator becomes 1 − α 2 /w 2 ∂ 2 x , one can see that a power law magnitude orderings between w and α might exist, for which the time scale when solutions of the regularized Burgers equation can be expected to behave like those of the Hopf equation can be made arbitrarily large. On the other hand, for the right edge of the nearly square-wave initial data (59), the particle-antiparticle example predicts that since m decays, this edge (resembling to the Riemann problem u r < u l as w → 0) does not exhibit any growth, and the solution appears to be smeared out around the corners, just like solutions of the Burgers equation.

5.
Conclusion. We introduce a particle algorithm for studying the regularized Burgers equation. The equation consisting of the Hopf equation and two thirdorder derivative terms is a model that regularizes the shock formation observed in the Hopf equation. The small parameter α in front of the third-order derivative terms is analog to the viscosity (diffusion coefficient) in the Burgers equation. We illustrate that dynamics of the momentum-like variable p in the particle method is related to solution behaviors of the PDE. When p has a particle-antiparticle structure (resembling to shock waves in the Hopf equation), solutions of the regularized Burgers equation behave like the particle-antiparticle collision. Namely, p approaches to zero as time approaches to infinity, and the wave amplitude decreases with time. On the contrary, when p has an antiparticle-particle structure (resembling to rarefaction waves in the Hopf equation), the wave amplitude first grows and then decays 1 . Such behavior is different from that of the Burgers equation, where viscosity dissipates both shock waves and rarefaction waves immediately.