{'fraction_non_alphanumeric': 0.05363296607290201, 'fraction_numerical': 0.016065621511430376, 'mean_word_length': 3.870921985815603, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equal-area method for scalar conservation laws

We study one-dimensional conservation law. We develop a simple numerical method for computing the unique entropy admissible weak solution to the initial problem. The method basis on the equal-area principle and gives the solution for given time directly.


Introduction
We study the one-dimensional scalar conservation law (1) u where u : R × [0, ∞) → R represents the conserved quantity (such as mass, density, momentum) while f : R → R is the flux. We equip the equation (1) with the initial condition (2) u(x, 0) = h(x), x ∈ R.
It is well known that classical solutions of (1)- (2), even for smooth initial conditions, do not always exist. Usually, the concept of solutions is generalized to the so-called weak solutions which are not unique and, in order to obtain the proper solution, one has to impose some entropy condition. There is a rich mathematical theory on this topic, see for example monographs [1,5,6,7,8].
We proceed more directly. In order to allow discontinuities we consider the integral form of the conservation law for all a < b and t > 0 for which u is continuous in points (a, t) and (b, t). It says that the area b a u(x, t) dx changes in time according to the flux at the boundaries. If both u and f are continuously differentiable, (3) implies (1). We search for the solutions of (3) inside a class of functions Υ that is defined by a characteristic surface associated to the problem (1)- (2). This functions may have jumps along some locally smooth curves in the (x, t)-plane (see Definition 2.2). In Section 2 we prove uniqueness of solutions to (3) within the class Υ.
Taking the integral over all R, (3) implies Hence the area under the graph of the solution does not change in time. Based on this observation we propose in Section 3 a simple method which we call an equal-area method. We show that, under suitable conditions, it yields the unique solution of (3) within the class Υ.
The idea of using the equal-area principle is not new. It has already been exploited in the classical textbook by Whitham [12,]. However, the method there is developed only analytically in a complicated way that is not suitable for explicit computations. Recently, Farjoun and Seibold [3] suggested a conservative particle method that uses equal-area principle. They use the Lagrangean approach, representing the solution as a cloud of particles which move with the flow. Particles carry function values and move according to their characteristic velocities. When the characteristic curves collide, the particles are merged in such a way that the total area under the function is conserved. So far we are not aware of any other numeric method using thr equal-area principle.
In Section 4 we describe a numerical algorithm we used to implement the method. We discuss numerical results obtained by our method in Section 5 and compare it to the finite volume method Clawpack [2] as well as to the newly suggested conservative particles method Particleclaw [9].
Our method performs well. It is exact and the shocks are formed when necessary, with high accuracy. In contrast to other known methods, where the solution at selected time is obtained by evolution from initial conditions, it gives the solution for any given time directly.

Characteristic surface and integral solutions
A well-known method for treating the initial value problem (1)-(2) is the method of characteristics. The characteristics are curves in the (x, t)-plane along which the function u is constant. In our case they are lines of the form (5) x ξ (t) = ξ + f (h(ξ)) t, ξ ∈ R, see e.g. [8,Sec. 2.2]. It is easy to see that, if a C 1 -solution u(x, t) of (1)-(2) exists, its graph in R 3 is given by We shall call Γ the characteristic surface to the problem (1)- (2). An example of such a surface is seen in Figure 1(a).
Out  We can form Γ a-priori, before investigating the solvability of (1)-(2), whereby it might happen that it represents a multivalued function which cannot be the solution of our problem (see Figure 1(b)). Indeed, this problem occurs whenever the characteristics collide in the (x, t)-plane. In this case the proper solution has jumps along some curves in the (x, t)-plane. Its graph, however, is still a subset of Γ. Starting with Γ, the solution can thus be obtained by finding the appropriate position of the jumps.
In the case when the initial function h is not continuous the above defined characteristic surface Γ is not connected. Before proceeding we shall hence modify the definition of characteristic surface to correct this. For some fixed t ≥ 0 we define plane transformation G t as Denote by γ 0 the graph of the initial function h together with vertical line segments joining discontinuities and which is continuous curve for all t ≥ 0, see Figure 2. Definition 2.1. The bounding characteristic surface to the problem (1)-(2) is defined as Note thatΓ contains the characteristic surface Γ and that the two surfaces agree whenever h is continuous.
We are now ready to define the appropriate domicile for our solutions.
Definition 2.2. Let f ∈ C 2 (R) and h is a piecewise C 1 -function with compact support. We say that u ∈ Υ = Υ(f, h) if the following conditions are satisfied: for all x ∈ R, (iii) the graph of u in R 3 is a subset of the bounding characteristic surfaceΓ defined by f and h, (iv) the boundary of the graph of u is a finite union of C 1 -curves, and the projections of these curves on (x, t)-plane intersect any line t = t 0 only finitely many times, Our aim is to find the solutions to the integral form of the conservation law (3) within the class Υ. First we demonstrate that every u ∈ Υ solves (1) on its areas of smoothness.
Proof. If the graph of u(D) lies on the original characteristic surface Γ, then u solves (1) on D by the method of characteristics (see e.g. [5, §3]). Therefore we may assume that the graph of u(D) is contained inΓ \ Γ.
We will show that the added vertical lines and their convolutions in time also yield a solution in a similar way as the original characteristics do. For any fixed ξ we parametrize the complemented surface inΓ \ Γ as Now by implicit derivation of the x-coordinate, by the equality τ = τ (x, t) = u(x, t), and by the Implicit function theorem (see [10,Theorem 9.28]) one obtains and it is easy to see that u, implicitly defined this way, is a local solution of (1).
If, however, u ∈ Υ is not smooth for some t > 0 and x ∈ R, then it has discontinuities called shocks which are positioned along piecewise smooth curves x = s(t) in the (x, t)-plane, called shock paths. By 2.2(iv), there are finitely many shock paths, which are all locally smooth. These paths may have singular points, they may cross or collide, but we can exclude these singularities without loss of generality.
We now continue our treatment by showing that the well-known Rankine-Hugoniot condition holds in our setting.

Lemma 2.4. A function u ∈ Υ is a solution of (3) if and only if the following
Rankine-Hugoniot condition is satisfied at all the shocks: by u + and u − we denote the one-sided limits of the solution u(x, t) from the left and from the right side of the shock, respectively, i. e. u + (x 0 , t) = lim x x 0 u(x, t)).
Proof. For the sake of simplicity we shall omit the arguments of functions whenever they are clear from the context. Take any a < b. If u ∈ Υ is smooth for x ∈ (a, b) and t > 0, it solves (1) on (a, b) and we have Now assume that u has a shock on (a, b) × {t} with shock path x = s(t). From the condition (iv) of Definition 2.2 it follows that x = s(t) is a locally C 1 -function. Hence, we can compute Thus, u is a solution of (3) on (a, b) if and only if the Rankine-Hugoniot condition (10) is fullfiled at the shock. If u has more then one shock on (a, b) × {t}, we divide the interval according to the shocks and proceed in the same manner. Note that Definition 2.2(iv) implies that u has only finitely many shocks.
We can prove uniqueness of the solutions to (2)-(3) in the class Υ, provided that the flux is a concave function and the solutions only have jumps upwards. The same is true for convex flux and downwards jumps.
Proposition 2.5. Let f ∈ C 2 (R) be strictly concave. Let u, v ∈ Υ be two solutions of (3) with the properties (a) u − < u + and v − < v + at every shock for u and v, respectively, and Then the L 1 -norm Proof. We slightly modify the idea Lax used for proving [5,Theorem 3.4]. Let u, v ∈ Υ be two solutions to (3) Since they both lie on the same characteristic surface their difference u − v can change the sign only in the shocks of either (or both) of these two functions. Denote these sign-changing curves in (x, t)-plane by y k (t) and order them as We work on maximal open intervals of t where the curves y k (t) do not intersect and the number of these curves does not change. For almost every t ≥ 0 we can then write Now choose any interval (y k (t), y k+1 (t)), 1 ≤ k ≤ n. Without loss of generality we may assume that u(x, t) > v(x, t) on this interval, hence the absolute value under the integral defining p k (t) can be omitted.
Note that p k (t) is a (continuous) piecewise differentiable function of t. If either u or v has some shocks inside the interval, in each of them the Rankine-Hugoniot condition is satisfied by Lemma 2.4. Dividing the interval according to these shocks and computing the derivative according to this division, similarly as it was done in (11), we see that the values around the shocks cancel out and only the values in y k (t) and y k+1 (t) are important. Therefore we shall from now on assume that none of u and v has shocks inside the interval (y k (t), y k+1 (t)).
Thus we may assume that function p k (t) is differentiable and we will now compute its derivative. As usual we will omit the arguments of functions whenever possible.
Since u and v solve (1) inside the interval and since the shock paths are piecewise differentiable we have From now on we shall explain the calculations only for the left endpoint of the interval, the right endpoint can be treated symmetrically. By assumption, u > v inside the interval, hence u − v changes sign at the endpoints and all the jumps are upwards. Taking all this into account we see that u has a shock in y k while v may have a shock or not (in the latter case we take v − = v + ). Furthermore we have Applying the Rankine-Hugoniot condition for the speed of shock y k for u we see that (14) equals since by concavity of f and condition (16) the first factor is smaller or equal to 0 while u + − v + > 0.
Following the same line of arguments for the right endpoint y k+1 we obtain the same conclusion for (15).
Note that the condition (a) of Proposition 2.5 is actually an entropy condition (see e. g. [6, p. 36]).

The equal-area method
We now describe the equal-area method for obtaining the solutions starting from the bounding characteristic surfaceΓ defined in (9). First note that the transformation G t defined in (7) preserves area, since its Jacobian equals 1. Hence the area bounded by the curves γ t and the x-axis remains unchanged in time and equals the initial area given by ∞ −∞ h(x) dx (compare the shaded areas in Figure  2). IntersectingΓ with the plane t = t 0 for some fixed time t 0 yields γ t 0 defined in (8), see also Figure 2. Our strategy is to insert vertical cuts to γ t 0 in such a way that the areas of the cut-off lobes coincide. Thus the initial area will be preserved.
Carrying out this procedure for all t we obtain a bounded, piecewise continuous function u(x, t), whose graph, without the added vertical surfaces, is contained in Γ, see  Lemma 3.1. Let f ∈ C 2 (R) and let h be a piecewise C 1 -function with compact support. Then for the solution u obtained by the above described equal-area method the following holds.
(a) At every shock u satisfies the Rankine-Hugoniot condition (10) and has only jumps upwards: u − < u + . (b) The shock paths are piecewise C 1 -curves.
For ξ 1 < ξ 2 we close the "S-curve" γ t ([ξ 1 , ξ 2 ]) with line segment between endpoints γ t (ξ 1 ) and γ t (ξ 2 ), see Figure 4. The signed area defined by this closed curve by Green's Theorem equals We define the mapping F : In points where the determinant of the 2 × 2 Jacobian matrix ∂F ∂ξ is equal to and is nonzero in jumps given by our method.
If in addition to (18), we have the equal-area condition p t 0 (ξ 0 1 , ξ 0 2 ) = 0, then the Implicit function theorem (see [10,Theorem 9.28]) gives the existence of two C 1functions ξ 1 (t), ξ 2 (t), such that for all t in some neighborhood of t 0 the following holds: This means that the equal-area condition holds for all parameters t in this neighborhood. Moreover, the Implicit function theorem gives us the formula Using symbolic computation (Mathematica [13]) one obtains . By above, the shock path is locally a C 1 -function, which proves (b). Differentiating this function at the point t = t 0 , we finally get the Rankine-Hugoniot condition (10) In the case when f is concave, f is a decreasing function and the curves γ t are for t > 0 inclined to the left (regarding x-axis). It then follows from the construction that the solution u only has jumps upwards, hence also the assertion (a) is proved.
We have proved the local smoothness of the shock paths, a nice property which is often assumed in advance and rarely verified. Assuming a finite number of shocks we are finally able to prove that the equal-area method yields the unique solution of our problem.
Theorem 3.2. Let f ∈ C 2 (R) be a strictly concave function and h a piecewise C 1 -function with compact support such that the function has only finitely many local extrema for every t ≥ 0. Then the above described equal-area method yields the unique solution u ∈ Υ to (3).
Proof. First observe that conditions (i), (ii), (iii), and (v) of Definition 2.2 are trivially fulfilled for a function u obtained by the equal-area method. Since the number of local extrema of the function (20) is finite, so is the number of shocks obtained by the equal-area method. Hence the condition (iv) of Definition 2.2 is satisfied by Lemma 3.1. Moreover, using Lemmas 3.1, 2.3, 2.4, and 2.5 we see that u ∈ Υ is the unique solution of (3).
A brief comment is in order here. For our theoretical approach we need to assume that there are only finitely many shocks. We are not aware of any explicit conditions in terms of functions f and h to meet this assumption (some generic conditions are given in [4,11]). In praxis however, it is not difficult to check the finiteness of the number of local extrema of the function x t (ξ) in (20) for any given f , h, and t (note that our procedure gives a solution for any fixed time t!). Moreover, numerically this condition is always satisfied, since we use polygonal approximation for continuous curves.

The algorithm
We shall now describe the algorithm based on the procedure introduced in the previous section. The solution at some fixed time t 0 is obtained directly, without need to proceed in time.
We start by taking a polygonal approximation K 0 for the continuous curve (21) γ t 0 := G t 0 (γ 0 ), see (7)- (8). Traveling along the curve we gradually 'equalize' the areas. The obtained graph of solution is a subset of K 0 (see Figure 6).

Iterative
Step of the Algorithm 4.1. Let the parametrization (x(τ ), y(τ )), τ ∈ R, of the curve K i on the i-th step be such that x(τ ) is increasing on the far ends of the interval (−∞, ∞). First we define three significant points for the curve K i (see Figure 5): (1) β = x(τ 1 ) is the first local maximum of x(τ ) in the direction of the increasing parameter τ . If such a maximum does not exist, we are done and K i is the graph of the weak solution (with redundant vertical lines in the jumps). (2) α = x(τ 2 ) is is the first local minimum of x(τ ) from τ 1 onwards. Since the function x(τ ) in not bounded from above, such a minimum always exists. (3) γ is the minimum of β and the first next local maximum of x(τ ). If such a maximum does not exist, let γ = β. p 1 p 2 Figure 5. Significant points α, β, γ, and δ on the curve K i and the areas p 1 (δ) and p 2 (γ).
Now we compute the areas. For any x 0 ∈ (α, β) denote by p 1 (x 0 ) the area bounded by the line x = x 0 and the part of the curve K i that contains (x(τ 1 ), y(τ 1 )). For any x 0 ∈ (α, γ) denote by p 2 (x 0 ) the area bounded by the line x = x 0 and the part of the curve K i that contains (x(τ 2 ), y(τ 2 )). For x ∈ (α, γ) let Then p 1 (x) is continuously decreasing while p 2 (x) and p(x) are continuously increasing functions and p(α) = 0 − p 1 (α) < 0. We distinguish two cases: (1) If p(γ) ≥ 0, let δ ∈ (α, γ] be the only zero of the function p(x), therefore p 1 (δ) = p 2 (δ). The curve K i+1 is obtained from K i where the parts of K i that determine p 1 (δ) and p 2 (δ) are replaced by the vertical line. (2) If p(γ) < 0, then p 1 (γ) > p 2 (γ) and by continuity and monotonicity of the function p 1 (x) there exists only one δ ∈ (γ, β) which satisfies p 1 (δ) = p 2 (γ) (note that p 1 (β) = 0). The point δ together with the areas p 1 (δ) and p 2 (γ) is marked on the Figure 5. The new curve K i+1 is obtained from K i by replacing those parts of K i that determine p 1 (δ) and p 2 (γ) by a vertical line.
In Figure 6 there is an example of the resulting steps of the above algorithm applied to the function given in (25) in time t 0 = 4.25. We shall briefly describe the method we use to compute the areas needed on each step of Algorithm 4.1. Let D be a polygon, determined by the points T 1 (x 1 , y 1 ), . . . , T n (x n , y n ) (we orient them in counterclockwise direction) and let T 0 (x 0 , y 0 ) be any point in the plane. Then the signed area of the triangle T 0 T i T i+1 can be computed by (23) p 0,i,i+1 := 1 2 By Green's Theorem one can easily see that the area of the polygon D then equals where T n+1 = T 1 .

Numerical results
We have programmed our equal-area method in Mathematica [13] and first compared the results with the basic Godunov method (which we have also implemented in Mathematica). Further we have compared our method to an advanced Godunov method, which is a basis of the widely-used software package Clawpack [2]. Finally, we have made a comparison to the very recent software package Particleclaw [9], which uses a Langrangean particle method and some information on the characteristics. Both software packages are freely available on the web.
The results of these tests are very good. Our algorithm performs favorably both in terms of time efficiency and accuracy. The solutions agreed and ours turned out to be even more accurate. Figure 7 contains graphs of the solution obtained by our method and both above mentioned software packages for the initial condition   Time complexity of our algorithm depends mostly on finding zeros of a function, obtained by the computation of areas of polygons. We used the secant method and typically 7 to 12 iterations (9 on average) were needed for 10 −14 accuracy. The number of necessary steps of the Algorithm 4.1 is bounded above by the number of stationary points of the function x t 0 (ξ) defined in (20).
We can approximate the error of the position of the shock. Assuming that the original curve γ t 0 lies in an ε-neighborhood of polygonal approximation line, we have an approximation of the area between the "S-curves" as in Figure 5: where l is the length of the "S-curve", ∆x is the displacement of the true shock and s is the height of the shock. This gives an approximation of the displacement ∆x: ∆x . = ε l s .
In Figure 6 polygonal approximation with 1000 points has ε less than 6.10 −4 and ∆x is approximated with 3.10 −3 . The method is quadratical, i.e. doubling the number of points would result in decreasing the value of ε by factor 4.
We see the following advantages of our equal-area method.
• The method is by its nature exactly conservative.
• The solution is computed for any given fixed time. Hence, the errors do not accumulate in time. • The method is accurate -the quality of the approximation relies only on the quality of the starting approximation of the curve γ t 0 with a polygonal line K 0 . • There is no need for separate treatment of the rarefaction waves, they are created on the way where appropriate. • The obtained shocks are sharp and propagate with correct speed. Their position is obtained automatically by equalizing the appropriate areas. Moreover, the shock paths are obtained easily by computing the solution for some selected times and then simply projecting the shocks from the surface to the (x, t)-plane (see Figure 8).