Approximate solutions to the initial value problem for some compressible flows in presence of shocks and void regions

For the natural initial conditions $L^1$ in the density field (more generally a positive bounded Radon measure) and $L^\infty$ in the velocity field we obtain global approximate solutions to the Cauchy problem for the 3-D systems of isothermal and isentropic gases, the 2-D shallow water equations and the 3-D system of collisionnal self-gravitating gases. We obtain a sequence of functions which are differentiable in time and continuous in space and tend to satisfy the equations in the sense of distributions in the space variables and in the strong sense in the time variable. The method of construction relies on the study of a specific family of two ODEs in a classical Banach space (one for the continuity equation and one for the Euler equation). Standard convergent numerical methods for the solution of these ODEs can be used to provide concrete approximate solutions. It has been checked in numerous cases in which the solutions of systems of fluid dynamics are known that our constuction always gives back the known solutions. It is also proved it gives the classical analytic solutions in the domain of application of the Cauchy-Kovalevska theorem.

Continuying a previous study of self gravitating pressureless fluids [11] we construct sequences of approximate solutions for 3-D compressible isothermal gases and isentropic gases in presence of shocks and void regions; this applies also to the shallow water equations and to the system of collisionnal selfgravitating isothermal gases. The need of the search of mathematical solutions to the initial value problem for compressible fluids in several dimensions when shocks show up has been recently pointed out by P. D. Lax and D. Serre [22,31].
The starting point of our attempt [7,8,9] was the search of a numerical scheme permitting mathematical proofs of partial results aiming to connect the functions produced by the scheme and the system one tries to solve: indeed, in absence of well defined solutions, it is natural to try to prove that the "approximate solutions" from the scheme tend to satisfy the equations in the sense of distributions, which makes sense in absence of properly defined mathematical solution, see the appendix. Sequences of functions that tend to satisfy the equations could be a provisional substitute of solutions, since one has to cope with an absence of known usual weak solutions in 2-D and 3-D in presence of shocks [22] p. 143, [31] p. 143 and [26] p. 150. These sequences permit in particular to explain results observed from computing [22,23]. In order to obtain full proofs, we replaced the original scheme [7,8,9] by ODEs in Banach spaces whose solutions provide the sequence of "approximate solutions", at the same time as one retains the possibility of computing the solutions of these ODEs to check that the method gives the correct results.
From a physical viewpoint the equations of fluid dynamics are mared with some imprecision since they do not take into account some minor effects and the molecular structure of matter. It is natural to expect these equations and their imprecision should be stated in the sense of distributions in the space variables. Sequences [(x, y, z, t) −→ U (x, y, z, t, )] of approximate solutions in the sense of distributions enter into this imprecision for > 0 small enough. Therefore these sequences appear to provide a natural representation of physical solutions. Of course one has checked that the sequences of approximate solutions we construct always give the correct solution in all the numerous tests that were performed. Uniqueness of the limit when → 0 for a general class of sequences of approximate solutions containing those we construct has been proved in the case of linear systems with analytic coefficients [13]. For the equations of fluid dynamics considered here, in absence of a known similar uniqueness result in a suitable family of approximations, one has to content with numerical confirmations.
These sequences of approximate solutions have been introduced under the name of weak asymptotic solutions by Danilov, Omel'yanov and Shelkovich [15], as an extension of Maslov asymptotic analysis, and they have proved to be an efficient mathematical tool to study creation and superposition of singular solutions to various nonlinear PDEs [15,16,17,18,27,32,33,34]. A weak asymptotic solution for the system is a sequence (U 1, , . . . , U p, ) of functions such that ∀i ∈ {1, . . . , p}, ∀ψ ∈ C ∞ c (R 2 ) ∀t when → 0 if we consider the strong derivative in t and weak derivatives in x, y as this will be the case in this paper. In short the sequence (U 1, , . . . , U p, ) tends to satisfy the equations. Of course the U i, s are also chosen so as to satisfy the initial condition stated in a natural sense.
The approximate solutions that we construct satisfy uniform L 1 bounds in density (consequence of mass conservation) but, in presence of pressure, one does not obtain a bounded velocity; this should not be too much unexpected: indeed various instances are known in which the ideal inviscid equations give infinite velocity [20] sections 6.6.2 , 7.4.3, [14] section 14.9.1. The concept of sequences of approximate solutions in this paper permits infinite velocity at the limit → 0.
We construct a weak asymptotic method for the 3-D systems of isothermal and isentropic gases that we state in 2-D for convenience (immediate extension to 3-D) in the form The notation is: ρ = ρ(x, y, t) is the gas density, (u, v) = (u(x, y, t), v(x, y, t)) is the velocity vector in the x, y directions respectively, p = p(x, y, t) is the pressure and K is a constant. The problem is studied on the n-dimensional torus T n = R n /(2πZ) n , n = 1, 2, 3. We combine the system of isothermal gases and self-gravitation to obtain the classical system of collisionnal selfgravitating gases which models Jeans gravitational instability [4,6]: a large enough cloud of gas possibly at rest collapses gravitationally besides pressure untill an equilibrium is reached with pressure, forming a star or a planet. A mathematical proof in this paper asserts that the numerical result observed is a depiction of an approximate solution of the equations with arbitrary precision, only limited by the cost of calculations.
Then, since the equations are similar, we construct weak asymptotic solutions for the 2-D shallow water equations stated in the form where h = h(x, y, t) is the water elevation, (u, v) = (u(x, y, t), v(x, y, t)) is the velocity vector in the x, y directions respectively, a = a(x, y) is the bottom elevation assumed to be of class C 2 and g = 9.8. The problem is studied on the n-dimensional torus T n , n = 1, 2.
The method of proof is the statement and study of a system of ODEs (one for the continuity equation, one for the Euler equation) in the Banach space C(T n ), n = 1, 2, 3, of continuous functions. We expose the proof in 1-D since the extension to several dimensions is straightforward. One uses the theory of ODEs in the Lipschitz case and a priori estimates needed to prove existence of a global in positive time solution. It has been inspired by the method introduced in [11] in the pressureless case and in [10] to put in evidence continuations of the analytic solutions after the time of the analytic blow up.
Numerical calculations of these approximate solutions can be done easily by solving the above ODEs in Banach space by classical convergent numerical methods for ODEs such as the explicit Euler order one method or the Runge Kutta RK4 method. This has given the known solutions in all tests, which could have been expected since the method in this paper is issued from the numerical scheme in [7,8,9], for which a large amount of numerical verifications for initial conditions in which the solution is known has been done on classical and demanding tests of different nature (Sod [35], Woodward-Colella [39], Toro [36,37], Lax [22,23], Bouchut-Jin-Li [2], LeVeque [24], Cherkov-Kurganov-Rykov [5], Evje-Flatten [19]).
2. Approximate solutions of the system of isothermal gases and numerical confirmations.
We state the classical system of isothermal gases in 1-D in the form: Therefore the last term in the Euler equation (11) is in the familiar form ρ ∂ ∂x Φ = ∂ ∂x p, with p = Kρ. Here Φ is the density of the body force caused by the pressure: ρ ∇Φ = ∇p. Formulation (12) excludes void regions. This apparent defect will be repaired by the term ρ ∇Φ in the Euler equation: it will be proved theoretically and checked numerically that the construction works even in presence of void regions. We use the notation We study the system in 1-D since the extensions to 2-D and 3-D are straightforward, following these extensions in section 6 of [11]. The 2-D extensions of formulas (15,16) below are given as (44, 45) below. In 1-D we state the method as the system of ODEs (15,16) complemented by the formulas (17,18): for some β > 0 to be made precise later, for which we will prove that ρ(x, t, ) > 0, thus permitting division. The term β in (15) is needed because of the specific form of the state law of isothermal gases and is not needed for the isentropic gases and shallow water equations.
We approximate the state law (12) in the form in which the term N has been introduced to permit void regions. We introduce an auxiliary function , so that the family {φ µ } µ tends in the sense of distributions to the Dirac measure when µ → 0 + . We use a real number α, 0 < α < 1, to be made more precise later. The convolution in (18) is justified by the fact (12) is a state law and, as such, is physically valid only in space regions larger than those in which the basic conservation laws (10,11) are valid.
We first establish a priori inequalities to prove existence of a global solution to (15)(16)(17)(18). For fixed > 0 we assume existence of a solution continuously differentiable on [0, δ( )[ (with a right hand-side derivative at t = 0) having the following properties for each fixed : Note that m and M depend on .
The aim of the following a priori inequalities is to obtain bounds independent on m and M to replace (20,21) for fixed > 0, depending on the initial condition and on , in order to prove that the solution can be extended for t > δ( ). All constants, denoted const in the proposition below, depend on the initial condition and possibly on a bounded time interval, but not on the values m and M in the a priori assumption (20,21) and not on . The independence on will be basic when they will be used later at the limit → 0 to prove that the solutions of the ODEs tend to satisfy the equations.
The dependence of const in t is explicitely stated when this appears clearer.
We assume ρ 0 and u 0 are given initial conditions with the properties ρ 0 ∈ L 1 (T) and u 0 ∈ L ∞ (T) and that ρ 0 and u 0 are regularizations of ρ 0 and u 0 ∈ C(T), with uniform L 1 and L ∞ bounds respectively (independent on ), and ρ 0 (x) > 0 ∀x ∈ T.
Inequality (23) is proved as follows: (18) implies Using (22) when ρ(x − y, t, ) ≥ 1 (then log(ρ + N ) < ρ), and, using that the above log is bounded in absolute value by const.log( 1 N ) when ρ(x − y, t, ) < 1, one obtains which gives the desired result (const does not depend on t from (22) when t ranges in a bounded interval). Now let us prove inequality (24). From (15) and the assumption that the solution of the ODE is of class C 1 on [0, δ( )[, valued in the Banach space C(T), one obtains for fixed > 0 and for dt > 0 small enough (t + dt < δ( )) that (20,21). Therefore since ρu ± ≥ 0 one can invert (27): where the new o has still the property that o(., t, Applying the analog of (27) with ρu in place of ρ, with the supplementary term ρ ∂ ∂x Φ from (16) and absence of the term β , one obtains (from 20, 21): where the new o has the same property as in (27) for fixed . For dt > 0 small enough the first term in the second member is a barycentric combination of u(x − , t, ), u(x, t, ) and u(x + , t, ) (which are in numerator in factor of ρ inside ρu), dropping the term β dt. For fixed the quotient ρ(x,t+dt, ) ρ(x,t, ) tends to 1 when dt → 0 (use (20, 21 and 27). Finally it follows from (23) and (28) that with uniform bound of o when t ranges in a compact set in [0, δ( )[. One obtains the bound (24) by dividing the interval [0, t] into n intervals [ it n , (i+1)t n ], 0 ≤ i ≤ n − 1, and applying (29) in each subinterval: summing on i and using that o( t n ) → 0 when n → ∞, see more details in [11]. Now let us prove inequalities (26). From (15), since ρ, u + and u − are positive (13,20). Therefore, from (24,25), Let v(x, t, ) = ρ(x, 0, )exp(− k( ) t). Then, using assumption (20) to divide by ρ, By integration, since ρ and v have same initial condition and are positive, which is the left hand-side inequality (26). Now let us prove the right handside inequality (26).
From the positiveness of the two terms ρu ± in (15) From (24,25) Since this holds for all x Gronwall's inequality implies For fixed > 0, if 0 < λ < 1 and Ω λ : (15)(16)(17)(18) with variables X = ρ, Y = ρu have the Lipschitz property on Ω λ with values in C(T), with Lipschitz constants uniform in Ω λ , see [11] section 4 for details. The existence of a unique global solution to (15)(16)(17)(18) for fixed is obtained from the a priori estimates in proposition 1 from classical arguments of the theory of ODEs in Banach spaces in the Lipschitz case as exposed in section 4 of [11].
To check (38) one has to prove from (18) that when → 0. To this end we share the integral (39) into the two parts (40, 41) below and we prove that each of them tends to 0 when → 0. Let and Now Since ρ(x, t, ) ≥ 0 from (26), using (22) in the case ρ(x, t, ) > 1 and using N in the case ρ(x, t, ) ≤ 1, as in the proof of (23), one has |I| ≤ const.log( 1 ) α ; therefore I → 0 when → 0. Now (41) and the mean value theorem give if min(ρ) denotes the inf of ρ(x, t, ) for fixed t and when x ranges in T.
The problem is to obtain a suitable inf. bound of min(ρ). The term β in (15) has been introduced for this purpose. Indeed, from (15), i.e. setting A := const 1 1+3α and B := β , the fact from (24) that u ∞ ≤ const 3α on any bounded time interval implies that on such interval For fixed t > 0 and for > 0 small enough in which we state the constant as const(t) since this constant depends on t through (24). We have From (42) |J| ≤ const(t). N −1−β−3α and it suffices to choose β + 3α < N − 1 to have that J → 0 when → 0.
It follows from the proof that the possible presence of void regions, approximated in the initial conditions by ρ 0 (x) ≥ ∀x, does not cause a problem although they appear a priori excluded by the formulation (12). Note also that the presence of in a denominator in (26) allows concentrations of matter, that can also be accepted in initial conditions by choosing ρ 0 (x) ≤ 1 ∀x: the initial condition ρ 0 can even be a bounded Radon measure. One could notice that formula (24) allows the possibility of infinite velocity at the limit = 0. This should not be troublesome: the remark that the ideal equations (2-6) could lead to solutions with infinite velocity in certain circumstances has been known since long time [14]  Finally we have proved that under the initial conditions ρ 0 ∈ L 1 (T), u 0 ∈ L ∞ (T), ρ 0 ≥ 0, and approximating the initial conditions by a family ( The result extends easily to 2-D and 3-D following section 6 in [11]. We state the 2-D equations (2,3): for some β > 0 to be made precise later, Numerical confirmations. The discretization in time is standard from numerical schemes for ODEs: for fixed > 0 the Lipschitz properties of the ODEs (15)(16) permit to prove convergence of the explicit Euler order one method. The discretization in space is as follows: the space R n is discretized by cells of length in each direction. The physical variables are constant in the cells. Numerical tests from the explicit Euler order one method and the RK4 Runge Kutta method for the solution of the ODEs (15,16) have shown that the weak asymptotic method always gives the correct solutions, even in presence of void regions, as this is the case in figures 1, 2 and 3.
In figure 1 we present a demanding test from [2] where an explicit solution is given. The numerical solution coincides with the explicit solution. It is the Riemann problem ρ g = 0, ρ d = 1, u g = 0 = u d at time T = 0.5 and with K = 0.04, which has been considered in [2] p. 154 and p. 157: this test is difficult due to the void region on the left in the initial conditions. The test has been done with dt = 0.00002, = 0.001, 2000 space steps.

Convergence to the analytic solution.
In this section we present a slight modification of the weak asymptotic method in section 2, which has better properties, in particular one can prove that when the initial data are analytic the weak asymptotic method gives the classical analytic solutions.
First one can check that the proof of the weak asymptotic method in section 2 holds without change: this is based on the fact that the presence of any fixed value δ > 0 does not affect significantly the bounds and that one uses only the formula u = u + − u − once one has replaced v(u) by u + + u − in (48) and in the ODE satisfied by ρu. Second, numerical tests (figure 2) show that the presence of a fixed δ > 0 affects the numerical results in the same way as a viscosity.
If we state D(u) = v(u)−|u| 2 and u ± = |u|±u 2 , i.e. u ± are the values u ± used in section 2, we obtain from (47, 48) The first term in the second member gives the formula used in section 2. The last term is some kind of vanishing viscosity.
Numerical confirmation. In figure 2 the discretization is exactly the same as in section 2. The difference is the use of v(u) = (u 2 + δ 2 ) 1 2 instead of |u|. We observe that δ has the same influence as a viscosity coefficient: for δ = 0.001 (top panel) and δ = 0.1 one obtains the exact solution; for δ = 1 one observes some viscosity effects which become quite important for δ = 2 (bottom panel). Now from the replacement of |u| by an analytic regularization v(u) we prove that the weak asymptotic method gives the classical analytic solution at the limit → 0, which was proved in [10] in a linear case when the function u has a fixed sign to avoid the singularity of the function absolute value at 0. The proof is given in the 1-D case since the multidimensional case is identical. The proof consists in applying the abstract nonlinear Cauchy-Kovalevska theorem of Nirenberg and Nishida [25] for each > 0 small enough, with results uniform in .
Recall of an abstract nonlinear Cauchy-Kovalevska theorem [25]. By definition a scale of Banach spaces is a family of Banach spaces (E s ) s , 0 < s ≤ s 0 , such that ∀s, s ∈]0, s 0 ], s > s ⇒ E s ⊂ E s with inclusion of norm ≤ 1. Let v 0 ∈ E s 0 be given. If R > 0 we denote by B s (v 0 , R) the open ball in the Banach space E s of center v 0 and radius R.
Let (E s ) 0<s≤s 0 be a scale of complex Banach spaces. Consider the Cauchy problem dv dt Assume the existence of R > 0, C > 0, K > 0 such that properties i) and ii) hold i) ∀s, s / 0 < s < s < s 0 the map u −→ G(u) is holomorphic from B s (v 0 , R) into E s , and satisfies a Lipschitz property in the sense The abstract Cauchy-Kovalevska theorem states: Then ∃ a number a > 0 and a unique holomorphic function t −→ v(t) which ∀s < s 0 maps {t ∈ C/|t| < a(s 0 − s)} into B s (v 0 , R) and is solution of (50).
The domain of v depends on its range through the number s: one understands that the functions v relative to various values of s stick together. The proof is a holomorphic form of the theorem in [25], setting there u = v − v 0 , F (u, t) = G(v). Property (51) implies the Lipschitz property stated in [25] if u, v ∈ B s (v 0 , R): The method of proof is the iteration method with adequate bounds, see [25]. The successive iterates lie in B s (v 0 , R). The number a > 0 depends only on R, K, C (formulas 13 p. 630 and end of p. 632 in [25]). Now we can prove the following coherence result.
The real number s 0 > 0 is chosen small enough so that the initial conditions v 0 := v 0 = (ρ 0 , ρ 0 u 0 ), which are independent on and are holomorphic extensions of the given real analytic initial data, are elements of the space E s 0 .
Then for any fixed > 0 small enough we consider the map G defined from (15,16) For fixed we apply the abstract Cauchy-Kovalevska theorem with the analytic initial condition v 0 = (ρ 0 , ρ 0 u 0 ) independent on and with the map G . We obtain a solution v = (ρ , (ρu) ) defined for |t| < a(s 0 − s) taking values in B s (v 0 , R) ∀s < s 0 (the domain of v depends of s: the same abuse of language has been done in the statement of the theorem above). The key of the proof lies in that the assumptions of the abstract Cauchy-Kovalevska theorem are satisfied uniformly in , therefore the properties (domains and bounds) of the solution will be also independent of .
The detailed formulas are more complicated due to the presence of u = ρu ρ in u ± inside the formula of G , and from the definition of u ± in (46) through the function v. To clarify we give the formulas in the particular case u ≥ 0 and v(u) = |u| = u which implies u + = u and u − = 0 and avoids longer formulas. Then from (58) To simplify the notation we set X = ρ, Y = ρu; then (60) Therefore If (X, Y ) ∈ B s (v 0 , R) ⊂ E s then the values X(x + iy), Y (x + iy) are defined for x ∈ T, |y| < s and are bounded in sup norm by v 0 s + R. Further X(x + iy) > b > 0; the terms N and β in (60, 61) are not used to check the validity of (54) for G uniformly in .
The abstract Cauchy-Kovalevska theorem asserts the existence of a solution v : t −→ (ρ (., t), (ρu) (., t)) ∈ B s (v 0 , R) to the ODEs (15, 16) defined on a domain of time {|t| < a(s 0 −s)} which is independent of since domains and bounds on the data are uniform in . Since the domains and bounds of the solutions are independent of the set of the functions v is a normal family of holomorphic functions on {|t| < a(s 0 − s)} × T×] − s, +s[ for some a > 0 independent of . From any sequence (v n ) n one can extract a convergent subsequence. This subsequence converges to the classical analytic solution of the system of PDEs (10,11). Therefore the whole family (v ) converges to the classical analytic solution.
4. Sequence of approximate solutions to the system of isentropic gas equations.
The system of isentropic gas equations (2,3,4,6) is stated in the form We state the system of ODEs as follows: we state the ODE (15) without the β term which is no more needed because the state law (64) is different from (12), i.e.
(in which u± are given by (13,14), but could also be given by (46, 47) with the replacement of |u| by v(u)). Then we state the ODE (16), the formula (17), and we replace the formula (18) by As usual the convolution in (66) is justified by the fact that the state law is obtained from measurements which are always done in a space region which is not too small. As in section 2 we assume initial conditions ρ 0 and u 0 defined on T. To obtain the a priori inequalities we assume (19)(20)(21). Then we obtain the same a priori inequalities as in proposition 1: the statements (22 without the β term), (23,24,25 with 2α in place of 3α) and (26) hold as in proposition 1.
proof. The proof of (22) in the present case is identical to the proof in proposition 1 without the β term. For the proof of (23) here we have Since 0 < γ − 1 ≤ 1 and since from (22) ρ ∈ L 1 (T) with ρ L 1 (T) independent of t and , then a fortiori ρ γ−1 ∈ L 1 (T) with bounds independent of t and . Therefore The proofs of (24) and (26) in the present context are identical to those in proposition 1.
From the a priori estimates one obtains existence and uniqueness of a global solution in time t ∈ [0, +∞[ and in space x ∈ T. The weak asymptotic method for the analogs of (36, 37) is proved as in section 2; here it suffices to have α < 1 4 since (23) is stated with 2α. For the state law (64) we have to prove that ∀ψ ∈ C ∞ c (R), ∀t
Finally we have obtained: let the initial conditions ρ 0 ∈ L 1 (T), more generally a positive bounded Radon measure, u 0 ∈ L ∞ (T) and ρ 0 ≥ 0. 5. Sequence of approximate solutions to the system of selfgravitating collisionnal gases.
Since the 1-D system is far simpler and the 3-D system quite analogous to the 2-D system, we state the system of isothermal collisionnal selfgravitating gases ([4] p. 49, [6] p. 207, [28] p. 460, [29] p. 231) in 2-D for convenience: where ρ, u, v are as usual the density and the components of the velocity vector, Φ press and Φ grav are respectively the density of body force per unit mass caused by the pressure and the gravitation. The system is similar in 3-D with (71) replaced by the usual 3-D elementary solution of the Poisson equation, const r , r = x 2 + y 2 , instead of const.logr in 2-D.
The weak asymptotic solutions in 1-D, 2-D, 3-D for the system (67-71) are adaptations of the ones in section 2 in absence of selfgravitation, in which Φ(x, y, t, ) there is replaced by the sum Φ press (x, y, t, ) + Φ grav (x, y, t, ) in which Φ press is given by (18) and in which The proofs (on the torus T n , n = 1, 2, 3, always considered in this paper for simplification) take place in the Banach space C(T n ), n = 1, 2, 3, and are not modified relatively to the isothermal case above, using the multidimensional pattern in section 6 of [11]. One obtains again a weak asymptotic solution.
Remark. For the modelling of large scale structure formation of the universe one uses periodic assumptions to model the cosmological principle that the universe is isotropic at large scales [6] p. 305. Therefore, for this problem, the study on the torus T 3 is more realistic than on the whole space R 3 .
6. Sequence of approximate solutions to the shallow water equations.
In this section we construct a weak asymptotic solution for the 2-D shallow water equations stated in the form where h = h(x, y, t) is the water elevation, (u, v) = (u(x, y, t), v(x, y, t)) is the velocity vector in the x, y directions respectively, a = a(x, y) is the bottom elevation assumed to be of class C 2 and g = 9.8. The initial condition and the function a are assumed to be periodic in the x, y variables (with period 2π for convenience). We give a proof in 1-D since the 2-D extension is straightforward following the pattern exposed in section 6 of [11].
We state h = ρ for convenience because of the similarity with the other systems considered before. Then the 1-D shallow water equations are: i.e.
. The proofs of the analogs of (24,26) are similar to those in proposition 1.
The existence of a global solution for fixed is obtained from these a priori estimates.
It remains to prove that the solution of the system of ODEs (80-83) provides a weak asymptotic solution for system (77-79) when → 0. For (77, 78) the proof in the isothermal case applies. For the state law (79) we have to check that where Φ is given by (82) and where f ( ) → 0 when → 0. To check (84) it suffices from (82) to consider the integral The last term is O( α ) from the L 1 loc property (22) of ρ since integration takes place on a compact set.
Finally we have proved: let the initial conditions ρ 0 ∈ L 1 (T), u 0 ∈ L ∞ (T), ρ 0 ≥ 0, a a function of class C 2 on T. Approximate the initial conditions by a family (ρ 0 , u 0 ) ∈ (C(T)) 2 , such that ρ 0 (x) > 0 ∀x, The method extends to 2-D following [11] section 6. One obtains coherence with the classical analytic solution as in proposition 2. Numerical tests from [36] have given results similar to those in section 2 and 3, with the possible use of v(u) in place of u. In figure 3 we have tested the three Riemann problem tests 1, 2 and 3 (from top to bottom) in [36] pp. 109-124. Top: we obtain the exact solution in test 1 with 500 space steps, dt=0.0001, an averaging on 3 cells to represent the convolution (82) on the state law with coefficients α, 1 − 2α, α with α = 0.1, no averaging needed on the initial condition and δ = 0 in the formula v(u) = (u 2 + δ 2 ) 1 2 . Middle: we obtain the exact solution in test 2 with 2000 space steps, dt=0.00004, an averaging for the state law (82) with α = 0.1, a similar averaging with α = 0.1 for the initial condition and δ = 1. Bottom: we obtain the exact solution in test 3 with 5000 space steps, dt = 8.10 −6 , averagings with α = 0.1 for the state law (82) and the initial condition, δ = 0.5. Due to a void region in initial condition (ρ = 1 on the left and 0 on the right) one has replaced the right value ρ = 0 by 10 −10 since the proof of weak asymptotic method requests ρ 0, (x) > 0 ∀x.

Conclusion.
In view of providing a substitute of solutions that could explain the observed numerical results such as [22,23], we have constructed approximate solutions (up to any given accuracy in the sense of distributions) to the general Cauchy problem on the n-dimensional torus T n , n = 1, 2, 3, for some standard equations of fluid dynamics in presence of shocks and void regions.
Up to our knowledge, sequences of global approximate solutions to the general initial value problem for the equations of fluid dynamics with a complete mathematical proof that they tend to satisfy the equations had not been constructed previously. The statement of the method is based on a family of two nonlinear ordinary differential equations in a classical Banach space of functions, for which a priori estimates permit to prove existence-uniqueness of global solutions (in space and in positive time). The continuity equation (mass conservation) gives a L 1 control on density which, from its use inside the state law giving pressure, permits from the Euler equations some control on velocity.
The method of proof allows numerical calculations from standard convergent numerical methods for ODEs such as the explicit Euler order 1 method and the RK4 method. They have given the known solutions in all tests.  [36] pp. 109-124; left: water elevation, right water velocity. The results were obtained with the explicit Euler order one method applied to the ODEs (80, 82). One observes coincidence with the exact solutions given in [36], although these tests have been chosen on purpose to be very demanding from the numerical viewpoint due to the presence of void regions in the water elevation.
Uniqueness of an "admissible" class of such sequences of approximate solutions, containing those constructed in this paper ("existence part"), such that all sequences in this class give same values for all time at the limit → 0 ("uniqueness part") has been obtained so far only for some linear equations [13,3]. Passage to the limit by some kind of weak compactness has not been considered so far.
We recall known examples which suggest the possible absence of distribution solutions also in some instances of fluid dynamics.
• The system ∂u ∂t produces shock waves involving a power δ n of the Dirac delta measure as a continuation of analytic solutions after their blow up. To prove this fact we produce a sequence of approximate solutions with the same properties as those in this paper and we compute them from a convergent numerical scheme for ODEs. One observes numerically the powers δ n [10]. System (85-86) is not strictly hyperbolic but this fact has only served to create a δ-wave in v to insert into (87) and such δ-wave can be created as well in strictly hyperbolic systems: see v in system (88-89) below.
• The Keyfitz-Kranzer system shows a more subtle phenomenon. The approximate solutions u , v explicitely calculated in formula (8) in [21] from the Dafermos-Di Perna viscosity, and obtained again numerically in [30] from the usual viscosity have a limit u, v in the sense of distributions when → 0: u shows a simple discontinuity and v shows a Dirac mass located on the discontinuity; they move with constant speed. When one inserts such u and v in u 2 − v the Dirac measure in v subsists untouched in u 2 − v and therefore the term ∂(u 2 −v) ∂x has a nonzero δ term. This term cannot be compensated in ∂u ∂t because u is a simple discontinuity. Equation (88) cannot be satisfied. Similarly in (89) ∂v ∂t shows a δ term that cannot be compensated in . At the level of the approximate sequence (u , v ) with = 0 the compensations take place because the function u has a term (the last term in the first formula (8) in [21]) which could be intuitively refered to has a "square root" of a Delta measure, i.e. an object whose square tends in the sense of distributions to the Dirac measure when → 0: then u 2 has a δ-term that compensates the δ term of v in u 2 − v; a similar verification holds for equation (89). This shows that when one passes to the limit → 0 in the approximate solutions in the sense of distributions one obtains a result (u, v) that, when plugged into the equations, does not always satisfy them.