Solitary wave formation in the compressible Euler equations

Let us consider the Cauchy problem for the 1D Euler equations (see (2) below). The background density is given by

	$\displaystyle\hat{\rho}(\chi)$	$\displaystyle=\begin{cases}1/4&0\leq\chi-\lfloor\chi\rfloor<1/2\\ 7/4&1/2\leq\chi-\lfloor\chi\rfloor<1.\end{cases}$		(1a)
and we define the material coordinate $x=\int\rho(\chi)d\chi$. Then the initial data is given by
	$\displaystyle p(x,0)$	$\displaystyle=1+\frac{3}{20}\exp(-x^{2}/25)$		(1b)
	$\displaystyle\rho(x,0)$	$\displaystyle=p(x,0)^{1/\gamma}\hat{\rho}(\chi(x))$		(1c)
	$\displaystyle u(x,0)$	$\displaystyle=0.$		(1d)

This initial data corresponds to a gas in which the initial density and temperature vary periodically in space, such that the pressure is constant over most of the domain. Near $x=0$, there is an initial positive pressure perturbation that will split into two opposite-traveling pulses.

Figure 1 shows the resulting pressure solution at four different times. The initial pulse first steepens, but at later times, rather than the usual N-wave formation typical of 1D hyperbolic conservation laws, we observe the emergence of a series of solitary waves, reminiscent of the behavior of dispersive nonlinear wave equations such as Korteweg-de Vries (KdV).

By way of contrast, in Figure 2 we show the behavior of the same disturbance when the background density is taken to be constant and equal to the average of the oscillating density from the previous example: $\hat{\rho}(x)=1/2$. In this case we see the N-wave formation, and the solution will asymptotically decay (as $t\to\infty$) to a spatially-uniform state. This is typical of solutions to 1D hyperbolic conservation laws. Note also that the solution propagates significantly faster.

We start from the one-dimensional Euler equations

	$\displaystyle\rho_{t}+(\rho u)_{\chi}$	$\displaystyle=0,$
	$\displaystyle(\rho u)_{t}+(\rho u^{2}+p)_{\chi}$	$\displaystyle=0,$		(2a)
	$\displaystyle\left(\frac{1}{2}\rho u^{2}+\rho e\right)_{t}+\left(\frac{1}{2}\rho u^{3}+\rho eu+up\right)_{\chi}$	$\displaystyle=0,$

representing conservation of mass, momentum, and energy. Here $\rho(\chi,t)$ is the mass per unit volume, $u(\chi,t)$ the gas velocity, $p(\chi,t)$ the gas pressure, and $e$ is the internal energy per unit mass. We consider a polytropic gas, for which

$$e=\frac{1}{\gamma-1}\frac{p}{\rho},$$

(3)

where $\gamma=c_{p}/c_{v}$ is the polytropic constant, given by the ratio of specific heats respectively at constant pressure and constant volume. The Euler equations can be conveniently written in Lagrangian form by adopting the mass coordinate:

$$x=x(\chi,t)=x_{0}+\int_{\chi_{0}}^{\chi}\rho(\tilde{\chi},t)\,d\tilde{\chi}.$$

(4)

The difference $x-x_{0}$ measures the mass of the gas (per unit area) between position $x_{0}$ and $x$ at time $t$. Notice that for any $t$ the relation between the mass Lagrangian coordinate $x$ and the Eulerian coordinate $\chi$ relation is invertible.

In the new coordinates, the Euler equations can be written as

	$\displaystyle v_{t}-u_{x}$	$\displaystyle=0,$
	$\displaystyle u_{t}+p_{x}$	$\displaystyle=0,$		(5)
	$\displaystyle\left(\frac{1}{2}u^{2}+e\right)_{t}+(up)_{x}$	$\displaystyle=0.$

where $v=1/\rho$ is the specific volume.

In this paper we are interested in studying the propagation of weakly nonlinear waves. In particular, we plan to show that such waves that propagate on an unperturbed state with constant pressure and zero velocity and periodic density will never form shocks. For such a reason we shall also consider a different form of the equations. Manipulating the energy equation, making use of the first principle of thermodynamics and of the first two equations, one obtains that, for smooth solutions, the entropy density $s=S/c_{v}$ does not depend on time (see, for example, [27], Chapter 6). The last equation of (5) can then be replaced by

$$s_{t}=0$$

(6)

The entropy for a polytropic gas is given by

$$s=\log(pv^{\gamma})+{\rm const},$$

which can be written as

$$s-s_{*}=\log\left(\frac{p}{p_{*}}\left(\frac{v}{v_{*}}\right)^{\gamma}\right),$$

where $p_{*}$ and $v_{*}$ denote a particular reference state, and $s_{*}$ the corresponding entropy density. Solving for $p$, the functional relation $p=p(\rho,s)$ becomes

$$p=p_{*}e^{s-s_{*}}(v_{*}/v)^{\gamma}.$$

(7)

From Eq. (6) it follows that the entropy density is a function of space only,

$$s=s(x)$$

In particular, $s(x)$ is determined by considering an initial condition which is a perturbation of a stationary state

$$u_{0}=0,\>p_{0}=p_{*},\>v_{0}=v_{0}(x),$$

so that

$$e^{s-s_{*}}=(v_{0}(x)/v_{*})^{\gamma},\quad s-s_{*}=\gamma\log(v_{0}(x)/v_{*}).$$

In the classical statistical theory of gases, the entropy is defined up to an additive constant. Without loss of generality we take $s_{*}=0$.

Using Eq. (7), and considering that $s=s(x)$, the $3\times 3$ system (5) reduces to the $2\times 2$ system with space-dependent coefficients

	$\displaystyle v_{t}-u_{x}$	$\displaystyle=0$		(8a)
	$\displaystyle u_{t}+p(v,s(x))_{x}$	$\displaystyle=0.$		(8b)

We consider an unperturbed situation in which the initial density $v_{0}(x)$ is a periodic function of $x$, with period $\delta$²²2In the examples in this work we use initial data with average density equal to 1, which conveniently means that the period is the same in both the Eulerian and Lagrangian coordinates. and assume that the initial condition is a perturbation of amplitude $O(\delta)$, with a typical wavelength large compared to $\delta$, see Figure 3.

Given that the specific volume has oscillations whose amplitude is $O(1)$, we prefer to use $p$ in place of $v$ as dependent variable, since the pressure is expected to be much less oscillatory than $v$.

We therefore reformulate the problem as

	$\displaystyle\left(\frac{\partial v}{\partial p}\right)_{x={\rm const}}p_{t}-u_{x}$	$\displaystyle=0$		(9)
	$\displaystyle u_{t}+p_{x}$	$\displaystyle=0$		(10)

Since the entropy depends only on space and not on time, we have

$$\left(\frac{\partial v}{\partial p}\right)^{-1}_{x={\rm const}}=\left(\frac{\partial p}{\partial v}\right)_{x={\rm const}}=\left(\frac{\partial p}{\partial v}\right)_{s={\rm const}}=-\frac{\gamma p}{v}=-c^{2}$$

where $c$ denotes the sound speed in Lagrangian coordinates. It represents the mass (per unit surface) crossed by a small pressure perturbation per unit time. It is related to the classical Eulerian sound speed $c_{E}$ by the relation

$$c_{E}=vc$$

In the field variables $(p,u)$ the system can be written as

	$\displaystyle p_{t}+c^{2}u_{x}$	$\displaystyle=0$		(11a)
	$\displaystyle u_{t}+p_{x}$	$\displaystyle=0$		(11b)

where

$$c^{2}=\frac{\gamma p}{v}=c_{*}^{2}\left(\frac{p}{p_{*}}\right)^{1+1/\gamma}K(x)=c_{*}^{2}\left(\frac{p}{p_{*}}\right)^{1+1/\gamma}v_{*}/v_{0}(x)$$

and we introduced the quantity

$$K(x)\equiv e^{-s(x)/\gamma}=v_{*}/v_{0}(x).$$

In this work we consider the Cauchy problem for propagation of waves in a compressible gas in which the entropy field is initially periodic in space. We consider the problem in Lagrangian coordinates and suppose that no shocks form, so that $s=s(x)$ and (11) holds. We then derive a constant-coefficient (i.e. effective medium) equation that approximates the p-system (11), using multiple-scale perturbation analysis. An analysis of the system (11) with a different constitutive relation and physical interpretation was carried out by LeVeque & Yong, who found that the resulting equations admit traveling solitary wave solutions, in good agreement with numerical experiments [10]. In that work the authors considered a special exponential equation of state and focused on piecewise constant materials. Here we revisit the analysis in the context of compressible gas dynamics and consider general periodic variation of $s(x)$. We also make a connection with related analysis in [7] and propose a criterion for shock stability in a compressible gas with a spatially periodic entropy.

In all numerical experiments in this work, we take $V_{*}=p_{*}=1$.

In what follows we make use of certain averaging operators, which appear in the literature but are recalled here for the convenience of the reader.

First, the integral (or average) of $f$ over one period, denoted by $\braket{f}\in\mathbb{R}$, is defined as

$$\braket{f}:=\int_{0}^{1}f(y)\,dy.$$

Second, the fluctuating part of the function $f$, denoted by $\{f\}$, is defined as

$$\{f\}(y):=f(y)-\braket{f}.$$

Finally, the fluctuating part of the antiderivative of the fluctuating part, denoted by $\llbracket f\rrbracket$, is defined for any $y$ as

$$\llbracket f\rrbracket(y):=\left\{\int_{0}^{y}\left\{{f(x)}\right\}\,dx\right\},$$

that is,

$$\llbracket f\rrbracket(y)=\int_{0}^{y}\left\{f\right\}(x)dx-\int_{0}^{1}\int_{0}^{\tau}\left\{f\right\}(x)dxd\tau.$$

(12)

Clearly, we have

$$\braket{\{f\}}=0\quad\quad\text{and}\quad\quad\braket{\llbracket f\rrbracket}=0.$$

Some useful properties of the $\llbracket\cdot\rrbracket$ operator, which will be used in this work, are provided in [8, Appendix A] .

LeVeque & Yong [10] studied the propagation of elastic waves in a periodically-varying medium, using the model

	$\displaystyle\sigma_{t}-K(x)G(\sigma)u_{x}$	$\displaystyle=0$		(13a)
	$\displaystyle\rho(x)u_{t}-\sigma_{x}$	$\displaystyle=0.$		(13b)

The p-system (11) is a special case of (13), with $\rho(x)=1$ and the notational changes

	$\displaystyle\rho(x)$	$\displaystyle\leftrightarrow 1$		(14a)
	$\displaystyle u(x,t)$	$\displaystyle\leftrightarrow u(x,t)$		(14b)
	$\displaystyle\sigma(x,t)$	$\displaystyle\leftrightarrow-p(x,t)$		(14c)
	$\displaystyle K(x)$	$\displaystyle\leftrightarrow e^{-s(x)/\gamma}=:K(x)$		(14d)
	$\displaystyle G(\sigma)$	$\displaystyle\leftrightarrow c_{*}^{2}\left(\frac{p}{p_{*}}\right)^{1+1/\gamma}=:G(p).$		(14e)

Here we have overloaded the definitions of $K$ and $G$ in order to facilitate comparison between the analysis here and in [10]. The correspondence (14) allows us to directly transfer results of the analysis conducted in [10] to the present setting. To do so, we introduce a fast spatial variable $y=x/\delta$, formally independent of $x$, such that partial derivatives are transformed as

$$\frac{\partial}{\partial x}\to\frac{\partial}{\partial x}+\delta^{-1}\frac{\partial}{\partial y}.$$

We suppose there exist power series for $p,u$ in terms of $\delta$:

	$\displaystyle p(x,y,t)$	$\displaystyle=p^{0}(x,t)+\delta p^{1}(x,y,t)+\delta^{2}p^{2}(x,y,t)+\cdots$		(15)
	$\displaystyle u(x,y,t)$	$\displaystyle=u^{0}(x,t)+\delta u^{1}(x,y,t)+\delta^{2}u^{2}(x,y,t)+\cdots.$		(16)

We can also expand $G(p)$ as a power series:

	$\displaystyle G(p)$	$\displaystyle=G(p^{0})+G^{\prime}(p^{0})(p-p^{0})+\frac{1}{2}G^{\prime\prime}(p^{0})(p-p^{0})^{2}+\cdots$		(17)
		$\displaystyle=\frac{c_{*}^{2}}{p_{*}^{1+1/\gamma}}\left((p^{0})^{1+1/\gamma}+(1+1/\gamma)(p^{0})^{1/\gamma}(p-p^{0})+\frac{1+1/\gamma}{2\gamma}(p^{0})^{1/\gamma-1}(p-p^{0})^{2}+\cdots\right).$		(18)

Following the analysis in [10] we obtain

	$\displaystyle\overline{p}_{t}+\frac{G}{\braket{K^{-1}}}\overline{u}_{x}$	$\displaystyle+\delta^{2}\frac{\mu\braket{K^{-1}}}{G}\left(\frac{4\overline{p}_{t}\overline{p}_{tt}G^{\prime}}{G}+\overline{p}_{t}^{3}\left(\frac{G^{\prime\prime}}{G}-\frac{3G^{\prime 2}}{G^{2}}\right)-\overline{p}_{ttt}\right)$
	$\displaystyle+\delta^{4}\frac{\zeta}{\braket{K^{-1}}}$	$\displaystyle\left(\alpha_{1}\overline{p}_{tttt}\overline{p}_{t}-\alpha_{2}\overline{p}_{tt}\overline{p}_{ttt}+\alpha_{3}\overline{p}_{t}\overline{p}_{tt}^{2}+\alpha_{4}\overline{p}_{t}^{2}\overline{p}_{ttt}+\alpha_{5}\overline{p}_{t}^{3}\overline{p}_{tt}+\alpha_{6}\overline{p}_{t}^{5}+\alpha_{7}\overline{p}_{ttttt}\right)$	$\displaystyle={\mathcal{O}}(\delta^{5}),$		(19a)
	$\displaystyle\overline{u}_{t}+\overline{p}_{x}$	$\displaystyle=0,$			(19b)

where $G=G(\overline{p})$ and

	$\displaystyle\mu$	$\displaystyle=\frac{\braket{\llbracket K^{-1}\rrbracket^{2}}}{\braket{K^{-1}}^{2}}$		(20a)
	$\displaystyle\zeta$	$\displaystyle=\braket{K^{-1}(\llbracket K^{-1}\rrbracket)^{2}}$		(20b)
	$\displaystyle\alpha_{1}$	$\displaystyle=-9\frac{G^{\prime}}{G}$		(20c)
	$\displaystyle\alpha_{2}$	$\displaystyle=-15\frac{G^{\prime}}{G^{3}}$		(20d)
	$\displaystyle\alpha_{3}$	$\displaystyle=\frac{62G^{\prime 2}-\frac{37}{2}GG^{\prime\prime}}{G^{4}}$		(20e)
	$\displaystyle\alpha_{4}$	$\displaystyle=\frac{46G^{\prime 2}-13GG^{\prime\prime}}{G^{4}}$		(20f)
	$\displaystyle\alpha_{5}$	$\displaystyle=\frac{-11G^{2}G^{\prime\prime\prime}-160G^{\prime 3}+108GG^{\prime}G^{\prime\prime}}{G^{5}}$		(20g)
	$\displaystyle\alpha_{6}$	$\displaystyle=\frac{-G^{3}G^{(4)}+9G^{2}G^{\prime\prime 2}+75G^{\prime 4}+13G^{2}G^{\prime\prime\prime}G^{\prime}-\frac{165}{2}GG^{\prime 2}G^{\prime\prime}}{G^{6}}$		(20h)
	$\displaystyle\alpha_{7}$	$\displaystyle=\frac{1}{G^{2}}.$		(20i)

Note that we always have $\mu\geq 0$ and $\zeta\geq 0$. Here we have assumed that $K(y)$ is translation-even³³3We say a function defined on a periodic domain is translation-even if it can be made even by a shift. and applied [8, Proposition 5] in order to simplify this expression somewhat; for a more general function $K(y)$ (19) includes several additional terms of order $\delta^{4}$. In the examples below we consider piecewise-constant or sinusoidal functions $K(y)$, both of which are shift-periodic. Additionally, several terms that are present in [10] instead vanish here because of the correspondence $\rho(x)\leftrightarrow 1$.⁴⁴4Note also that there is a typo in [10, eqn. (5.18)]: the expressions for $C_{12}$ and $C_{22}$ are reversed.

System (19) contains high order derivatives in time. This form of the equations is not convenient for two reasons. First, the original system of PDE’s is a first order hyperbolic system of two equations, therefore the initial value problem is well posed when two initial profiles are assigned. In contrast, system (19) requires four initial conditions, namely the initial values of $\overline{u}$, $\overline{p}$, $\overline{p}_{t}$, and $\overline{p}_{tt}$. Furthermore, system (19) is linearly unstable (see Appendix A.1 ). For these reasons, we shall convert some time derivatives to space derivatives, keeping the same order of consistency ${\mathcal{O}}(\delta^{5})$ in the asymptotic expansion.

It is possible to convert all of the $t$-derivatives in the higher-order terms in (19) to $x$-derivatives, in order to write it as a first order system of evolution equations. However, as noted in [8], this approach leads to a system that is linearly unstable for large wavenumber modes. As in [8], we instead keep one $t$-derivative in the higher-order linear terms, in order to obtain a system that is linearly stable for all wavenumbers. This yields

	$\displaystyle\overline{p}_{t}+\frac{G(\overline{p})}{\braket{K^{-1}}}\overline{u}_{x}-\delta^{2}\mu\left(\overline{p}_{xxt}+\frac{G^{\prime}(\overline{p})}{\braket{K^{-1}}}\overline{p}_{xx}\overline{u}_{x}\right)+\delta^{4}\left(\frac{\zeta}{\braket{K^{-1}}^{3}}-\mu^{2}\right)\overline{p}_{xxxxt}$	$\displaystyle=\delta^{4}{\mathcal{N}}(\overline{p},\overline{u})+{\mathcal{O}}(\delta^{5})$		(21a)
	$\displaystyle\overline{u}_{t}+\overline{p}_{x}$	$\displaystyle=0.$		(21b)

Here ${\mathcal{N}}(\overline{p},\overline{u})$ is a function of $\overline{p}$ and $\overline{u}$ and their derivatives, in which every term is nonlinear:

	$\displaystyle{\mathcal{N}}(\overline{p},\overline{u})=$	$\displaystyle\frac{1}{\braket{K^{-1}}}\left(\frac{\zeta}{\braket{K^{-1}}^{3}}-\mu^{2}\right)\Bigg{[}\beta\left(2\frac{G^{\prime}}{G}\overline{p}_{xx}-\frac{G^{\prime\prime}}{G}(\overline{p}_{x})^{2}\right)-6G^{\prime}\overline{p}_{xxx}\overline{u}_{xx}-G^{\prime}\overline{p}_{xx}\overline{u}_{xxx}-6\frac{(G^{\prime})^{2}}{G}(\overline{p}_{x})^{2}\overline{u}_{xxx}$
		$\displaystyle-\left(8\frac{(G^{\prime})^{2}}{G}+6G^{\prime\prime}\right)\overline{p}_{x}\overline{p}_{xx}\overline{u}_{xx}-6G^{\prime\prime}\overline{p}_{x}\overline{p}_{xxx}\overline{u}_{x}-6\frac{G^{\prime}G^{\prime\prime}}{G}(\overline{p}_{x})^{3}\overline{u}_{xx}$
		$\displaystyle+\frac{\beta}{\braket{K^{-1}}}\left(\frac{G^{\prime\prime}}{G}-2(G^{\prime})^{2}\right)(\overline{u}_{x})^{2}+\frac{1}{\braket{K^{-1}}}\left(2(G^{\prime})^{2}-6GG^{\prime\prime}\right)\overline{u}_{x}(\overline{u}_{xx})^{2}$
		$\displaystyle+\frac{1}{\braket{K^{-1}}}\left(2(G^{\prime})^{2}-GG^{\prime\prime}\right)(\overline{u}_{x})^{2}\overline{u}_{xxx}+\frac{1}{\braket{K^{-1}}}\left(2\frac{(G^{\prime})^{3}}{G}-G^{\prime}G^{\prime\prime}\right)\overline{p}_{xx}(\overline{u}_{x})^{3}$
		$\displaystyle+\left(-9\frac{G^{\prime}G^{\prime\prime}}{G}-3G^{(3)}\right)(\overline{p}_{x})^{2}\overline{p}_{xx}\overline{u}_{x}-2\frac{G^{\prime}G^{(3)}}{G}(\overline{p}_{x})^{4}\overline{u}_{x}$
		$\displaystyle+\frac{1}{\braket{K^{-1}}}\left(4\frac{(G^{\prime})^{3}}{G}-4G^{\prime}G^{\prime\prime}-6GG^{(3)}\right)\overline{p}_{x}(\overline{u}_{x})^{2}\overline{u}_{xx}$
		$\displaystyle+\frac{1}{\braket{K^{-1}}}\left(2\frac{(G^{\prime})^{2}G^{\prime\prime}}{G}-2(G^{\prime\prime})^{2}-2G^{\prime}G^{(3)}-GG^{(4)}\right)(\overline{p}_{x})^{2}(\overline{u}_{x})^{3}+G^{\prime}\overline{p}_{xxxx}\overline{u}_{x}$
		$\displaystyle-2G^{\prime}\overline{p}_{x}\overline{u}_{xxxx}-2\frac{(G^{\prime})^{2}}{G}(\overline{p}_{xx})^{2}\overline{u}_{x}\Bigg{]}$

where for shortness we set

$\displaystyle+\delta^{4}\frac{\zeta G^{\prime\prime}}{2\braket{K^{-1}}^{4}}(\overline{p}_{xx})^{2}\overline{u}_{x}\beta$

$\displaystyle=G\overline{u}_{xxx}+G^{\prime}\overline{p}_{xx}\overline{u}_{x}+2G^{\prime}\overline{p}_{x}\overline{u}_{xx}+G^{\prime\prime}\overline{u}_{x}(\overline{p}_{x})^{2}.$

The number of terms appearing here in the evolution equation for $p$ is much larger than in [10], because in that work the equation of state was chosen such that $G^{\prime\prime}=0$. On the other hand, the evolution equation for $u$ here is much simpler, due to the fact that the variable coefficient $\rho(x)$ present in [10] is replaced by a constant here.

In this section we study the well-posedness, for small initial data, of the initial value problem for the model equations derived in the previous section. As mentioned above, the linearization of system (19) is unstable for all wave numbers, and the initial value problem is ill-posed (see Appendix A.1).

Next we consider the system (21). We linearize around an equilibrium configuration $(\overline{u},\overline{p})=(0,p_{*})$. Neglecting terms of $O(\delta^{5})$, the linearized equations read

	$\displaystyle\overline{u}_{t}+\overline{p}_{x}$	$\displaystyle=0$		(22a)
	$\displaystyle\overline{p}_{t}+c^{2}\overline{u}_{x}-\delta^{2}\mu\overline{p}_{xxt}+\delta^{4}\nu\overline{p}_{xxxxt}$	$\displaystyle=0$		(22b)

where

$$c^{2}\equiv\frac{G(p_{*})}{\braket{K^{-1}}},\quad\nu\equiv\frac{\zeta}{\braket{K^{-1}}^{3}}-\mu^{2}$$

(23)

We look for solutions of the system (22) of the form

$$\overline{u}=\hat{u}e^{i(kx-\omega t)},\quad\overline{p}=e^{i(kx-\omega t)}$$

Plugging this ansatz into system (22) we obtain a homogeneous system for $\hat{p}$ and $\hat{u}$:

	$\displaystyle-i\omega\hat{u}+ik\hat{p}$	$\displaystyle=0$
	$\displaystyle-i\omega\hat{p}+ic^{2}k\hat{u}-i\delta^{2}k^{2}\omega\mu\hat{p}-i\delta^{4}\nu k^{4}\omega\hat{p}$	$\displaystyle=0.$

Non-trivial solutions of the above system exist provided the determinant of the coefficient matrix of the system vanishes. This gives the following dispersion relation

$$c^{2}k^{2}-\omega^{2}(1+\mu\delta^{2}k^{2}+\nu\delta^{4}k^{4})=0,$$

(24)

which can be solved for $\omega$:

$$\omega=\pm\frac{ck}{\sqrt{1+\mu\delta^{2}k^{2}+\nu\delta^{4}k^{4}}}=\pm ck\left(1-\frac{1}{2}\mu\delta^{2}k^{2}+\frac{1}{8}(3\mu^{2}-4\nu)\delta^{4}k^{4}+{\mathcal{O}}(\delta^{6}k^{6})\right).$$

For all the profiles $K(x)$ we have tested, $\nu>0$, and we conjecture that $\nu>0$ in general, which means that the systems admits linearly dispersive waves for all wave numbers $k\in\mathbb{R}$.

We now compare solutions obtained by solving the Euler equations, the p-system, and the homogenized equations. The Euler equations in Eulerian form and the variable-coefficient p-system are solved using Clawpack [13], specifically with the high-order SharpClaw algorithm based on fifth-order WENO reconstruction and a ten-stage, fourth-order strong stability preserving Runge-Kutta method [6, 5, 4]. The homogenized equations (21), including terms up to ${\mathcal{O}}(\delta^{4})$ are solved using a Fourier pseudospectral collocation method in space and the three-stage, third-order strong stability preserving Runge-Kutta method of Shu & Osher [22]. The code used to produce the results reported here is available online⁵⁵5ADD link to Github repository..

In the previous section we showed numerical evidence that under certain circumstances, nonlinear waves in the Euler equations do not break into shocks, and no entropy is produced. This suggests that the solution remains smooth for all times, presumably maintaining the regularity of the initial data. If that is the case, accurate numerical solutions do not require the use of a shock capturing scheme, and discretizations can be based on a non-conservative form of the equations. In the case of smooth initial data it is possible to adopt, for example, a pseudospectral discretization, similar to the one adopted for the numerical solution of the homogenized equations.

Here we write the Euler system as

	$\displaystyle\rho_{t}$	$\displaystyle=-(\rho u)_{\chi}$		(26a)
	$\displaystyle u_{t}$	$\displaystyle=-uu_{\chi}-p_{\chi}/\rho,$		(26b)
	$\displaystyle p_{t}$	$\displaystyle=-up_{\chi}-\gamma pu_{\chi}.$		(26c)

All spatial derivatives are computed by Fourier pseudospectral differentiation; for example, $u_{x}\approx{\rm IFFT}(ik{\rm FFT}(u))$. The system is integrated in time using the standard four-stage, fourth-order Runge-Kutta method. In the computation we use a Courant number 0.9. As usual, we integrate the equations in the domain $[-L,L]$ with periodic boundary conditions, and plot the solution only in the interval $[0,L]$. Note that if the initial density and pressure are even functions, and the initial velocity is an odd function, then the symmetry of the solution is maintained for all times.

As an example we compare the solution obtained by the high-order finite volume method described in the previous section, with the one obtained by the pseudospectral method. We take as initial condition

	$\displaystyle\rho(\chi,0)$	$\displaystyle=1+\cos(2\pi\chi)$		(27a)
	$\displaystyle u(\chi,0)$	$\displaystyle=0$		(27b)
	$\displaystyle p(\chi,0)$	$\displaystyle=\frac{3}{20}\exp(-\chi^{2}/16).$		(27c)

The result of the comparison is reported in Fig. 6. In order to get a fully resolved calculation we used 144,000 cells with Clawpack in the interval $[-600,600]$ and only 32K points with the pseudospectral code in the domain [-512,512]. It is not our intent here to make a detailed performance comparison, and the solvers are implemented in different languages, but we remark that the pseudospectral simulation takes about one tenth of the time of the FV simulation (about 10 minutes versus about 100 minutes), with both running on a M1 Max Macbook Pro with 32 GB of RAM.

We now look for traveling waves for system (21), keeping the term proportional to $\overline{p}_{xxxxt}$, and neglecting only the right hand side of (21a). Let us denote by $\nu$ the coefficient of $\overline{p}_{xxxxt}$, as in (23). Using as before the approximation $G^{\prime}(p)\approx G^{\prime}(p_{*})$ in system (21), and proceeding as above, we obtain the following fourth-order equation for $u$:

$$F[u]:=\frac{\delta^{2}\nu}{\mu}u^{{}^{\prime\prime\prime\prime}}-\left(u^{\prime\prime}-\frac{G^{\prime}(p_{*})}{2V\braket{K^{-1}}}(u^{\prime})^{2}+\frac{{\mathcal{G}}(p(u))-{\mathcal{G}}(p_{*})}{\delta^{2}\mu V^{3}\braket{K^{-1}}}-\frac{u}{\delta^{2}\mu}\right)=0.$$

(32)