# On hyperbolic approximations for a class of dispersive and diffusive-dispersive equations

Rahul Barthwal, Firas Dhaouadi, Christian Rohde

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rahul Barthwal, Firas Dhaouadi, Christian Rohde. On hyperbolic approximations for a class of dispersive and diffusive-dispersive equations. 2025. hal-05405026

HAL Id: hal-05405026

https://hal.science/hal-05405026v1

Preprint submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On hyperbolic approximations for a class of dispersive and diffusive-dispersive equations

Rahul Barthwal*¹, Firas Dhaouadi†², and Christian Rohde‡¹

¹Institute of Applied Analysis and Numerical Simulation, University of Stuttgart, Pfaffenwaldring 57, D-70569 Stuttgart, Germany

²Institut de Mathématiques de Bordeaux, Université de Bordeaux, CNRS UMR 5251, Bordeaux INP, INRIA, F-33400, Talence, France

## Abstract

We introduce novel approximate systems for dispersive and diffusive-dispersive equations with nonlinear fluxes. For purely dispersive equations, we construct a first-order, strictly hyperbolic approximation. Local well-posedness of smooth solutions is achieved by constructing a unique symmetrizer that applies to arbitrary smooth fluxes. Under stronger conditions on the fluxes, we provide a strictly convex entropy for the hyperbolic system that corresponds to the energy of the underlying dispersive equation. To approximate diffusive-dispersive equations, we rely on a viscoelastic damped system that is compatible with the found entropy for the hyperbolic approximation of the dispersive evolution. For the resulting hyperbolic-parabolic approximation, we provide a global well-posedness result. Using the relative entropy framework [1], we prove that the solutions of the approximate systems converge to solutions of the original equations.

The structure of the new approximate systems allows to apply standard numerical simulation methods from the field of hyperbolic balance laws. We confirm the convergence of our approximations even beyond the validity range of our theoretical findings on set of test cases covering different target equations. We show the applicability of the approach for strong nonlinear effects leading to oscillating or shock-layer-forming behavior.

**Keywords.** Hyperbolic balance laws, dispersive and diffusive-dispersive equations, relative entropy, generalized KdV equations, non-classical shock waves.

**Mathematics subject classification.** 35A35, 35A01, 35G20, 35L45, 65M08.

## 1 Introduction

We consider for some interval $I \subset \mathbb{R}$ scalar balance laws of the form

$$u_t + f(u)_x = \varepsilon u_{xx} + \gamma u_{xxx} \quad \text{in } I_T := I \times (0, T), \ T > 0, \qquad (1.1)$$

with associated initial data

$$u(\cdot, 0) = u_0. \qquad (1.2)$$

In (1.1), $u = u(x, t) \in \mathbb{R}$ is the unknown. $f: \mathbb{R} \to \mathbb{R}$ denotes a smooth flux function, $\gamma$ is a nonzero real parameter, $u_0: I \to \mathbb{R}$ indicates the initial function. We consider both cases, i.e., $\varepsilon = 0$ leading

<footnote>*rahul.barthwal@mathematik.uni-stuttgart.de</footnote>

<footnote>†firas.dhaouadi@bordeaux-inp.fr</footnote>

<footnote>‡christian.rohde@mathematik.uni-stuttgart.de</footnote>
to a purely dispersive evolution and the diffusive-dispersive equation with $\varepsilon > 0$.

Many mathematical models for nonlinear wave propagation in the low-frequency regimes lead to (1.1) with $\varepsilon = 0$. They have been analyzed extensively and we can only mention some examples: The model (1.1) represents the Korteweg-De Vries (KdV) equation for $f(u) = u^2/2$ [2], the modified Korteweg De Vries (mKdV) equation for $f(u) = u^3$ [3], the generalized KdV equation for $f(u) = u^p + u$, $p \in \mathbb{N}$ [4], and the Gardner equation for $f(u) = u^3 + \mu u^2$ with $\mu \in \mathbb{R}$ [5]. The model (1.1) with $\varepsilon > 0$ can be considered as a toy model to describe phase transition dynamics when viscous and capillary effects are taken into account [3]. We refer the interested readers to the review paper [6] for more details on dispersive and diffusive-dispersive equations and their physical significance; see also [2, 7].

Apart from its significance in physics, the model (1.1) comes with several mathematical challenges. In particular, detecting non-classical shock waves, such as nonlinear dispersive shocks or undercompressive shocks due to the coupling of the third-order term with the flux function, has gained attention, see, e.g., [8–10]. In recent years, a lot of crucial work has been done for the analysis of diffusive-dispersive scalar balance laws of the form (1.1). For the cubic-flux case $f(u) = u^3$ (the mKdV-Burgers model), the existence and structure of undercompressive traveling waves were rigorously characterized, with subsequent studies addressing stability; see, e.g., [3, 6, 11, 12] and references therein. Later, undercompressive dynamics in certain diffusive-dispersive systems were observed numerically, offering strong confirmation and additional understanding of wave selection mechanisms [13, 14]. On the dispersive side, the broader theory of dispersive shock waves offers a unifying modulation framework for oscillatory shock-like structures generated by equations in the KdV/mKdV family and their viscous perturbations; see [6, 15].

Notably, the equation (1.1) comes with an energy-like structure. Sufficiently smooth solutions of (1.1) satisfy

$$ \frac{d}{dt} E[u] + \frac{d}{dx} Q[u] = -\varepsilon (\operatorname{sgn}(\gamma) f'(u) u_x^2 + |\gamma| u_{xx}^2), \quad (1.3) $$

with the energy $E[u]$ and the energy flux $Q[u]$ given by

$$ \begin{aligned} E[u] &= \operatorname{sgn}(\gamma) F(u) + \frac{|\gamma| u_x^2}{2}, \quad F'(u) = f(u), \\ Q[u] &= \frac{f(u)^2}{2} + |\gamma|(u_x f(u)_x - f(u)u_{xx}) - \varepsilon \operatorname{sgn}(\gamma)u_x f(u) - \varepsilon|\gamma|u_x u_{xx} \\ & \quad + \frac{\operatorname{sgn}(\gamma)\gamma^2 u_{xx}^2}{2} - \operatorname{sgn}(\gamma)\gamma^2 u_x u_{xxx}. \end{aligned} \quad (1.4) $$

If $\operatorname{sgn}(\gamma)F(u)$ is strictly convex or in other words $\operatorname{sgn}(\gamma)f'(u) \ge \delta$ for some $\delta > 0$, then the energy $E$ becomes a strictly convex Lyapunov functional and (1.3) becomes an energy inequality, which then provides some $H^1$-control on $u$. In particular, for the dispersive case ($\varepsilon = 0$), the energy $E$ is conserved. Relying on this energy structure and higher-order energy estimates, classical results for (generalized) KdV/mKdV establish local and global well-posedness in $H^s$-spaces, which provide the analytical foundation for studying solutions emanating from sufficiently regular data; see, e.g., [16–20] and references cited therein.

Beyond direct analysis of (1.1), there have been attempts to provide lower-order approximations for (1.1). Such relaxations are typically tailored to specific examples of (1.1) so as to preserve key structural features (e.g., energy preservation or dissipation) while simplifying the dynamics. For the diffusive-dispersive law (1.1) with $\gamma > 0$, Corli&Rohde [21] proposed an elliptic relaxation in which an auxiliary variable is introduced and coupled to a reduced-order evolution through a screened Poisson equation. Their approximation takes the form.

$$ u_t^\alpha + f(u^\alpha)_x = \varepsilon u_{xx}^\alpha - \alpha (u^\alpha - c^\alpha)_x, \quad -\gamma c_{xx}^\alpha = \alpha (u^\alpha - c^\alpha). \quad (1.5) $$

Relying on the energy structure of the approximate system, it was proven that as $\alpha \to \infty$, the sequence of solutions $\{(u^\alpha, c^\alpha)\}_{\alpha>0}$ of the initial value problem for (1.5) converge to some function $(u, c)$, such
that $u$ solves (1.1), (1.2). Extending this idea, Corli et al. [22] introduced a parabolic relaxation by adding temporal dynamics for the auxiliary variable. The parabolic relaxation system for $\gamma > 0$ takes the form.

$$u_t^\alpha + f(u^\alpha)_x = \varepsilon u_{xx}^\alpha - \alpha (u^\alpha - c^\alpha)_x, \quad \beta c_t^\alpha - \gamma c_{xx}^\alpha = \alpha (u^\alpha - c^\alpha), \quad (1.6)$$

where the parameter $\beta$ is related with the relaxation parameter $\alpha$ such that $\beta(\alpha) \to 0$ as $\alpha \to \infty$. By higher-order energy estimates, convergence of the sequence of solutions $\{(u^\alpha, c^\alpha)\}_{\alpha>0}$ of (1.6) to the solutions of (1.1), (1.2) has been established as $\alpha \to \infty$. While both elliptic and parabolic approximations capture the diffusive-dispersive behavior effectively, they still result in second-order and mixed-type systems. These cannot be treated by standard numerical methods. In particular, the numerical simulation becomes demanding in the stiff limit regime for $\alpha \to \infty$. On the other hand, first-order hyperbolic or second-order hyperbolic-parabolic systems admit a broad class of high-resolution methods (like, e.g., finite volume methods with classical Riemann solvers and limiters), making them straightforward to implement and comparatively inexpensive; see, e.g., [23, 24].

In this article, therefore, we propose alternative hyperbolic and hyperbolic-parabolic systems that approximate solutions of (1.1), (1.2) while retaining essential structural features and enabling efficient computation. The proposed approximations also rely on the introduction of an auxiliary variable, like in (1.5), (1.6), but let it evolve by special wave equations. However, our proposed system does not rely on the sign of dispersion in contrast to the elliptic and parabolic relaxation systems and considers both the cases $\gamma > 0$ and $\gamma < 0$.

In Section 2, we first formulate the suggested approximation models for the purely dispersive case $\varepsilon = 0$ in (1.1). By transformation of the unknowns, we can recast the approximation in the form of a first-order $(4 \times 4)$-system of balance laws; see (2.7). This system turns out to be strictly hyperbolic, admits explicit Riemann invariants, and preserves the Hamiltonian structure of (1.1). For monotone fluxes, we can provide a convex entropy which resembles the energy structure of (1.1) as in (1.3). In the second part of Section 2, we provide a hyperbolic-parabolic approximation in (2.26) using a damped wave operator. Besides local well-posedness, we establish a global well-posedness result for initial data close to constant equilibrium in Theorem 2.7. Again, we exploit the entropy structure which is inherited from the hyperbolic approximation (2.7) and ensures dissipative behavior. The main novel contributions of Section 2 are the two hyperbolic and hyperbolic-parabolic approximations (2.7) and (2.26).

Section 3 is devoted to prove the convergence of the solutions of the proposed systems to solutions of (1.1) using the relative entropy framework. As our main analytical contributions, we prove in Theorem 3.1 the convergence towards solutions of (1.1) in the purely dispersive case and verify in Theorem 3.3 the convergence in the diffusive-dispersive case. Let us note that we consider in Theorem 3.1 weak entropy solutions of the approximate system (2.7), whereas Theorem 3.3 concerns classical solutions of (2.26) which exist globally under the assumptions of Theorem 2.7.

In Section 4, we provide several test cases for a variety of dispersive and diffusive-dispersive equations. We show that the approximations (2.7) and (2.26) can be solved numerically by a standard approach from the realm of hyperbolic balance laws. The method applies to all test cases, including those with strong nonlinear effects. The study confirms the theoretical convergence results from Section 3. Conclusions and future outlooks are provided in Section 5.

We conclude the introduction with some remarks on the available literature on first-order approximations for diffusive-dispersive equations. Such models have long been studied in diverse settings. In particular, many hyperbolic relaxation models for higher-order partial differential equations have been suggested in recent years, see, e.g., [24–29] and references cited therein. Numerical investigations in all these works suggest a good convergence towards the solutions of the original equations, but mostly lack a rigorous convergence analysis or do not exploit or rely on the energy structure of the
original problem. In particular, a rigorous justification of the Favrie-Gavrilyuk approximation [26] to the Serre-Green-Naghdi model was provided in [30]; however, most of the other relaxation systems still lack rigorous convergence proofs. A more recent contribution by Giesselmann&Ranocha [31] addresses this gap by proving the convergence for several hyperbolic relaxation systems via the relative entropy framework. Building on this idea, we show that the same approach applies to our proposed approximate system (2.7) using its convex entropy. We refer the interested readers to [1, 32, 33] and references cited therein for more details on the relative entropy framework.

## Notations

Throughout the article, we use the following standard notations. We denote by $C^k(D)$, the space of $k$-times continuously differentiable functions defined in some domain $D$ and $k \in \mathbb{N} \cup \{\infty\}$. By $C_0^{\text{Lip}}(D)$, we denote the space of Lipschitz continuous functions having a compact support. Further, for $p \in [1, \infty) \cup \{\infty\}$, we denote the Lebesgue space on $D$ as $L^p(D)$ and the $L^p(D)$-norm as $\|\cdot\|_{L^p(D)}$, which for a function $g \in L^p(D)$ is defined as

$$ \|g\|_{L^p(D)} = \left( \int_D |g(x)|^p dx \right)^{\frac{1}{p}} \quad \text{for } p \in [1, \infty), \quad \|g\|_{L^\infty(D)} = \esssup_{x \in D} \{|g|\} \quad \text{for } p = \infty. $$

We denote by $H^k(D)$ the Sobolev space of order $k$ consisting of square-integrable functions whose weak derivatives up to order $k$ are also square-integrable, that is

$$ H^k(D) := \{g \in L^2(D) : \partial^m g \in L^2(D) \text{ for all multi-indices } m \text{ with } |m| \le k\}, $$

equipped with the norm

$$ \|g\|_{H^k(D)}^2 := \sum_{|m| \le k} \|\partial^m g\|_{L^2(D)}^2, $$

where $\partial^m g$ denotes the weak derivative of $g$ of order $m$.

Finally, for a Banach space $X$ and $T > 0$, we denote by $L^p(0, T; X)$, the Bochner space of measurable functions $g : (0, T) \to X$ with an induced norm

$$ \|g\|_{L^p(0,T;X)} = \left( \int_0^T \|u(t)\|_X^p dt \right)^{1/p}, $$

for $p \in [1, \infty)$, while for $p = \infty$

$$ \|g\|_{L^\infty(0,T;X)} = \sup_{t \in [0,T]} \|g(t)\|_X. $$

# 2 Hyperbolic and hyperbolic-parabolic approximations

In this section, we propose the relaxation approximations for both the purely dispersive equation (1.1) with $\varepsilon = 0$ in Section 2.1 and the diffusive-dispersive equation (1.1) with $\varepsilon > 0$ in Section 2.3. The approximations can be seen as extensions of the elliptic approach from [21], but replacing the elliptic operator by different types of second-order hyperbolic wave operators. These will be reformulated by introducing first- and second-order systems, which are shown to be of hyperbolic and hyperbolic-parabolic structure, respectively. In the dispersive case, we show that the proposed system is strictly hyperbolic and is locally well-posed via a unique symmetrizer. Moreover, we succeed in showing that the first-order system is derivable from Hamilton's principle and conserves an associated entropy. The approach is then extended to the dissipative case, for which we succeed in providing a complete entropy framework (see [1]) that is compatible with the energy related to the limit problem (1.1), (1.2). We exploit the entropy structure and prove local and global well-posedness results for the corresponding initial value problems.
## 2.1 A first-order approximation of the purely dispersive equation

We start from (1.1) with $\varepsilon = 0$. Thus, we consider for $\gamma \in \mathbb{R} \setminus \{0\}$, the dispersive equation

$$u_t + f(u)_x = \gamma u_{xxx}, \quad \text{in } I_T, \qquad (2.1)$$

subject to periodic boundary conditions and the initial condition

$$u(\cdot, 0) = u_0 \quad \text{in } I. \qquad (2.2)$$

The existence of classical solutions for such dispersive equations with periodic boundary conditions has been analyzed in various settings. In particular, for the generalized KdV equation, local and global well-posedness for high regularity solutions in $H^s$-spaces is a well-established result. We refer the interested readers to the classical work of Kenig et al. [16] in this regard, where the local well-posedness of $H^s$ solutions for $s > 3/4$ was developed. For global well-posedness of smooth solutions in $H^s$-spaces for modified and generalized KdV equations in both $\mathbb{R}$ and $\mathbb{T}$ (periodic BCs), we refer the interested readers to Colliander et al. [19], Farah et al. [20], and references therein. For a more general flux function and variable dispersion, we refer to the recent work of Mietka [34]. In particular, for results related to the periodic setting, we also refer to [35, 36].

As outlined in the introduction, an important property of classical solutions of (2.1) is the energy conservation such that for the energy/energy flux pair $(E, Q)$ from (1.4), we have the relation

$$\frac{d}{dt}E[u] + \frac{d}{dx}Q[u] = 0. \qquad (2.3)$$

We search now for a first-order approximation of (2.1) that resembles this gradient flow structure. Analogously to the screened Poisson equation (1.5) and the screened heat equation (1.6), we relax the Laplacian term $u_{xx}$ in (2.1) by a screened linear wave operator. For some relaxation parameter $\alpha > 0$, we seek then for a pair $(u^\alpha, c^\alpha)$ that satisfies in $I_T$ the approximate system

$$\begin{aligned} u_t^\alpha + f(u^\alpha)_x &= -\operatorname{sgn}(\gamma)\alpha (u^\alpha - c^\alpha)_x, \\ \beta c_{tt}^\alpha - |\gamma| c_{xx}^\alpha &= \alpha (u^\alpha - c^\alpha), \end{aligned} \qquad (2.4)$$

subject to periodic boundary conditions, and the prepared initial data

$$u^\alpha(\cdot, 0) = c^\alpha(\cdot, 0) = u_0, \quad c_t^\alpha(\cdot, 0) = -f(u_0)_x - |\gamma|u_{0,xxx}, \quad \alpha(u^\alpha - c^\alpha)(\cdot, 0) = -|\gamma|u_{0,xx} \quad \text{in } I. \qquad (2.5)$$

The parameter $\beta$ in (2.4)$_2$ will be expressed in the sequel in terms of $\alpha$. In the limit $\alpha \to \infty$ we expect to recover (2.1) if $\beta = \beta(\alpha)$ vanishes. Naturally, the characteristic speed associated with the wave operator (2.4)$_2$ given by $v_\beta := \sqrt{|\gamma|/\beta}$ goes to $\infty$ as $\beta \to 0$, in line with the instantaneous signal speed of the elliptic operator in (1.5).

Now, in order to reduce the coupled mixed-order system (2.4) to an equivalent purely first-order form, let us consider a change of variables following [24, 37], i.e.,

$$w^\alpha = c_t^\alpha, \quad p^\alpha = c_x^\alpha, \quad \psi^\alpha = \alpha (u^\alpha - c^\alpha). \qquad (2.6)$$

We then substitute these quantities into (2.4) and supply the expressions for the time-derivatives of $p^\alpha$ and $\psi^\alpha$ as closure equations to obtain a system of first-order balance laws of the form

$$\begin{aligned} u_t^\alpha + (f(u^\alpha) + \operatorname{sgn}(\gamma)\psi^\alpha)_x &= 0, \\ \psi_t^\alpha + \alpha (f(u^\alpha) + \operatorname{sgn}(\gamma)\psi^\alpha)_x &= -\alpha w^\alpha, \\ w_t^\alpha - v_\beta^2 p_x^\alpha &= \beta^{-1} \psi^\alpha, \\ p_t^\alpha - w_x^\alpha &= 0 \end{aligned} \qquad \text{in } I_T. \qquad (2.7)$$

Note that if $w^\alpha \to u_t$, $\psi^\alpha \to -|\gamma|u_{xx}$ and $p^\alpha \to u_x$ as $\alpha \to \infty$, we recover the original equation (2.1). Thus, the system (2.7) appears to be a consistent approximation for (2.1).
To maintain the consistency with the initial data from (2.5), the prepared initial conditions

$$u^{\alpha}(\cdot, 0) = u_0, \psi^{\alpha}(\cdot, 0) = -|\gamma|u_{0,xx}, w^{\alpha}(\cdot, 0) = -f(u_0)_x - |\gamma|u_{0,xxx}, p^{\alpha}(\cdot, 0) = u_{0,x} \text{ in } I \quad (2.8)$$

are supposed to hold.

In Theorem 3.1 we prove that solutions of the periodic initial value problem (2.7), (2.8) converge for $\alpha \to \infty$ to solutions of (2.1), (2.2) if $\beta = \alpha^{-1}$ is chosen.

## 2.2 Qualitative properties of the first-order approximation

### 2.2.1 Hyperbolicity, symmetrizability and local well-posedness

In this section we omit the superscript (.)$\alpha$ on the approximation unknowns and simply write $U = (u, \psi, w, p)^T$ for $U^\alpha = (u^\alpha, \psi^\alpha, w^\alpha, p^\alpha)^T \in \mathbb{R}^4$, if no confusion is caused.

In order to verify the hyperbolicity of the system (2.7), we consider a smooth solution $U$ and cast (2.7) into the balance law form

$$U_t + \mathbf{f}(U)_x = \mathbf{b}(U) \iff U_t + \mathbf{D}\mathbf{f}(U)U_x = \mathbf{b}(U), \quad (2.9)$$

where $\mathbf{D}\mathbf{f}(U)$ is the Jacobian of the flux $\mathbf{f} = \mathbf{f}(U)$. The flux and the linear source $\mathbf{b} = \mathbf{b}(U)$ are given by

$$\mathbf{f}(U) = \begin{pmatrix} f(u) + \operatorname{sgn}(\gamma)\psi \\ \alpha(f(u) + \operatorname{sgn}(\gamma)\psi) \\ -v_{\beta}^2 p \\ -w \end{pmatrix}, \quad \mathbf{b}(U) = \mathbf{S}U, \quad \mathbf{S} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & -\alpha & 0 \\ 0 & \beta^{-1} & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix},$$

such that the Jacobian of the flux is explicitly given by

$$\mathbf{D}\mathbf{f}(U) = \begin{pmatrix} f'(u) & \operatorname{sgn}(\gamma) & 0 & 0 \\ \alpha f'(u) & \alpha \operatorname{sgn}(\gamma) & 0 & 0 \\ 0 & 0 & 0 & -v_{\beta}^2 \\ 0 & 0 & -1 & 0 \end{pmatrix}.$$

The eigenvalues of $\mathbf{D}\mathbf{f}(U)$ compute as the real numbers $\lambda_1, \dots, \lambda_4$ given by

$$\lambda_1 = -v_{\beta}, \quad \lambda_2 = 0, \quad \lambda_3 = v_{\beta}, \quad \lambda_4 = f'(u) + \alpha \operatorname{sgn}(\gamma). \quad (2.10)$$

While the ordering $\lambda_1 < \lambda_2 < \lambda_3$ is direct, the placement of $\lambda_4$ within this ranking depends on the scaling of the free parameters $\alpha$ and $\beta$. In all what follows, we choose $\alpha$ such that $|\lambda_4| > \lambda_3$ and ensure $\lambda_4$ does not change sign under the sufficient condition

$$\alpha > \max |f'(u)|. \quad (2.11)$$

Under these conditions, depending on the dispersion sign one obtains

$$\lambda_1 < \lambda_2 < \lambda_3 < \lambda_4 \text{ for } \operatorname{sgn}(\gamma) = 1, \quad \lambda_4 < \lambda_1 < \lambda_2 < \lambda_3 \text{ for } \operatorname{sgn}(\gamma) = -1.$$

In both cases the approximate system (2.7) is strictly hyperbolic. The corresponding right eigenvectors $\mathbf{r}_1, \dots, \mathbf{r}_4 \in \mathbb{R}^4$ can be explicitly computed and are given by

$$\mathbf{R} := (\mathbf{r}_1 | \cdots | \mathbf{r}_4) = \begin{pmatrix} 0 & -\operatorname{sgn}(\gamma) & 0 & 1 \\ 0 & f'(u) & 0 & \alpha \\ v_{\beta} & 0 & -v_{\beta} & 0 \\ 1 & 0 & 1 & 0 \end{pmatrix}. \quad (2.12)$$

We proceed to prove that the hyperbolic system (2.7) can be symmetrized using a unique symmetrizer. To be precise, we prove the following.

**Lemma 2.1** (Symmetrizability and local well-posedness of the initial value problem for (2.7)). For $\alpha > \max |f'(u)|$, the hyperbolic system (2.7) can be converted into an equivalent symmetric hyperbolic form. In particular, the hyperbolic system (2.7) is a Friedrichs-symmetrizable system.
For

$$ \mathbf{U}_0 := (u_0, \psi_0, w_0, p_0)^\top \in (H^s(I))^4, s > 3/2, $$

there exists a time $T^* \in (0, \infty)$ such that the initial value problem for (2.7) and $\mathbf{U}(\cdot, 0) = \mathbf{U}_0$ in $I$ is solved in $I \times [0, T^*)$ by a unique function $\mathbf{U}$ such that

$$ \mathbf{U}(\cdot, t) \in C([0, T^*]; (H^s(I))^4) \cap C^1([0, T^*]; (H^{s-1}(I))^4) $$

holds for $s > 3/2$.

*Proof.* In order to prove this lemma, first consider the map

$$ \mathcal{T}: \mathbb{R}^4 \mapsto \mathbb{R}^4, \quad \mathbf{U} = \begin{pmatrix} u \\ \psi \\ w \\ p \end{pmatrix} \mapsto \mathbf{\Phi} = \begin{pmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \end{pmatrix}, \quad \text{with} \quad \begin{cases} \phi_1 = \operatorname{sgn}(\gamma)f(u) + \psi, \\ \phi_2 = \alpha u - \psi, \\ \phi_3 = v_\beta p, \\ \phi_4 = w. \end{cases} $$

The Jacobian of this transformation is given by

$$ D\mathcal{T}(\mathbf{U}) = \begin{pmatrix} \operatorname{sgn}(\gamma)f'(u) & 1 & 0 & 0 \\ \alpha & -1 & 0 & 0 \\ 0 & 0 & 0 & v_\beta \\ 0 & 0 & 1 & 0 \end{pmatrix}, $$

such that $\det(D\mathcal{T}(\mathbf{U})) = v_\beta\lambda_4 \neq 0$. This shows that the mapping $\mathcal{T}$ is bijective, and thus the system (2.9) can be converted into the equivalent form

$$ \mathbf{\Phi}_t + \tilde{\mathbf{A}}(\mathbf{\Phi})\mathbf{\Phi}_x = \tilde{\mathbf{b}}(\mathbf{\Phi}), $$

with the symmetric matrix $\tilde{\mathbf{A}} = \tilde{\mathbf{A}}(\mathbf{\Phi})$ given by

$$ \tilde{\mathbf{A}}(\mathbf{\Phi}) = \begin{pmatrix} \operatorname{sgn}(\gamma)\alpha + f'(u(\mathbf{\Phi})) & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -v_\beta \\ 0 & 0 & -v_\beta & 0 \end{pmatrix} $$

and $\tilde{\mathbf{b}}(\mathbf{\Phi}) = (0, -\alpha\phi_4, \psi(\phi_1, \phi_2), 0)^\top$ such that $\psi(\mathbf{\Phi})$ satisfies $\operatorname{sgn}(\gamma)f(u(\mathbf{\Phi})) + \psi = \phi_1$ and $\alpha u(\mathbf{\Phi}) - \psi = \phi_2$. Here, the quantity $u(\mathbf{\Phi})$ is simply the value of $u$, computed from $\mathbf{\Phi}$ through the inverse transform $\mathcal{T}^{-1}$.

Clearly, under this change of variables, the system (2.7) is a Friedrichs-symmetrizable system, and thus the local well-posedness of the Cauchy problem for classical solutions is a direct consequence of classical results on the theory of hyperbolic balance laws, see, e.g., [1]. $\square$

**Remark 2.1.** In Section 2.2.3 below we show that (2.7) is equipped with an entropy/entropy-flux pair with a convex entropy provided that $\operatorname{sgn}(\gamma)f'(u) \ge \delta > 0$ holds. This observation also leads to a local well-posedness result like that in Lemma 2.1. But note that Lemma 2.1 requires $f$ to be smooth only, without any assumption on the sign of $f'$.

### 2.2.2 Riemann invariants

We explore the mathematical structure of the system (2.7) and prove in this section that the system (2.7) possesses an entire set of linearly independent Riemann invariants, which allows us to convert the system (2.7) into an equivalent diagonal form as well.

**Lemma 2.2** (Riemann invariants for the system (2.7)). The system (2.7) possesses a full set of Riemann invariants, which form a coordinate system. The Riemann invariants $R_1 = R_1(\mathbf{U}), \dots, R_4 =$
$R_4(\mathbf{U})$ are given by

$$R_1(\mathbf{U}) = w + v_{\beta} p, \quad R_2(\mathbf{U}) = \alpha u - \psi, \quad R_3(\mathbf{U}) = w - v_{\beta} p, \quad R_4(\mathbf{U}) = \operatorname{sgn}(\gamma)f(u) + \psi.$$

*Proof.* It suffices to show that the row vectors $\ell_i$ defined as

$$\ell_i = \ell_i(\mathbf{U}) = (R_{i,u}(\mathbf{U}), R_{i,\psi}(\mathbf{U}), R_{i,w}(\mathbf{U}), R_{i,p}(\mathbf{U})) \quad (i = 1, \dots, 4, \mathbf{U} \in \mathbb{R}^4),$$

are each a left eigenvector associated with $\lambda_i$. A direct computation gives

$$\ell_1 = (0, 0, 1, v_{\beta}), \quad \ell_2 = (\alpha, -1, 0, 0), \quad \ell_3 = (0, 0, 1, -v_{\beta}), \quad \ell_4 = (\operatorname{sgn}(\gamma)f'(u), 1, 0, 0),$$

and straightforward algebra shows that

$$\ell_i \mathbf{D}f(\mathbf{U}) = \lambda_i \ell_i \quad (i = 1, \dots, 4),$$

which establishes the result.

<figure><img src="image_3.png" /></figure>

Note that the full set of explicit Riemann invariants allows us to make the system completely diagonal along the characteristic directions. This allows, e.g., to find exact solutions of the Riemann problem for the hyperbolic system (2.7) or to devise proper boundary conditions (2.7) beyond our periodic choice.

### 2.2.3 Entropy/entropy-flux pair and weak entropy solutions of (2.7)

In Section 2.2.1, we have already proven that the system (2.7) is a Friedrichs-symmetrizable system and that the associated initial value problem is locally well-posed for classical solutions. Since (2.7) is a system of nonlinear hyperbolic balance laws, we also require a notion of weak solvability. In fact, we show that (2.7) is equipped with an entropy/entropy-flux pair if the flux $f$ is monotone. Then we can naturally define weak entropy solutions, see, e.g., [1].

Let us first collect some assumptions on the behavior of $f$ and some growth conditions needed later on in Section 3.

**Assumption A.** Let $L > 0$ be a constant such that the smooth function $f$ from (1.1) satisfies

$$|f'(u)|, |f''(u)| \le L \quad (u \in \mathbb{R}).$$

Further, there is a constant $\delta > 0$ such that the lower bound $\operatorname{sgn}(\gamma)f'(u) \ge \delta > 0$ holds for all $u \in \mathbb{R}$.

Note that the latter condition ensures that the function

$$\operatorname{sgn}(\gamma)F(u) = \operatorname{sgn}(\gamma) \int_{0}^{u} f(v) dv$$

is uniformly convex in $\mathbb{R}$. It is satisfied for the choice $|f(u)| = au^{2n-1} + bu$, ($n \in \mathbb{N}$) with $a, b > 0$ which has been used in, e.g., the seminal works [4,38].

Under the Assumption A, we can provide an entropy/entropy-flux pair for the system (2.7). By an entropy/entropy-flux pair for the $\alpha$-parameterized system (2.7), we mean a pair $(E^{\alpha}, Q^{\alpha}) : \mathbb{R}^4 \to \mathbb{R}^2$ such that $E^{\alpha}$ is convex with positive definite Hessian $\nabla_{\mathbf{U}}^2 E^{\alpha}$ in $\mathbb{R}^4$ and such that the pair $(E^{\alpha}, Q^{\alpha})$ satisfies

$$\nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U})^{\top} \mathbf{D}f(\mathbf{U}) = \nabla_{\mathbf{U}} Q^{\alpha}(\mathbf{U})^{\top} \text{ for all } \mathbf{U} \in \mathbb{R}^4. \quad (2.13)$$

A specific entropy/entropy-flux pair for the system (2.7) is given in the following lemma.

**Lemma 2.3** (Entropy/entropy-flux pair for (2.7)). The pair $(E^{\alpha}, Q^{\alpha})$ with

$$E^{\alpha}(\mathbf{U}) = \operatorname{sgn}(\gamma)f(u) + \frac{|\gamma|}{2}p^2 + \frac{1}{2\alpha}\psi^2 + \frac{\beta}{2}w^2, \quad (2.14)$$

and

$$Q^{\alpha}(\mathbf{U}) = \frac{1}{2}f(u)^2 + f(u)\psi + \operatorname{sgn}(\gamma)\frac{\psi^2}{2} - |\gamma|pw \quad (2.15)$$
is an entropy/entropy-flux pair for (2.7).

*Proof.* In order to prove this lemma, we first compute the gradients of $E^\alpha$ and $Q^\alpha$ as

$$ \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U}) = (\operatorname{sgn}(\gamma)f(u), \psi/\alpha, \beta w, |\gamma|p)^{\top}, $$

$$ \nabla_{\mathbf{U}} Q^{\alpha}(\mathbf{U}) = (f'(u)(f(u) + \psi), f(u) + \operatorname{sgn}(\gamma)\psi, -|\gamma|p, -|\gamma|w)^{\top}. $$

Then, we directly compute

$$ \begin{aligned} \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U})^{\top} \mathbf{D} f(\mathbf{U}) &= (\operatorname{sgn}(\gamma)f(u), \psi/\alpha, \beta w, \gamma p) \begin{pmatrix} \operatorname{sgn}(\gamma)f'(u) & \operatorname{sgn}(\gamma) & 0 & 0 \\ \alpha f'(u) & \alpha \operatorname{sgn}(\gamma) & 0 & 0 \\ 0 & 0 & 0 & -v_{\beta}^2 \\ 0 & 0 & -1 & 0 \end{pmatrix} \\ &= (f'(u)(f(u) + \psi), f(u) + \operatorname{sgn}(\gamma)\psi, -|\gamma|p, -|\gamma|w) \\ &= \nabla_{\mathbf{U}} Q^{\alpha}(\mathbf{U})^{\top}. \end{aligned} $$

Hence, the pair $(E^\alpha, Q^\alpha)$ satisfies the condition (2.13). Due to the Assumption A, the entropy $E^\alpha$ is strictly convex, and the lemma is proven. □

In passing, we note that the entropy/entropy pair $(E^\alpha, Q^\alpha)$ is asymptotically consistent with the pair $(E, Q)$ for (2.1) from (2.3) in the sense that $(E^\alpha, Q^\alpha)$ collapses to $(E, Q)$ when we let $\alpha \to \infty$. One can readily check that the balance term $\mathbf{b}$ in (2.9) satisfies

$$ \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U}) \cdot \mathbf{b}(\mathbf{U}) = 0 $$

for all $\mathbf{U} \in \mathbb{R}^4$. This property and the compatibility condition (2.13) imply that any smooth solution $\mathbf{U}$ of (2.7) is entropy conservative, i.e.,

$$ E^{\alpha}(\mathbf{U})_t + Q^{\alpha}(\mathbf{U})_x = \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U}) \cdot \mathbf{b}(\mathbf{U}) = 0 \text{ in } I_T. $$

In other words, the balance term does not dissipate the energy $E^\alpha$. For this reason, we cannot directly apply any standard results like [39] to deduce global well-posedness of classical solutions (for at least appropriate initial data). Therefore, we have to account for discontinuous solutions and introduce a weak solution concept for the hyperbolic system (2.7).

**Definition 2.2.** We call a function $\mathbf{U} \in L_{\text{loc}}^{\infty}(I_T; \mathbb{R}^4)$ a weak solution of the system (2.7) with initial data $\mathbf{U}_0 \in L^{\infty}(I; \mathbb{R}^4)$, if we have for each vector-valued function $\varphi \in C_0^{\text{Lip}}(I_T; \mathbb{R}^4)$ the identity

$$ \int_0^T \int_I \left( \mathbf{U} \cdot \varphi_t + \mathbf{f}(\mathbf{U}) \cdot \varphi_x + \mathbf{b}(\mathbf{U}) \cdot \varphi \right) dx dt + \int_I \mathbf{U}_0 \cdot \varphi(\cdot, 0) dx = 0. $$

In view of explicit entropy/entropy-flux pairs of (2.7), we define weak entropy solutions of (2.7) as follows.

**Definition 2.3.** We call a function $\mathbf{U} \in L_{\text{loc}}^{\infty}(I_T; \mathbb{R}^4)$ a weak entropy solution of the system (2.7) associated with the entropy (2.14) if it is a weak solution and if

$$ \int_0^T \int_I \left( E^\alpha(\mathbf{U})\varphi_t + Q^\alpha(\mathbf{U})\varphi_x + \nabla_{\mathbf{U}} E^\alpha(\mathbf{U}) \cdot \mathbf{b}(\mathbf{U})\varphi \right) dx dt + \int_I E^\alpha(\mathbf{U}_0)\varphi(\cdot, 0) dx \geq 0. \quad (2.16) $$

holds for all non-negative functions $\varphi \in C_0^{\text{Lip}}(I_T)$.

### 2.2.4 Least action principle and energy conservation: Hamiltonian structure of the system (2.7)

In this section, we use the principle of classical mechanics to show that the entropy of the hyperbolic system (2.7) is a Hamiltonian with respect to an augmented Lagrangian. In particular, the entropy (2.14) can then be interpreted as the energy of the approximate system (2.7).
First, we recall that the dispersive equation (2.1) is derivable from a least action principle (see §11.7 in [7]), as it can be considered as a generalized KdV equation. Indeed we consider the scalar potential $\eta = \eta(x, t)$ such that $u(x, t) = \eta_x(x, t)$ and define the Lagrangian $\mathcal{L}$ as

$$ \mathcal{L}[\eta_t, \eta_x, \eta_{xx}] := \int_I \left( \operatorname{sgn}(\gamma) \left( \frac{1}{2}\eta_x\eta_t + F(\eta_x) \right) + \frac{|\gamma|}{2} (\eta_{xx})^2 \right) dx, \quad (2.17) $$

and then apply Hamilton's principle for $t_1, t_2 \in [0, T]$ to the action

$$ S[\eta_t, \eta_x, \eta_{xx}] = \int_{t_0}^{t_1} \mathcal{L}[\eta_t, \eta_x, \eta_{xx}] dt. \quad (2.18) $$

This yields the Euler-Lagrange equation

$$ \left(\frac{\partial \mathcal{L}}{\partial \eta_t}\right)_t + \left(\frac{\partial \mathcal{L}}{\partial \eta_x}\right)_x - \left(\frac{\partial \mathcal{L}}{\partial \eta_{xx}}\right)_{xx} = \operatorname{sgn}(\gamma) (\eta_{xt} + (f(\eta_x))_x) - |\gamma|\eta_{xxx} = 0, \quad (2.19) $$

which corresponds to (2.1) after substituting $u = \eta_x$.

We point out that the analogous least action principle can be obtained for the hyperbolic approximation in its wave propagation form (2.4). We prove this in the following lemma.

**Proposition 2.4.** The hyperbolic system (2.7) can be formulated as a Hamiltonian system with an augmented Lagrangian of the form

$$ \mathcal{L}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x] := \int_I \left( \operatorname{sgn}(\gamma) \left( \frac{1}{2}\eta_x\eta_t + F(\eta_x) \right) + \frac{|\gamma|c_x^2}{2} + \frac{\alpha}{2}(\eta_x - c)^2 - \frac{\beta}{2}c_t^2 \right) dx. \quad (2.20) $$

Moreover, the Hamiltonian density with respect to the augmented Lagrangian (2.20) satisfies

$$ \mathcal{H}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x] = - \int_I E^{\alpha}(\mathbf{U}) dx, \quad (2.21) $$

where $E^{\alpha}(\mathbf{U})$ is the entropy of the hyperbolic system (2.7) defined in (2.14).

*Proof.* Indeed, let us consider the so-called augmented Lagrangian

$$ \mathcal{L}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x] = \int_I \left( \operatorname{sgn}(\gamma) \left( \frac{1}{2}\eta_x\eta_t + F(\eta_x) \right) + \frac{|\gamma|c_x^2}{2} + \frac{\alpha}{2}(\eta_x - c)^2 - \frac{\beta}{2}c_t^2 \right) dx, $$

which is built from the Lagrangian (2.17) by substituting $\eta_{xx}$ by $c_x$ and adding two relaxation terms proportional to $\alpha$ and $\beta$, respectively. Then, applying Hamilton's principle to the Lagrangian (2.20), under the variations of $\eta$ and $c$ respectively, yields the two Euler-Lagrange equations

$$ \left(\frac{\partial \mathcal{L}_{\alpha}}{\partial \eta_t}\right)_t + \left(\frac{\partial \mathcal{L}_{\alpha}}{\partial \eta_x}\right)_x = \operatorname{sgn}(\gamma) (\eta_{xt} + f(\eta_x)_x) + \alpha(\eta_x - c)_x = 0, \quad (2.22a) $$

$$ -\frac{\partial \mathcal{L}_{\alpha}}{\partial c} + \left(\frac{\partial \mathcal{L}_{\alpha}}{\partial c_t}\right)_t + \left(\frac{\partial \mathcal{L}_{\alpha}}{\partial c_x}\right)_x = -\beta c_{tt} + |\gamma|c_{xx} + \alpha(\eta_x - c) = 0, \quad (2.22b) $$

which again, after substituting $\eta_x = u$ reduce to the system (2.4).

Now the Hamiltonian density $\mathcal{H}_{\alpha}$ is defined by the partial Legendre transform of the Lagrangian density, and is given by

$$ \mathcal{H}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x] = \frac{\partial \mathcal{L}_{\alpha}}{\partial \eta_t} \eta_t + \frac{\partial \mathcal{L}_{\alpha}}{\partial c_t} c_t - \mathcal{L}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x]. $$

Using (2.20), $\eta_x = u$ and in view of (2.14), we have

$$ \mathcal{H}_{\alpha}[\eta_t, \eta_x, c, c_t, c_x] = - \int_I E^{\alpha}(\mathbf{U}) dx. $$

The lemma is then proven.

Lemma 2.4 shows that the system (2.7) also preserves the Hamiltonian and Lagrangian structure of the dispersive equation (2.1) as $\alpha \to \infty$, which is not possible in many relaxation systems. In
particular, a Lagrangian interpretation is not available for the parabolic approximation from the work of Corli&Rohde [22]. Therefore, we consider the hyperbolic approximation (2.7), the more consistent choice if it is about the approximation of the purely dispersive equation (2.1).

## 2.3 Second-order approximations of the diffusive-dispersive equation (1.1)

We focus on the complete diffusive-dispersive equation, that is (1.1) with $\varepsilon > 0$ and $\gamma \in \mathbb{R} \setminus \{0\}$. Thus, we consider the equation

$$u_t + f(u)_x = \varepsilon u_{xx} + \gamma u_{xxx}, \quad \text{in } I_T, \qquad (2.23)$$

subject to periodic boundary conditions and the initial condition

$$u(\cdot, 0) = u_0 \text{ in } I. \qquad (2.24)$$

The energy structure of the (2.23) is similar to the energy as defined in (1.4). Moreover, classical results for local and global well-posedness of classical solutions are available in $H^s$-spaces. In particular, we refer the interested reader to the work of Gallego&Pazato [40] for a systematic treatment of generalized KdV-Burgers equations.

To approximate the diffusive-dispersive equation (2.23), we include a viscoelastic damping operator [41] instead of the pure wave operator in (2.4). Precisely, we propose the following approximate system in $I_T$ for the diffusive-dispersive equation (2.23):

$$\begin{gathered} u_t + f(u)_x + \operatorname{sgn}(\gamma)\alpha(u-c)_x = \varepsilon u_{xx}, \qquad (2.25) \\ \beta(c_{tt} - \varepsilon c_{xxt}) - |\gamma|c_{xx} = \alpha(u-c) \end{gathered}$$

With the change of variables

$$w = c_t - \varepsilon c_{xx}, \quad p = c_x, \quad \psi = \alpha(u - c),$$

one can easily convert the system (2.25) into the mixed- but lower-order system

$$\begin{gathered} u_t + (f(u) + \operatorname{sgn}(\gamma)\psi)_x = \varepsilon u_{xx}, \\ \frac{1}{\alpha}\psi_t + (f(u) + \operatorname{sgn}(\gamma)\psi)_x = -w + \frac{\varepsilon}{\alpha}\psi_{xx}, \quad \text{in } I_T. \qquad (2.26) \\ \beta w_t - |\gamma|p_x = \psi, \\ p_t - w_x = \varepsilon p_{xx} \end{gathered}$$

The system (2.26) is complemented with the prepared initial data

$$u^{\alpha}(\cdot, 0) = u_0, \quad \psi^{\alpha}(\cdot, 0) = -|\gamma|u_{0,xx}, \quad w^{\alpha}(\cdot, 0) = -f(u_0)_x - |\gamma|u_{0,xxx}, \quad p^{\alpha}(\cdot, 0) = u_{0,x} \text{ in } I. \qquad (2.27)$$

In what follows, we prove that the entropy/entropy-flux pairs for the first-order system (2.7) are compatible with the system (2.26) as well.

## 2.4 Qualitative properties of the second-order approximation

### 2.4.1 Hyperbolic-parabolic structure, energy dissipation, and local well-posedness

Let us again omit in this section the superscript $(\cdot)^\alpha$ on the unknowns and write $U = (u, \psi, w, p)^T$ for $U^\alpha = (u^\alpha, \psi^\alpha, w^\alpha, p^\alpha)^T \in \mathbb{R}^4$.

With the notations from (2.9) and the positive semi-definite matrix

$$\mathbf{D} = \begin{pmatrix} \varepsilon & 0 & 0 & 0 \\ 0 & \varepsilon & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \varepsilon \end{pmatrix}$$
we write the system (2.26) for some smooth solution $\mathbf{U}$ in the form

$$ \mathbf{U}_t + \mathbf{f}(\mathbf{U})_x = \mathbf{S}\mathbf{U} + \mathbf{D}\mathbf{U}_{xx} \iff \mathbf{U}_t + \mathbf{D}\mathbf{f}(\mathbf{U})\mathbf{U}_x = \mathbf{S}\mathbf{U} + \mathbf{D}\mathbf{U}_{xx}. \quad (2.28) $$

From (2.28) we see that (2.26) is a (degenerate) hyperbolic-parabolic system. In the next step we show that for Assumption A imposed, smooth solutions of (2.26) dissipate the entropy $E^\alpha$ from (2.14). Recall that smooth solutions of the hyperbolic system (2.7) conserve $E^\alpha$.

**Lemma 2.4** (Entropy dissipation for (2.26)). Let Assumption A hold and let $\mathbf{U}$ be a smooth solution of the hyperbolic-parabolic system (2.26) in $I_T$.

Then we have for $E^\alpha$ from (2.13) in $(0,T)$ the estimate

$$ \begin{aligned} \frac{d}{dt} \int_I E^\alpha(\mathbf{U}) dx &= \frac{d}{dt} \int_I \left( \operatorname{sgn}(\gamma) F(u) + \frac{|\gamma|p^2}{2} + \frac{\psi^2}{2\alpha} + \frac{\beta w^2}{2} \right) dx \\ &= -\varepsilon \int_I \left( \operatorname{sgn}(\gamma) f'(u) u_x^2 + |\gamma| p_x^2 + \frac{\psi_x^2}{\alpha} \right) dx \quad (2.29) \\ &\le 0. \end{aligned} $$

*Proof.* We directly take the time derivative of $E^\alpha(\mathbf{U})$ and compute

$$ \frac{d}{dt} E^\alpha(\mathbf{U}) = \operatorname{sgn}(\gamma) f(u) u_t + \beta w t + \frac{\psi}{\alpha} \psi_t + |\gamma| p t. $$

Then, in view of (2.26), we have

$$ \frac{d}{dt} E^\alpha(\mathbf{U}) = -\frac{d}{dx} \left( \frac{f(u)^2}{2} + f(u)\psi + \operatorname{sgn}(\gamma) \frac{\psi^2}{2} - |\gamma|pw \right) + \varepsilon \left( \operatorname{sgn}(\gamma) f(u) u_{xx} + |\gamma| p p_{xx} + \frac{\psi}{\alpha} \psi_{xx} \right). $$

Integrating both sides with respect to $x$ and using the periodic boundary conditions, we get the statement of the lemma. $\square$

In view of Assumption A, the entropy in (2.14) is strictly convex and therefore the local well-posedness of the hyperbolic-parabolic system (2.26) for the class of flux functions satisfying Assumption A can be obtained using the classical theory of symmetrizable hyperbolic-parabolic systems due to Kato [42] or Kawashima [43]. To be precise, we have the following.

**Lemma 2.5** (Local well-posedness of the initial value problem for system (2.26)). Let the function $\mathbf{U}_0 = (u_0, \psi_0, w_0, p_0)^\top$ be given such that

$$ \mathbf{U}_0 \in (H^s(I))^4, \quad s > 3/2, $$

holds. Suppose that the flux $f$ satisfies Assumption A.

Then, there exists for all $\alpha, \beta > 0$ and $\gamma \in \mathbb{R} \setminus \{0\}$ a finite time $T^* \in (0, \infty)$ depending on $\|\mathbf{U}_0\|_{(H^s(I))^4}$, such that the initial value problem for (2.26) with $\mathbf{U}(\cdot, 0) = \mathbf{U}_0$ admits a unique classical solution $\mathbf{U}$, which satisfies

$$ \mathbf{U}(\cdot, t) \in C([0, T^*]; (H^s(I))^4) \cap C^1([0, T^*]; (H^{s-1}(I))^4), \quad s > 3/2. $$

**Remark 2.6.** 1. Note that the entropy dissipation for our system (2.26) in Lemma 2.4 coincides with the energy dissipation for the elliptic approximation (1.5) (see Lemma 3.3 in [21]). In particular, this shows that the proposed system (2.26) does not introduce any additional diffusion into the dynamics. In contrast, the parabolic approximation in [22] exhibits stronger diffusive effects than the proposed system (2.26) (see Lemma 3.5 in [22]).

2. The condition on the convexity of $\operatorname{sgn}(\gamma)F(u)$ for the diffusive-dispersive case can be further relaxed. Assume that $\operatorname{sgn}(\gamma)f'(u) \ge -\bar{f}$ for some $\bar{f} > 0$. Then one can choose a number $\kappa$ such that $\kappa + \bar{f} > \delta > 0$. In that case, we can have a similar energy inequality with a shift in $\operatorname{sgn}(\gamma)F(u)$ given by $\operatorname{sgn}(\gamma)F(u) + \frac{\kappa u^2}{2}$, which still makes the energy convex. However, this
can only work for the diffusive-dispersive equation, as it gives enough dissipation to absorb the residual terms for a fixed $\epsilon > 0$.

3. Note that the system (2.7) preserves the energy of the limit equation (2.23) regardless of the choice of flux functions. In this work, we restrict ourselves to the class of flux functions satisfying Assumption A, as the convergence analysis in Section 3 relies on the uniform convexity of the entropy $E_\alpha$. Nevertheless, the numerical examples in Section 4 indicate that the approximate system (2.7) provides accurate approximations even for equations, where the flux function does not satisfy Assumption A.

### 2.4.2 Global well-posedness of classical solutions of (2.26) for prepared initial data

In this section, we prove that the local solutions developed in Lemma 2.5 can be extended to $(0, \infty)$ as long as the initial data is in a sufficiently small neighborhood of a constant equilibrium point $U^* = (u^*, 0, 0, 0)^\top$. In particular, we prove that the system (2.26) belongs to the class of partially diffusive systems satisfying the Shizuta-Kawashima (SK) condition [43, 44]. Precisely, we prove the following.

**Theorem 2.7.** Let $\alpha, \beta > 0$ and $\gamma \in \mathbb{R} \setminus \{0\}$ be given. Suppose that the flux $f$ satisfies Assumption A. Then there exists a $\Lambda > 0$ such that for all $U_0 = (u_0, \psi_0, w_0, p_0)^\top \in (H^s(I))^4$, $s \ge 2$, $s \in \mathbb{N}$, with

$$ \left( \|u_0 - u^*\|_{H^s(I)}^2 + \frac{1}{\alpha} \|\psi_0\|_{H^s(I)}^2 + \beta \|w_0\|_{H^s(I)}^2 + |\gamma| \|p_0\|_{H^s(I)}^2 \right) < \Lambda. $$

the initial value problem for (2.26) has a unique classical solution $U$ with $U(\cdot, 0) = U_0$. It satisfies

$$ U - U^* \in C([0, \infty); (H^s(I))^4) \cap C^1([0, \infty); (H^{s-1}(I))^4). \quad (2.30) $$

Furthermore, we have for a positive constant $C = C(\Lambda)$ the estimate

$$ \begin{aligned} & \sup_{0 \le t < \infty} \left( \|u - u^*\|_{H^s(I)}^2 + \frac{1}{\alpha} \|\psi\|_{H^s(I)}^2 + \beta \|w\|_{H^s(I)}^2 + |\gamma| \|p\|_{H^s(I)}^2 \right) \\ & + \varepsilon \int_0^t \left( \|u - u^*\|_{H^{s+1}(I)} + |\gamma| \|p\|_{H^{s+1}(I)} + \frac{1}{\alpha} \|\psi\|_{H^{s+1}(I)} \right) d\tau \quad (2.31) \\ & \le C(\Lambda) \left( \|u_0 - u^*\|_{H^s(I)}^2 + \frac{1}{\alpha} \|\psi_0\|_{H^s(I)}^2 + \beta \|w_0\|_{H^s(I)}^2 + |\gamma| \|p_0\|_{H^s(I)}^2 \right). \end{aligned} $$

*Proof.* Since the system (2.26) is endowed with an entropy/entropy-flux pair, it can be converted into a symmetric hyperbolic-parabolic form by using the Hessian of the entropy $E^\alpha(U)$ as its symmetrizer. By multiplying the equation (2.28) with $\nabla^2 U E^\alpha(U)$, we obtain the symmetric form of (2.28) as

$$ S_0(U)U_t + S_1(U)U_x - S_2(U)U_{xx} = S_3(U)U, \quad (2.32) $$

where

$$ S_0(U) = \nabla_U^2 E^\alpha(U) = \begin{pmatrix} \operatorname{sgn}(\gamma) f'(u) & 0 & 0 & 0 \\ 0 & 1/\alpha & 0 & 0 \\ 0 & 0 & \beta & 0 \\ 0 & 0 & 0 & |\gamma| \end{pmatrix}, $$

$$ S_1(U) = \nabla_U^2 E^\alpha(U) D f(U) = \begin{pmatrix} (f'(u))^2 & f'(u) & 0 & 0 \\ f'(u) & \operatorname{sgn}(\gamma) & 0 & 0 \\ 0 & 0 & 0 & -|\gamma| \\ 0 & 0 & -|\gamma| & 0 \end{pmatrix}, $$
$$ \mathbf{S}_2(\mathbf{U}) = \nabla_{\mathbf{U}}^2 E^{\alpha}(\mathbf{U})\mathbf{D} = \begin{pmatrix} \varepsilon \operatorname{sgn}(\gamma) f'(u) & 0 & 0 & 0 \\ 0 & \varepsilon/\alpha & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & |\gamma|\varepsilon \end{pmatrix}, $$

$$ \mathbf{S}_3(\mathbf{U}) = \nabla_{\mathbf{U}}^2 E^{\alpha}(\mathbf{U})\mathbf{S} = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. $$

Clearly, the system (2.32) belongs to the general class of symmetric partially diffusive hyperbolic-parabolic systems with a skew-symmetric source matrix.

We first verify that the flux and the diffusion matrix of the system (2.32) satisfy the Shizuta-Kawashima (SK) conditions [43,44]. The kernel of the diffusion matrix $\mathbf{S}_2(\mathbf{U})$ is given by

$$ \ker \mathbf{S}_2(\mathbf{U}) = \{(0, 0, w, 0)^\top : w \in \mathbb{R}\} = \operatorname{span} \{(0, 0, 1, 0)^\top\}. $$

A straightforward computation shows that no eigenvectors of $\mathbf{S}_1(\mathbf{U})$ lies in $\ker \mathbf{S}_2(\mathbf{U})$ and thus we have for all eigenvalues $\lambda$ of $\mathbf{S}_1(\mathbf{U})$

$$ \ker \mathbf{S}_2(\mathbf{U}) \cap \ker(\mathbf{S}_1(\mathbf{U}) - \lambda \mathbf{I}) = \{0\}. $$

This shows that the system (2.26) or equivalently (2.28) belongs to the partially diffusive systems satisfying the SK coupling conditions. Moreover, the source term does not contribute to the higher-order energy dissipation (or increment). To be precise, it is easy to see that the source term vanishes after multiplying by energy multipliers for $H^s$-estimates ($s = 0, 1, 2, \dots$), or in other words

$$ \partial_x^s (\mathbf{U} - \mathbf{U}^*)^\top \mathbf{S}_3(\mathbf{U}) \partial_x^s (\mathbf{U} - \mathbf{U}^*) = 0, \quad s = 0, 1, 2, \dots. $$

This, in particular, allows us to follow the classical theory for the global well-posedness of smooth solutions of partially diffusive systems (see e.g. Theorem 4.3 in [43]). Therefore, the results of the theorem are a direct consequence of the classical theory due to Kawashima [43]. □

**Remark 2.8.** There are recent results on the local and global well-posedness for partially diffusive systems where the initial data need less regularity and, e.g., belong to some Besov space instead. We refer the interested reader to the recent work of Adogbo&Danchin [45] for more details. However, for our analysis, we work with $H^s$-spaces.

## 3 Hyperbolic approximations and asymptotics

In this section, we first apply a Chapman-Enskog expansion on the variables of the system (2.4) and (2.25) to provide a formal justification of the consistency of the proposed system(s) as $\alpha \to \infty$. The Chapman-Enskog expansion provides an asymptotic expansion for the solutions of (2.4) and (2.25) and identifies (2.1) and (2.23) as the limit equation(s), respectively. However, this argument is not rigorous at the level of entropy solutions. To obtain a rigorous convergence result, we complement the formal expansion with a relative entropy framework to the systems (2.7) and (2.26). Precisely, we prove the convergence of the solutions of approximate systems to the classical solutions of (2.23) for $\varepsilon \ge 0$ and the flux functions satisfying Assumption A.

Throughout this section, we denote the solutions of the approximate systems (2.7) and (2.26) as $\mathbf{U}^\alpha = (u^\alpha, \psi^\alpha, w^\alpha, p^\alpha)^\top$ and the solution of the diffusive-dispersive equation (2.23) by $\mathbf{U} = (u, \psi, w, p)^\top = (u, -|\gamma|u_{xx}, u_t, u_x)^\top$ for $\varepsilon \ge 0$. Based on the proposed systems (2.7) and (2.26), we divide our convergence analysis into two different sections.
## 3.1 Dispersive asymptotics for the first-order hyperbolic system (2.7)

We discuss the asymptotics for $\alpha \to \infty$ of the sequence $\{\mathbf{U}^\alpha\}_{\alpha>0}$ solving (2.7) for prepared initial data (2.27) and periodic boundary conditions. We will show in Theorem 3.1 that the sequence converges to a solution of the dispersive equation (2.1) subject to the initial datum from (2.2).

To motivate a scaling of $\beta$ in terms of the parameter $\alpha$ we consider Chapman-Enskog expansions. Precisely we choose $\beta = \alpha^{-1}$ and assume that the solutions of the system (2.4) admit a smooth expansion in the parameter $\alpha$ such that

$$ 
\begin{aligned}
u^{\alpha}(x,t) &= u^{0}(x,t) + \alpha^{-1}u^{1}(x,t) + \alpha^{-2}u^{2}(x,t) + \mathcal{O}(\alpha^{-3}), \\
c^{\alpha}(x,t) &= c^{0}(x,t) + \alpha^{-1}c^{1}(x,t) + \alpha^{-2}c^{2}(x,t) + \mathcal{O}(\alpha^{-3}).
\end{aligned}
 $$

Inserting these expansions in (2.4)$_2$ allows to write

$$ \alpha(u^0 - c^0) + (u^1 - c^1 + |\gamma|c_{xx}^0) + \alpha^{-1}(u^2 - c^2 - c_{tt}^0 + |\gamma|c_{xx}^1) = \mathcal{O}(\alpha^{-2}). \quad (3.1) $$

Then, successively equating the series coefficients to zero in (3.1) leads to

$$ u^0 - c^0 = 0, \quad u^1 - c^1 = -|\gamma|c_{xx}^0, \quad u^2 - c^2 = c_{tt}^0 - |\gamma|c_{xx}^1, \quad (3.2) $$

and therefore

$$ 
\begin{aligned}
\alpha(u^{\alpha} - c^{\alpha}) &= -|\gamma|c_{xx}^{0} + \alpha^{-1}(c_{tt}^{0} - |\gamma|c_{xx}^{1}) + \mathcal{O}(\alpha^{-2}) \\
&= -|\gamma|u_{xx}^{0} + \alpha^{-1}(u_{tt}^{0} - |\gamma|(u^{1} + |\gamma|u_{xx}^{0})_{xx}) + \mathcal{O}(\alpha^{-2}) \\
&= -|\gamma|(u^{0} + \alpha^{-1}u^{1})_{xx} + \alpha^{-1}(u_{tt}^{0} - \gamma^{2}u_{xxxx}^{0} + \mathcal{O}(\alpha^{-2})) \\
&= -|\gamma|u_{xx}^{\alpha} + \alpha^{-1}(u_{tt}^{\alpha} - \gamma^{2}u_{xxxx}^{\alpha}) + \mathcal{O}(\alpha^{-2}).
\end{aligned}
 $$

Using the last relation in (2.4)$_1$ finally yields

$$ u_t^\alpha + f(u^\alpha)_x = \operatorname{sgn}(\gamma)|\gamma|u_{xxx}^\alpha + \alpha^{-1}(\gamma^2 u_{xxxxx}^\alpha - u_{ttx}^\alpha) + \mathcal{O}(\alpha^{-2}). \quad (3.3) $$

Thus, given the suggested scaling $\beta = \alpha^{-1}$, the leading order error term with respect to the dispersive equation (2.1) is first-order in $\alpha^{-1}$.

In what follows, we prove that weak entropy solutions of (2.7)-(2.8) converge to the solutions of the dispersive equation (2.1) when the parameter $\alpha \to \infty$ as our first main result.

**Theorem 3.1** (Main result I). Let $u_0 \in H^6(I)$, and suppose that Assumption A holds. We assume that $u \in L^2(0,T; H^3(I))$ with $u_{tt} \in L^2(I_T)$ and $\|u\|_{L^\infty(0,T; H^4(I))} < \infty$ is a classical solution of (2.1), (2.2). Consider for $\beta = \alpha^{-1}$, a sequence $\{\mathbf{U}^\alpha = (u^\alpha, \psi^\alpha, w^\alpha, p^\alpha)\}_{\alpha>0}$ of weak entropy solutions of (2.7), (2.8) in the sense of Definition 2.3.

Then there is a constant $C = C(L, \delta, u) \ge 0$ that is independent of $\alpha$ such that

$$ \|u^{\alpha} - u\|_{L^{\infty}(0,T;L^2(I))} + \|p^{\alpha} - u_x\|_{L^{\infty}(0,T;L^2(I))} \le C\sqrt{\frac{1}{\alpha}} \quad (3.4) $$

holds.

*Proof.* With $\psi = -|\gamma|u_{xx}$, $w = u_t$ and $p = u_x$, we rewrite the dispersive limit equation (2.1) in terms of the approximate system (2.7), which leads to the system

$$ 
\begin{aligned}
& u_t + f(u)_x + \operatorname{sgn}(\gamma)\psi_x = 0, \\
& \frac{1}{\alpha}\psi_t + f(u)_x + \operatorname{sgn}(\gamma)\psi_x = -w - \frac{1}{\alpha}|\gamma|u_{xxt}, \\
& \beta w_t - |\gamma|p_x = \psi + \beta u_{tt}, \\
& p_t - w_x = 0.
\end{aligned}
 \quad (3.5) $$

Define $\mathbf{U} = (u, \psi, w, p)^\top$ and fix the initial datum $\mathbf{U}_0 = (u_0, -|\gamma|u_{0,xx}, -f(u_0)_x - |\gamma|u_{0,xxx}, u_{0,x})^\top$. Since $u$ is a classical solution of the initial value problem for (2.1) with $u(\cdot, 0) = u_0$, the function $\mathbf{U}$ is
a classical solution of (3.5) with $\mathbf{U}(\cdot, 0) = \mathbf{U}_0$.

It is then straightforward to check that $\mathbf{U}$ satisfies for the entropy/entropy-flux pair $(E^\alpha, Q^\alpha)$ from Lemma 2.3 the integral identity

$$ \int_0^T \int_I (E^\alpha(\mathbf{U})\varphi_t + Q^\alpha(\mathbf{U})\varphi_x + \nabla_{\mathbf{U}}E^\alpha(\mathbf{U}) \cdot \bar{\mathbf{b}}(\mathbf{U})\varphi) \, dx \, dt + \int_I E^\alpha(\mathbf{U}_0)\varphi(\cdot, 0) \, dx = 0 \quad (3.6) $$

for all $\varphi \in C_0^{\text{Lip}}(I \times [0, T))$, $\varphi \ge 0$. In (3.6) we used $\bar{\mathbf{b}}(\mathbf{U}) = (0, -\alpha w - |\gamma|u_{xxt}, \beta^{-1}\psi + u_{tt}, 0)^\top$. Furthermore, since $\mathbf{U}^\alpha$ is a weak entropy solution of (2.7), (2.8), we have

$$ \int_0^T \int_I (E^\alpha(\mathbf{U}^\alpha)\varphi_t + Q^\alpha(\mathbf{U}^\alpha)\varphi_x + \nabla_{\mathbf{U}}E^\alpha(\mathbf{U}^\alpha) \cdot \mathbf{b}(\mathbf{U}^\alpha)\varphi) \, dx \, dt + \int_I E^\alpha(\mathbf{U}_0^\alpha)\varphi(\cdot, 0) \, dx \ge 0. \quad (3.7) $$

Now, we recall that the entropy $E^\alpha$ from (2.14) is almost quadratic and define the associated relative entropy $\mathcal{E}^\alpha$ for the approximate solution $\mathbf{U}^\alpha$ and the function $\mathbf{U}$ according to

$$ \begin{aligned} \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}) &= E^\alpha(\mathbf{U}^\alpha) - E^\alpha(\mathbf{U}) - \nabla_{\mathbf{U}}E^\alpha(\mathbf{U})(\mathbf{U}^\alpha - \mathbf{U}) \\ &= \operatorname{sgn}(\gamma)(F(u^\alpha) - F(u) - f(u)(u^\alpha - u)) + \frac{|\gamma|}{2}(p^\alpha - p)^2 + \frac{\beta}{2}(w^\alpha - w)^2 + \frac{1}{2\alpha}(\psi^\alpha - \psi)^2. \end{aligned} \quad (3.8) $$

Recall that $\operatorname{sgn}(\gamma)F$ is a uniformly convex function by Assumption A. Thus, there exists an $\alpha$-independent constant $C_1 > 0$ such that

$$ C_1 \left( (u^\alpha - u)^2 + \frac{1}{\alpha}(\psi^\alpha - \psi)^2 + \beta(w^\alpha - w)^2 + (p^\alpha - p)^2 \right) \le \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}). \quad (3.9) $$

Now, we subtract (3.6) from (3.7) to obtain

$$ \begin{aligned} & \int_0^T \int_I ((E^\alpha(\mathbf{U}^\alpha) - E^\alpha(\mathbf{U}))\varphi_t + (Q^\alpha(\mathbf{U}^\alpha) - Q^\alpha(\mathbf{U}))\varphi_x) \, dx \, dt \\ & \qquad + \int_0^T \int_I (\nabla_{\mathbf{U}}E^\alpha(\mathbf{U}^\alpha) \cdot \mathbf{b}(\mathbf{U}^\alpha) - \nabla_{\mathbf{U}}E^\alpha(\mathbf{U}) \cdot \mathbf{b}(\mathbf{U})) \, \varphi \, dx \, dt \\ & \qquad + \int_I (E^\alpha(\mathbf{U}_0) - E^\alpha(\mathbf{U}_0)) \, \varphi(x, 0) \, dx \ge 0. \end{aligned} \quad (3.10) $$

Moreover, since $\mathbf{U}$ is a classical solution of (3.5), it must also be a solution of the system (3.5) in the distributional sense. If we choose for $\varphi \in C_0^{\text{Lip}}(I \times [0, T))$ the test function $\varphi = \varphi \nabla_{\mathbf{U}}E^\alpha(\mathbf{U}) \in C_0^{\text{Lip}}(I \times [0, T); \mathbb{R}^4)$ like in Definition 2.2, we obtain from the fact that $\mathbf{U}^\alpha$ is a weak solution, the identity

$$ \begin{aligned} & \int_0^T \int_I (\mathbf{U}^\alpha - \mathbf{U}) \cdot (\varphi \nabla_{\mathbf{U}} E^\alpha(\mathbf{U})_t + (\mathbf{f}(\mathbf{U}^\alpha) - \mathbf{f}(\mathbf{U})) \cdot (\varphi \nabla_{\mathbf{U}} E^\alpha(\mathbf{U}))_x) \, dx \, dt \\ & \qquad + \int_0^T \int_I (\mathbf{b}(\mathbf{U}^\alpha) - \mathbf{b}(\mathbf{U})) \cdot (\varphi \nabla_{\mathbf{U}} E^\alpha(\mathbf{U})) \, dx \, dt + \int_I (\mathbf{U}_0^\alpha - \mathbf{U}_0) \cdot (\varphi \nabla_{\mathbf{U}} E^\alpha(\mathbf{U}))(x, 0) \, dx = 0. \end{aligned} \quad (3.11) $$

From the relations

$$ (\nabla_{\mathbf{U}} E^\alpha(\mathbf{U}))_t = \mathbf{U}_t^\top \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}), \quad \text{and} \quad (\nabla_{\mathbf{U}} E^\alpha(\mathbf{U}))_x = \mathbf{U}_x^\top \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}), $$
we get after subtracting (3.11) from (3.10) the identity

$$
\begin{align}
\mathcal{I} = & \int_0^T \int_I (\mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}))\varphi_t + (\mathcal{Q}^\alpha(\mathbf{U}^\alpha|\mathbf{U}))\varphi_x \, dx \, dt \nonumber \\
& + \underbrace{\int_0^T \int_I (\mathbf{U}^\alpha - \mathbf{U})^\top \mathbf{U}_t^\top \nabla_\mathbf{U}^2 E^\alpha(\mathbf{U}) + (\mathbf{f}(\mathbf{U}^\alpha) - \mathbf{f}(\mathbf{U}))^\top \mathbf{U}_x^\top \nabla_\mathbf{U}^2 E^\alpha(\mathbf{U})) \, \varphi \, dx \, dt}_{=: \mathcal{I}_1} \tag{3.12} \\
& + \underbrace{\int_0^T \int_I (\nabla_\mathbf{U} E^\alpha(\mathbf{U}^\alpha) - \nabla_\mathbf{U} E^\alpha(\mathbf{U})) \cdot \mathbf{b}(\mathbf{U}^\alpha) \, \varphi \, dx \, dt}_{=: \mathcal{I}_2} + \int_I \mathcal{E}^\alpha(\mathbf{U}_0|\mathbf{U}_0)\varphi(x, 0) \, dx \ge 0. \nonumber
\end{align}
$$

Here we used the relative entropy flux

$$
\mathcal{Q}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U}) = \mathcal{Q}^{\alpha}(\mathbf{U}^{\alpha}) - \mathcal{Q}^{\alpha}(\mathbf{U}) - \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U}) \cdot (\mathbf{f}(\mathbf{U}^{\alpha}) - \mathbf{f}(\mathbf{U})),
$$

which satisfies for an α-independent constant $C_2 \ge 0$ the inequality $|\mathcal{Q}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U})| \le C_2|\mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U})|$.

Now we focus on the integrals $\mathcal{I}_1$ and $\mathcal{I}_2$ in (3.12). First, consider the integral $\mathcal{I}_1$, which, in view of (3.5) writes as

$$
\begin{align}
\mathcal{I}_1 &= \int_0^T \int_I \left( (\mathbf{U}_t \cdot (\mathbf{U}^\alpha - \mathbf{U}) + \mathbf{U}_x \cdot (\mathbf{f}(\mathbf{U}^\alpha) - \mathbf{f}(\mathbf{U})))^\top \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}) \right) \varphi \, dx \, dt \nonumber \\
&= \int_0^T \int_I \left( (\mathbf{U}_x \cdot (\mathbf{f}(\mathbf{U}^\alpha) - \mathbf{f}(\mathbf{U}) - \mathbf{D}\mathbf{f}(\mathbf{U})(\mathbf{U}^\alpha - \mathbf{U})))^\top \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}) \right) \varphi \, dx \, dt \nonumber \\
&\quad + \int_0^T \int_I (\mathbf{b}(\mathbf{U}) \cdot (\mathbf{U}^\alpha - \mathbf{U}))^\top \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}) \varphi \, dx \, dt \tag{3.13} \\
&= \int_0^T \int_I \left( \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}) (\mathbf{f}(\mathbf{U}^\alpha) - \mathbf{f}(\mathbf{U}) - \mathbf{D}\mathbf{f}(\mathbf{U})(\mathbf{U}^\alpha - \mathbf{U})) \cdot \mathbf{U}_x \right) \varphi \, dx \, dt \nonumber \\
&\quad + \underbrace{\int_0^T \int_I \nabla_{\mathbf{U}}^2 E^\alpha(\mathbf{U}) (\mathbf{U}^\alpha - \mathbf{U})) \cdot \mathbf{b}(\mathbf{U}) \varphi \, dx \, dt}_{=: \mathcal{I}_3} \nonumber
\end{align}
$$

where we have used the symmetry of the Hessian matrix.

Now in view of (3.5) and the definition of $E^\alpha$, it is easy to obtain that

$$
\begin{align}
& \nabla_{\mathbf{U}}^2 E^{\alpha}(\mathbf{U}) (\mathbf{f}(\mathbf{U}^{\alpha}) - \mathbf{f}(\mathbf{U}) - \mathbf{D}\mathbf{f}(\mathbf{U})(\mathbf{U}^{\alpha} - \mathbf{U})) \cdot \mathbf{U}_x \tag{3.14} \\
& \qquad = \operatorname{sgn}(\gamma) (f(u^{\alpha}) - f(u) - f'(u)(u^{\alpha} - u)) (f'(u)u_x + \psi_x). \nonumber
\end{align}
$$

Moreover, $\mathcal{I}_2$ from (3.12) and $\mathcal{I}_3$ from (3.13) together become

$$
\begin{align}
\mathcal{I}_2 + \mathcal{I}_3 &= \int_0^T \int_I (\nabla_{\mathbf{U}}^2 E^{\alpha}(\mathbf{U})(\mathbf{U}^{\alpha} - \mathbf{U})) \cdot \mathbf{b}(\mathbf{U}) + (\nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U}^{\alpha}) - \nabla_{\mathbf{U}} E^{\alpha}(\mathbf{U})) \cdot \mathbf{b}(\mathbf{U}^{\alpha}) \varphi \, dx \, dt \tag{3.15} \\
&= \int_0^T \int_I \left( \frac{|\gamma|}{\alpha} (\psi^{\alpha} - \psi) u_{xxt} - \beta (w^{\alpha} - w) u_{tt} \right) \varphi \, dx \, dt. \nonumber
\end{align}
$$

Thus, using (3.13), (3.14) and (3.15) in (3.12), the integral $\mathcal{I}$ in (3.12) reduces to

$$
\begin{align}
\mathcal{I} = & \int_0^T \int_I (\mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U})\varphi_t + \mathcal{Q}^\alpha(\mathbf{U}^\alpha|\mathbf{U})\varphi_x) \, dx \, dt \nonumber \\
& + \int_0^T \int_I \operatorname{sgn}(\gamma) (f(u^\alpha) - f(u) - f'(u)(u^\alpha - u)) (f'(u)u_x + \psi_x) \varphi \, dx \, dt \tag{3.16} \\
& + \int_0^T \int_I \left( \frac{|\gamma|}{\alpha} (\psi^\alpha - \psi) u_{xxt} - \beta(w^\alpha - w) u_{tt} \right) \varphi \, dx \, dt + \int_I \mathcal{E}^\alpha(\mathbf{U}_0^\alpha|\mathbf{U}_0)\varphi(\cdot, 0) \, dx \ge 0. \nonumber
\end{align}
$$

Since we are using periodic boundary conditions for the bounded domain $\mathcal{I}$, we get for the test function
$\varphi(x,t) = \vartheta(t) \in C_0^{\text{Lip}}([0,T))$ the derivatives $\vartheta_t = \vartheta'(t)$ and $\vartheta_x = 0$. In such a situation, the integral identity (3.16) simplifies to the following form

$$ 
\begin{aligned}
& \int_0^T \int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U})\vartheta' dx dt + \int_I (\mathcal{E}^\alpha(\mathbf{U}_0^\alpha|\mathbf{U}_0))\vartheta(\cdot, 0) dx \\
& \geq -\int_0^T \int_I \operatorname{sgn}(\gamma)(f(u^\alpha) - f(u) - f'(u)(u^\alpha - u))(f'(u)u_x + \psi_x)\vartheta dx dt \\
& \quad - \int_0^T \int_I \left(\frac{|\gamma|}{\alpha}(\psi^\alpha - \psi)u_{xxt} - \beta(w^\alpha - w)u_{tt}\right)\vartheta dx dt.
\end{aligned} 
\quad (3.17) $$

Following the proof of Theorem 5.2.1 in [1], we fix the test function $\vartheta$ for some time $\sigma \in (0, T)$. The test function $\vartheta$ for some $\epsilon > 0$ is given by

$$ \vartheta(t) = \begin{cases} 1, & 0 \le t < \sigma, \\ \frac{\sigma - t}{\epsilon} + 1, & \sigma \le t \le \sigma + \epsilon, \\ 0, & t \ge \sigma + \epsilon. \end{cases} \quad (3.18) $$

such that

$$ \vartheta'(t) = \begin{cases} 0, & 0 \le t < \sigma, \\ -\frac{1}{\epsilon}, & \sigma < t < \sigma + \epsilon, \\ 0, & t \ge \sigma + \epsilon. \end{cases} \quad (3.19) $$

Clearly, (3.17) then reduces to

$$ \frac{1}{\epsilon} \int_{\sigma}^{\sigma+\epsilon} \int_I \mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U}) dx dt \le \int_I \mathcal{E}^{\alpha}(\mathbf{U}_0^{\alpha}|\mathbf{U}_0) dx + \underbrace{\int_0^{\sigma} \int_I \mathcal{R}\vartheta(t) dx dt}_{=: \mathcal{I}_4} + \underbrace{\int_{\sigma}^{\sigma+\epsilon} \int_I \mathcal{R}\vartheta(t) dx dt}_{=: \mathcal{I}_5}, \quad (3.20) $$

where $\mathcal{R} = \operatorname{sgn}(\gamma)(f(u^\alpha) - f(u) - f'(u)(u^\alpha - u))(f'(u)u_x + \psi_x) + \left(\frac{|\gamma|}{\alpha}(\psi^\alpha - \psi)u_{xxt} + \beta(w^\alpha - w)u_{tt}\right)$.

Now for $\mathcal{I}_4$, in view of Young's inequality, we have

$$ \begin{aligned} \mathcal{I}_4 \le & \|(f'(u)u_x + \psi_x)\|_{L^\infty(I \times (0, \sigma))} \int_0^\sigma \int_I |(f(u^\alpha) - f(u) - f'(u)(u^\alpha - u))| dx dt \\ & + \int_0^\sigma \int_I \left(\frac{1}{2\alpha}(\psi^\alpha - \psi)^2 + \frac{\beta}{2}(w^\alpha - w)^2\right) dx dt + \int_0^\sigma \int_I \left(\frac{\gamma^2}{2\alpha}u_{xxt}^2 + \frac{\beta}{2}u_{tt}^2\right) dx dt. \end{aligned} \quad (3.21) $$

In view of Assumption A, we have $\|f''\|_\infty \le L$ and thus

$$ (f(u^\alpha) - f(u) - f'(u)(u^\alpha - u)) \le \frac{L}{2}(u^\alpha - u)^2. \quad (3.22) $$

Moreover, the uniform convexity of $\operatorname{sgn}(\gamma)F$ with $\operatorname{sgn}(\gamma)f'(u) \ge \delta > 0$ ensures

$$ \frac{L}{2}(u^\alpha - u)^2 \le \frac{L \operatorname{sgn}(\gamma)}{2\delta} (F(u^\alpha) - F(u) - f(u)(u^\alpha - u)). \quad (3.23) $$

Therefore, using (3.22) and (3.23) in (3.21) and in view of the relative entropy from (3.8), we have

$$ \mathcal{I}_4 \le C_3 (L, \delta) \|u_{xxx}\|_{L^\infty(I \times (0, \sigma))} \int_0^\sigma \int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}) dx dt + \int_0^\sigma \int_I \left( \frac{\gamma^2}{2\alpha} u_{xxt}^2 + \frac{\beta}{2} u_{tt}^2 \right) dx dt, $$

where the constant $C_3 = C(L, \delta)$ does not depend on the parameters $\alpha$ and $\beta$.

Furthermore, for $\mathcal{I}_5$ in (3.20), one can apply the Cauchy-Schwarz inequality to obtain the inequality

$$ \mathcal{I}_5 \le \sqrt{\epsilon|I|} \|\mathcal{R}\|_{L^2(I \times (\sigma, \sigma+\epsilon))} \to 0, \quad \text{as } \epsilon \to 0. $$
Thus, by letting $\epsilon \to 0$ in (3.20), we obtain using the Sobolev inequality

$$ \int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha(x, \sigma)|\mathbf{U}(x, \sigma)) dx \le \int_I \mathcal{E}^\alpha(\mathbf{U}_0^\alpha|\mathbf{U}_0) dx + \int_0^\sigma \int_I \left(\frac{\gamma^2}{2\alpha}u_{xxt}^2 + \frac{\beta}{2}u_{tt}^2\right) dx dt \\ + C_3(L, \delta) \|u\|_{L^\infty((0,\sigma)\times H^4(I))} \int_0^\sigma \int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}) dx dt. $$

Therefore, arguing as in Theorem 5.2.1 in [1] and using Gronwall's inequality, one can conclude that

$$ \sup_{t \in (0,T)} \int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}) dx \\ \le \exp \left( C_3 T \|u\|_{L^\infty(0,T;H^4(I))} \right) \left( \int_I \mathcal{E}^\alpha(\mathbf{U}_0^\alpha|\mathbf{U}_0) dx + \frac{\gamma^2}{2\alpha} \|u_{xxt}\|_{L^2(I_T)}^2 + \frac{\beta}{2} \|u_{tt}\|_{L^2(I_T)}^2 \right). \quad (3.24) $$

From our choice of prepared initial data $\mathbf{U}_0^\alpha = \mathbf{U}_0$ (see (2.8)) we deduce

$$ \int_I \mathcal{E}^\alpha(\mathbf{U}_0^\alpha|\mathbf{U}_0) dx \\ = \operatorname{sgn}(\gamma) (F(u_0^\alpha) - F(u_0) - f(u_0)(u_0^\alpha - u_0)) + \frac{|\gamma|}{2}(p_0^\alpha - p_0)^2 + \frac{\beta}{2}(w_0^\alpha - w_0)^2 + \frac{1}{2\alpha}(\psi_0^\alpha - \psi_0)^2 = 0. $$

Then if $\alpha \to \infty$, we get with $\beta = \alpha^{-1}$ the estimate $\int_I \mathcal{E}^\alpha(\mathbf{U}^\alpha|\mathbf{U}) dx \le C_4\alpha^{-1}$. The convergence rate in (3.4) is a consequence of (3.9). $\square$

**Remark 3.2.** To prove the relative entropy estimates, the condition $f'' \in L^\infty(I)$ can be relaxed. What is really needed is that there exists a constant $K > 0$ such that

$$ |f(u^\alpha) - f(u) - f'(u)(u^\alpha - u)| \le K|u^\alpha - u|^2 \quad \text{uniformly in } I_T. $$

## 3.2 Diffusive-dispersive asymptotics for the hyperbolic-parabolic system (2.26)

In Theorem 2.7, we already proved the global well-posedness of the classical solutions for the system (2.26). In this section, we proceed to analyze the convergence of solutions for the initial value problem (2.26), (2.27) with $\epsilon > 0$ to the solutions of (2.23), (2.24). We use again the scaling $\beta = \alpha^{-1}$, and we employ again relative entropy estimates for our second main result in Theorem 3.3. Precisely, we prove the following.

**Theorem 3.3 (Main result II).** Let $u_0 \in H^6(I)$, and suppose that Assumption A holds. We assume that $u \in L^2(0,T; H^3(I))$ with $u_{tt} \in L^2(I_T)$ and $\|u\|_{L^\infty(0,T;H^4(I))} < \infty$ is a classical solution of (2.23), (2.24). Consider for $\beta = \alpha^{-1}$, a sequence $\{\mathbf{U}^\alpha = (u^\alpha, \psi^\alpha, w^\alpha, p^\alpha)^\top\}_{\alpha>0}$ of smooth solutions of (2.26), (2.27).

Then there is a constant $C = C(L, \delta, u) \ge 0$ that is independent of $\alpha$ such that

$$ \|u^\alpha - u\|_{L^\infty(0,T;L^2(I))} + \|p^\alpha - u_x\|_{L^\infty(0,T;L^2(I))} \le C\sqrt{\frac{1}{\alpha}} \quad (3.25) $$

holds.

*Proof.* The proof follows the lines of the proof for Theorem 3.1 but exploits the fact that $u$ and the elements of the sequence $\{\mathbf{U}^\alpha\}_{\alpha>0}$ are smooth functions.

Define $\psi = -|\gamma|u_{xx}$, $w = u_t$ and $p = u_x$. Then we deduce from the fact that $u$ solves the diffusive-
dispersive equation (2.23) the validity of the system

$$
\begin{align}
u_t + f(u)_x + \operatorname{sgn}(\gamma)\psi_x &= \varepsilon u_{xx}, \nonumber \\
\frac{1}{\alpha}\psi_t + f(u)_x + \operatorname{sgn}(\gamma)\psi_x &= -w + \frac{\varepsilon}{\alpha}\psi_{xx} - \frac{1}{\alpha}|\gamma|u_{xxt}, \tag{3.26} \\
\beta w_t - |\gamma|p_x &= \psi + \beta u_{tt}, \nonumber \\
p_t - w_x &= \varepsilon p_{xx}. \nonumber
\end{align}
$$

Subtracting (3.26) from (2.26) we obtain the error system

$$
\begin{align}
(u^\alpha - u)_t + (f(u^\alpha) - f(u))_x + \operatorname{sgn}(\gamma)(\psi^\alpha - \psi)_x &= \varepsilon(u^\alpha - u)_{xx}, \nonumber \\
\frac{1}{\alpha}(\psi^\alpha - \psi)_t + (f(u^\alpha) - f(u))_x + \operatorname{sgn}(\gamma)(\psi^\alpha - \psi)_x &= -(w^\alpha - w) + \frac{\varepsilon}{\alpha}(\psi^\alpha - \psi)_{xx} + \frac{1}{\alpha}|\gamma|u_{xxt}, \tag{3.27} \\
\beta(w^\alpha - w)_t - |\gamma|(p^\alpha - p)_x &= (\psi^\alpha - \psi) - \beta u_{tt}, \nonumber \\
(p^\alpha - p)_t - (w - w)_x &= \varepsilon(p - p)_{xx}. \nonumber
\end{align}
$$

We recall the definition of the relative entropy from (3.8). A straightforward computation using the lower bound $\delta$ of $\operatorname{sgn}(\gamma)f'$ and the upper bound $L$ of $f''$ from Assumption A yields then

$$
\begin{align*}
\frac{d}{dt} \int_I \mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U}) dx + \varepsilon \int_I \left( \delta(u_x^{\alpha} - u_x)^2 + \frac{1}{\alpha}(\psi_x^{\alpha} - \psi_x)^2 + |\gamma|(p_x^{\alpha} - p_x)^2 \right) dx \\
\le C_5(\delta, L) \|u_t\|_{\infty} \int_I \mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U}) dx + \frac{\gamma^2}{\alpha} \|u_{xxt}\|_{L^2}^2 + \beta \|u_{tt}\|_{L^2}^2.
\end{align*}
$$

The constant $C_5$ is independent of $\alpha$. Then, using Sobolev embedding and Gronwall's inequality, we have

$$
\begin{align}
\sup_{t \in (0,T)} \int_I \mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U}) dx &\le \exp \left( C_5(\delta, L) \left( \|u\|_{L^{\infty}(0,T;H^4(I))} \right) T \right) \nonumber \\
&\quad \times \left( \int_I \mathcal{E}^{\alpha}(\mathbf{U}_0^{\alpha}|\mathbf{U}_0) dx + \frac{\gamma^2}{\alpha} \|u_{xxt}\|_{L^2}^2 + \beta \|u_{tt}\|_{L^2}^2 \right) \tag{3.28}
\end{align}
$$

The estimate (3.25) is then a consequence of the lower bound (3.9) for the relative entropy $\mathcal{E}^{\alpha}(\mathbf{U}^{\alpha}|\mathbf{U})$ and the fact that $\mathcal{E}^{\alpha}(\mathbf{U}_0^{\alpha}|\mathbf{U}_0) = 0$ due to the choice of the prepared initial data (2.27). $\square$

For numerical evidence supporting the statement of Theorem 3.3 under the scaling choice $\beta = \alpha^{-1}$ we refer to Section 4.2.5 below.

# 4 Numerical results

In this section, we provide a series of numerical examples for our new approximate systems (2.7) and (2.26). For these hyperbolic and hyperbolic-parabolic systems, reliable numerical schemes are available. We show that our chosen approach can handle strong nonlinear effects leading to oscillating or shock-wave-like behavior. In particular, we provide test cases for a variety of diffusive-dispersive equations ranging from KdV, Gardner, and modified KdV-Burgers equations. We also provide test cases to capture non-classical shock waves such as dispersive shock waves (DSWs) or undercompressive shock waves with our proposed systems, which is very rare in the literature to the best of our knowledge. The results confirm the convergence results obtained in Section 3.

The relaxation approach we propose here applies to the rather general class of PDEs (1.1). Therefore, the numerical results we present cover several cases, for both convex and non-convex flux functions, in the purely dispersive or diffusive-dispersive cases; see also Remark 2.6(iii).
## 4.1 The numerical scheme

Our discretization approach is described using the abstract form

$$ \mathbf{U}_t^\alpha + \mathbf{f}(\mathbf{U}^\alpha)_x = \mathbf{S}(\mathbf{U}^\alpha)\mathbf{U}^\alpha + \mathbf{D}\mathbf{U}_{xx}^\alpha \quad (4.1) $$

for both approximate systems. For $\gamma \in \mathbb{R} \setminus \{0\}$, we refer to (2.7) for the first-order hyperbolic approximation of the dispersive equation (2.1) with $\varepsilon = 0$ and to (2.26) for the second-order hyperbolic-parabolic approximation of (2.1) with $\varepsilon > 0$.

The computational domain $I_c = [x_L, x_R]$ is subdivided into $N$ equidistant cells of length $\Delta x = (x_R - x_L)/N$. The $i^{th}$ cell is identified by its center $x_i = x_L + (i + 1/2) \Delta x$ and is delimited by the boundaries $x_{i \pm 1/2}$. The time $t \in [0, T]$ is discretized so that $t^n = n \Delta t$, where $\Delta t$ is subject to a CFL-like condition. We refer to the cell average of the numerical solution in cell $i$, at time $t^n$ by $\mathbf{U}_i^n$.

Let us point out that there are by now multiple well established numerical methods for (4.1), see, e.g., [46] for some recent retreat. We propose a semi-implicit scheme based on an explicit second-order two-step Lax-Wendroff method for the hyperbolic first-order operator of the system, combined with an implicit treatment of the stiff source term and a standard central discretization for the second-order derivatives. The choice of a Lax-Wendroff scheme is justified by its low numerical dissipation, especially since the considered system admits potentially fast characteristic speeds. An explicit diffusive scheme may require an excessively refined mesh to conserve a good quality of the numerical solution at large times. An operator splitting approach for the parabolic terms is employed for the parabolic terms for the test cases where they are present. In this setting, at each time step and in every computational cell, the new solution is computed by first evaluating the intercell values $\mathbf{U}_{i \pm 1/2}^{n+1/2}$ via

$$ \mathbf{U}_{i+1/2}^{n+1/2} = \frac{1}{2}(\mathbf{U}_i^n + \mathbf{U}_{i+1}^n) - \frac{\Delta t}{2\Delta x}(\mathbf{f}(\mathbf{U}_{i+1}^n) - \mathbf{f}(\mathbf{U}_i^n)) + \frac{\Delta t}{2}\mathbf{S}\mathbf{U}_{i+1/2}^{n+1/2}, \quad (4.2a) $$

$$ \mathbf{U}_{i-1/2}^{n+1/2} = \frac{1}{2}(\mathbf{U}_i^n + \mathbf{U}_{i-1}^n) - \frac{\Delta t}{2\Delta x}(\mathbf{f}(\mathbf{U}_i^n) - \mathbf{f}(\mathbf{U}_{i-1}^n)) + \frac{\Delta t}{2}\mathbf{S}\mathbf{U}_{i-1/2}^{n+1/2}, \quad (4.2b) $$

whose values can be computed explicitly as

$$ \mathbf{U}_{i \pm 1/2}^{n+1/2} = \left( \mathbf{I} - \frac{\Delta t}{2} \mathbf{S} \right)^{-1} \left( \frac{1}{2} (\mathbf{U}_i^n + \mathbf{U}_{i \pm 1}^n) \mp \frac{\Delta t}{2\Delta x} (\mathbf{f}(\mathbf{U}_{i \pm 1}^n) - \mathbf{f}(\mathbf{U}_i^n)) \right). \quad (4.3) $$

The solution of the hyperbolic first-order part of the system is then given by

$$ \mathbf{U}_i^* = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x} \left( \mathbf{f} \left( \mathbf{U}_{i+1/2}^{n+1/2} \right) - \mathbf{f} \left( \mathbf{U}_{i-1/2}^{n+1/2} \right) \right) + \frac{\Delta t}{2} \mathbf{S} \left( \mathbf{U}_{i+1/2}^{n+1/2} + \mathbf{U}_{i-1/2}^{n+1/2} \right). $$

Finally, the solution at the new time, accounting for the second-order terms (present for $\varepsilon > 0$) is

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^* + \frac{\Delta t}{\Delta x^2} \mathbf{D} (\mathbf{U}_{i+1}^* - 2\mathbf{U}_i^* + \mathbf{U}_{i-1}^*). $$

The time-step is restricted by a CFL-type condition of the form

$$ \Delta t = \text{CFL} \frac{\Delta x}{\lambda + \frac{2\varepsilon}{\Delta x}}, \quad \text{where} \quad \text{CFL} \le 1 \quad \text{and} \quad \lambda = \max_{i \in \{1, \dots, N\}} \max_{k=1, \dots, 4} |\lambda_k(\mathbf{U}_i^n)|, $$

where $\lambda$ is the maximum signal speed and CFL is the Courant-Friedrichs-Levy number, which we take equal to 0.9 in all the carried out simulations.

## 4.2 Numerical results

In an attempt to help the reader navigate this section, we recall at the beginning of each test case which evolution equation we are approximating and for which initial conditions. We shall refer to the latter as the *Original Initial Value Problem* (OIVP). The numerical computations are then carried out by using the previously proposed scheme to solve the *Approximate Initial Value Problem* (AIVP)
consisting of the PDE (2.7) (or (2.26)) and for which the well-prepared initial conditions are supplied as follows

$$u(x, 0) = u_0(x), \quad \psi(x, 0) = -|\gamma|u_{0,xx}, \quad w(x, 0) = -f(u_0)_x + |\gamma|u_{0,xxx}, \quad p(x, 0) = u_{0,x}.$$

Note that $u_0$ is the equivalent initial condition function associated with the original problem (1.1). Concerning boundary conditions, we use either periodic boundary conditions or mixed Neumann-Dirichlet boundary conditions. Although the former is rather straightforward, the latter is designed to mimic homogeneous Neumann boundary conditions $u_x(x_L, t) = u_x(x_R, t) = 0$. We shall refer to this as pseudo-Neumann boundary conditions. Since the unknown $c$ relaxes towards $u$ (see (2.4) and (2.25)), we impose the same boundary conditions for both. In terms of the evolved variables $u, \psi, w, p$ this translates to

$$u_x(x_L, t) = u_x(x_R, t) = 0, \quad \psi_x(x_L, t) = \psi_x(x_R, t) = 0, \qquad (4.4a)$$

$$w_x(x_L, t) = w_x(x_R, t) = 0, \quad p(x_L, t) = p(x_R, t) = 0. \qquad (4.4b)$$

### 4.2.1 One-Soliton solution (KdV equation)

In this example, the OIVP is given by

$$\left\{ \begin{array}{l} u_t + 6uu_x = \gamma u_{xxx}, \\ u(x, 0) = u_0(x) = \frac{V}{2} \mathrm{sech}^2 \left( \sqrt{\frac{V}{4|\gamma|}} x \right). \end{array} \right. \qquad (4.5)$$

The exact solution to (4.5) is a traveling wave, i.e., $u(x, t) = u_0(x - Vt)$. The computational domain is $I_c = [-2, 2]$ with $N = 1000$ cells. We take $\gamma = -10^{-2}$. The values of the relaxation parameters are $\alpha = 10^3$ and $\beta = 10^{-6}$. We take $V = 1$ and the final time is $T = 100$. Here, we use periodic boundary conditions so that the final time corresponds to 25 periods. The periodization of the solution naturally introduces negligible errors at the boundaries since the solution is exponentially decaying but not perfectly flat. The numerical results are provided in Figure 1. The results show that the soliton

<figure><img src="image_4.png" /><figcaption>Figure 1: Numerical results for the AIVP corresponding to (4.5). Left: The numerical solution $u(x)$ at different times up to the final time $T = 100$. Right: The time evolution of the amplitude error $e_a$ as well as the absolute $L^2$ error between the numerical solution and the exact solution to the OIVP (4.5). $I_c = [-2, 2]$ with $N = 1000$ cells, $\gamma = -10^{-2}$, $\alpha = 10^3$ and $\beta = 10^{-6}$.</figcaption></figure>

remains stable for long times and the soliton does not break, not even after 25 periods. The numerical solution seemingly behaves as a traveling wave, with the same speed and profile as the solution to
(4.5). The long-time behavior of the numerical solution is quantitatively illustrated in the right part of Figure 1, where the amplitude error $e_a(t)$ and the absolute $L^2$-error between the numerical solution of the AIVP and the exact solution of the OIVP $e_{L^2}(t)$, given by

$$e_a(t^n) = \left| \max_i (u_i^n) - V/2 \right|, \quad e_{L^2}(t^n) = \sqrt{\sum_i (u_i^n - u_0(x_i - Vt^n))^2},$$

are plotted over time. The evolution of the error over the time interval $[0, 100]$ shows no significant increase in the order of magnitude.

This shows the validity of our numerical approach even for large times. We conjecture that the long-time dynamics of the solution of the approximate system (2.7) is also governed by a stable soliton solution under small perturbations, due to both periodization and relaxation. Thus, the approximation resembles qualitative properties of the original equation (2.1).

### 4.2.2 Two-Soliton solution (KdV equation)

To address a more complex example, we consider an OIVP given by

$$\left\{ \begin{array}{l} u_t + 6uu_x = -u_{xxx}, \\ u(x, 0) = \frac{2(V_1 - V_2)(V_1 \cosh^2 \theta_2(x) + V_2 \sinh^2 \theta_1(x))}{((\sqrt{V_1} + \sqrt{V_2}) \cosh(\theta_1(x) - \theta_2(x)) + (\sqrt{V_1} - \sqrt{V_2}) \cosh(\theta_1(x) + \theta_2(x)))^2}, \end{array} \right. \qquad (4.6)$$

where $\theta_j(x) = \sqrt{V_j}(x - x_j)/2$, $j = 1, 2$. The computational domain is $I_c = [-15, 15]$ with $N = 1000$ cells. The values of the relaxation parameters are $\alpha = 4 \times 10^3$ and $\beta = 10^{-6}$. Periodic boundary conditions are used and we set the final time to $T = 60$. The solution consists of two solitons that move with different velocities and pass each other. Upon interacting, a phase shift occurs before each soliton resumes its course. Different aspects of the numerical solution are provided in Figure 2.

Again, we observe that our numerical approach provides valid results for large times. We point out that the approximate system (2.7) allows for crossings of solitons just like the dispersive target problem (2.1). In Figure 3, we show the evolution over time of the energy error with respect to its initial value, that is $E^\alpha - E_0^\alpha$. We deliberately do not use absolute value, only to emphasize that the total energy computed from the numerical solution, increases from its initial value and does not show monotonicity, given that the used numerical scheme is not entropy stable. Nevertheless, we report in Table 1 the values of $||E^\alpha - E_0^\alpha||_{L^\infty(0,T)}$, corresponding to the amplitude of the periodic peaks, and which show second order convergence rate.

<table><caption>Table 1: Convergence of $||E^\alpha - E_0^\alpha||_{L^\infty(0,T)}$, corresponding to the curves shown in Figure 3.</caption><thead><tr><td>N</td><td>||E<sup>α</sup> - E<sub>0</sub><sup>α</sup>||<sub>L<sup>∞</sup>(0,T)</sub></td><td>Order</td></tr></thead><tbody><tr><td>500</td><td>6.32 × 10<sup>-3</sup></td><td>—</td></tr><tr><td>1000</td><td>1.46 × 10<sup>-3</sup></td><td>2.11</td></tr><tr><td>2000</td><td>3.50 × 10<sup>-4</sup></td><td>2.06</td></tr><tr><td>4000</td><td>8.57 × 10<sup>-5</sup></td><td>2.02</td></tr><tr><td>8000</td><td>2.10 × 10<sup>-5</sup></td><td>2.02</td></tr></tbody></table>

### 4.2.3 Negative hump DSW (KdV equation)

The OIVP is given by

$$\left\{ \begin{array}{l} u_t + 6uu_x = \gamma u_{xxx}, \\ u(x, 0) = -\operatorname{sech}^2(x). \end{array} \right. \qquad (4.7)$$
<figure><img src="image_5.png" /><figcaption>Figure 2: Numerical results for the AIVP corresponding to (4.6), showing the main dynamics of the two-soliton solution. The parameters are $I_c = [-15, 15]$ with $N = 1000$ cells. $\gamma = -1$, $\alpha = 4 \times 10^3$, $\beta = 10^{-6}$. $T = 60$ and boundary conditions are periodic.</figcaption></figure>

This test case is adapted from [47], where the initial value problem (4.7) is analyzed. The computational domain is $I_c = [-5, 5]$ with $N = 8000$ cells. We take $\gamma = -10^{-4}$. The values of the approximation parameters are $\alpha = 10^3$ and $\beta = 10^{-7}$. The final time is $T = 0.4$. The solution dynamics are such that the smooth initial datum evolves up to a given finite time, where it develops a wave-train of rapid oscillations in a region of space expanding over time. This is usually referred to as a Collisionless [48] or Dispersive Shock Wave (DSW). The main dynamics of the numerical solution are shown in Figure 4. In particular, in Figure 4a, we show the numerical solution at different times $t \in \{0, 0.1, t_c, 0.4\}$ where $t_c = 0.216$ is the predicted time marking the onset of the DSW, see [47]. In Figure 4b, we compare $p^\alpha$ and $u_x^\alpha$ in the DSW region at the final time $t = 0.4$. Then, we also compare the numerical solution at the latter time with the corresponding Whitham envelope for the solution to the OIVP (4.7) in Figure 4c, which shows excellent agreement. Finally, the space-time evolution of the DSW region is summarized in the $x-t$ diagram of Figure 4d, superimposed with the asymptotic $(x, t)$ coordinates of the DSW edges. All these results show an excellent performance of the method, since the resulting numerical solution to the AIVP is seemingly indistinguishable from the exact solution of the OIVP (4.7), in terms of space-time evolution but also regarding the dispersion effects in a test case where they are not negligible.
<figure><img src="image_6.png" /><figcaption>Figure 3: Time evolution of the total energy error with respect to its initial value $E^{\alpha} - E_{0}^{\alpha}$ for different mesh sizes. The parameters are $I_{c} = [-15, 15]$, $\gamma = -1$, $\alpha = 2 \times 10^{3}$, $\beta = 10^{-8}$, $T = 100$, and boundary conditions are periodic.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(a) Numerical solution $u$ at $t \in \{0, 0.1, 0.216, 0.4\}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) $p^{\alpha}$ and $u_{x}$ in the DSW region at $t = 0.4$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Comparison of the numerical solution with the Whitham envelope for the OIVP (4.7).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) $x-t$ diagram of the numerical solution at the onset of the DSW formation.</figcaption></figure>

<caption>Figure 4: Numerical results for the AIVP corresponding to (4.7). The parameters are $\gamma = -10^{-4}$, $\alpha = 10^{3}$, $\beta = 10^{-7}$. The computational domain is $I_{c} = [-5, 5]$ with $N = 8000$ cells, and the final time is $T = 0.4$. Periodic boundary conditions are used.</caption>
### 4.2.4 Riemann Problem DSW (KdV equation)

The OIVP is given by

$$ egin{cases} u_t + uu_x = \gamma u_{xxx}, \\ u(x, 0) = \begin{cases} 1 & \text{if } x < 0, \\ 0 & \text{if } x \ge 0. \end{cases} \end{cases} \qquad (4.8) $$

The computational domain is $I_c = [-8, 2]$ with $N = 10000$ cells. We take $\gamma = -10^{-4}$. The values of the relaxation parameters are $\alpha = 10^3$ and $\beta = 10^{-7}$. The final time is $T = 3$. With $\gamma \neq 0$, the initial discontinuity develops a DSW (or undular bore). The derivative of the initial data in the sense of distributions is a minus Dirac delta function, which does not belong to $L^2(I)$. The total energy (2.14) would be formally infinite should $p$ be initialized as $p = u_{0,x}$ for the AIVP. Instead, we consider an analytic approximation of the initial step function given by

$$ \tilde{u}_0(x) = \frac{1}{2} (1 - \tanh(x/\omega)), \quad \omega = 10^{-3}. $$

The asymptotic solution to the OIVP (for $t \to +\infty$) was studied in [48] and is given by

$$ u(x, t) = 2 \operatorname{dn}^2 \left( \frac{1}{\sqrt{6}} \left( x - \frac{1+s^2}{3} t \right), s \right) - (1-s^2), \qquad (4.9) $$

where $dn$ is the Jacobi delta amplitude function and $s$ is the elliptic modulus, and which obeys the implicit equation

$$ \frac{1 + s^2}{3} - \frac{2}{3} \frac{s^2(1 - s^2)K(s)}{E(s) - (1 - s^2)K(s)} = \frac{x}{t}, \qquad (4.10) $$

with K and E, the complete elliptic integrals of first and second kind, respectively. The asymptotic bounds of the DSW region are then determined by taking the limits $s \to 0$ and $s \to 1$ in (4.10) and which yields [48]

$$ \frac{x}{t} = \tau^{-} = -1, \quad \frac{x}{t} = \tau^{+} = \frac{2}{3}. $$

Upon computing $s$ as a function of $x/t$ in the DSW region, the local extrema of (4.9) define the upper and lower asymptotic envelopes $A_{\pm}(x/t)$ defined as

$$ A_{\pm}(x/t) = 1 \pm s^2(x/t). \qquad (4.11) $$

The numerical results are given in Figure 5. The comparison shows that both the amplitude of the oscillations as well as their position and velocity match well with the asymptotic values of the original equation, confirming once again the quality of the approximation. Given the structure of the solution, the vanishing oscillations propagating from the trailing edge of the DSW to the leftmost boundary are usually a source of concern from the numerical point of view especially when the boundary conditions are not completely transparent [49]. Nevertheless, the numerical solution seems to perform well in that regard, even by simply implementing pseudo-Neumann boundary conditions. To further highlight the relevancy of this choice, we compare this same test case using pseudo-Neumann boundary conditions against classic Neumann boundary conditions, for which $p_x(x_L, t) = p_x(x_R, t) = 0$. The comparison is reported on Figure 6 and shows that classic Neumann boundary conditions fail to maintain the correct structure of the solution near the boundary. The consequences of this erroneous choice are visibly more significant on the variable $u^\alpha$ than on $p^\alpha$ itself.
<figure><img src="image_11.png" /><figcaption>Figure 5: Numerical results for the AIVP corresponding to (4.7). Left: Numerical solution at $t = 2$ compared with the asymptotic envelopes $A_±$. Right: $x-t$ diagram of the numerical solution, superimposed with the asymptotic DSW bounds. The parameters are $γ = -10^{-4}$, $α = 10^3$, $β = 10^{-7}$. The computational domain is $I_c = [-8, 2]$ with $N = 10000$ cells, and the simulation is run up to $T = 3$. Pseudo Neumann boundary conditions are used.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 6: Comparison of pseudo Neumann boundary conditions (4.4) with classic Neumann boundary conditions for which $p_x(x_L, t) = p_x(x_R, t) = 0$, for the DSW test case. The computational domain is $I_c = [-6, 2]$ with $N = 10000$ cells. The time is $T = 1.5$. The rest of the parameters are $α = 10^3$, $β = 10^{-8}$, $γ = -10^{-4}$.</figcaption></figure>

### 4.2.5 Traveling wave front of KdV-Burgers

The OIVP is given by

$$ \begin{cases} u_t + uu_x = \varepsilon u_{xx} + \gamma u_{xxx}, \\ u(x, 0) = u_0(x) = \frac{-3\varepsilon^2}{25\gamma} \left( 2 + 2 \tanh \left( \frac{\varepsilon x}{10\gamma} \right) + \operatorname{sech}^2 \left( \frac{\varepsilon x}{10\gamma} \right) \right). \end{cases} \quad (4.12) $$

The exact solution to this IVP is a traveling wave $u_0(x - Vt)$ where $V = -\frac{6\varepsilon^2}{25\gamma}$. We take for this test case $I_c = [-3, 1]$, discretized over $N = 2000$ computational cells. We take $γ = 10^{-4}$, $ε = 10^{-2}$ and the values of the relaxation parameters are $α = 2 \times 10^3$ and $β = 10^{-6}$. The final time is $T = 10$. Pseudo Neumann boundary conditions (4.4) are used. The numerical results are given in Figure 7.
In particular, we show the numerical solution $u^\alpha$ at different times as well as the time evolution of $E^\alpha$. To validate the decay rate of the solution, we compared it with the reference decay law (2.29), integrated over the computational domain and in time using an explicit Euler scheme. In Figure 8,

<figure><img src="image_13.png" /><figcaption>Figure 7: Left: The numerical solution to the AIVP corresponding to (4.12) at different times. Right: the evolution of the energy computed from the numerical solution, compared with the reference decay law, computed by solving the ODE in time. The parameters are $N = 2000$, $\gamma = 10^{-4}$, $\varepsilon = 10^{-2}$, $\alpha = 2 \times 10^3$, $\beta = 10^{-6}$ and the final time is $T = 10$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 8: Left: $L^2$-error as a function of time at fixed $\alpha = 10^4$ with variable $\beta$. Right: $L^2$-error as a function of time at fixed $\beta$ with variable $\alpha$. Other parameters are $N = 20000$, $\gamma = 10^{-4}$, $\varepsilon = 10^{-2}$ and $T = 0.1$.</figcaption></figure>

we display errors for different pairings of the approximation parameters $\alpha$ and $\beta$. The error grows linearly in time as expected. For fixed $\alpha$ or $\beta$ we observe a clear error decay if $\beta$ tends to zero or $\alpha$ increases, respectively.

The plots in Figure 9 support our analytical findings from Theorem 3.25. As expected, to decrease the error, both parameters $\alpha$ and $\beta$ have to be tuned. The linear error decay that can be observed from the right plot in Figure 9 is computed using the scaling $\beta = \gamma\alpha^{-1}$, which corresponds to the scaling $\beta \sim \alpha^{-1}$ from Theorem 3.25. In Table 2, we display corresponding numerical quantities from which we compute the experimental order of convergence $O_L^2$.
<figure><img src="image_15.png" /><figcaption>Figure 9: Convergence plots in different error norms. The choice of parameters is as in Figure 8 above. Left: Error decay and saturation for fixed $\beta$ and $\alpha \to \infty$. Middle: Error decay and saturation for fixed $\beta$ and $\alpha \to \infty$. Right: Linear error decay for $\beta = \alpha^{-1}$ and $\alpha \to \infty$.</figcaption></figure>

<table><caption>Table 2: Experimental order of convergence for the setting as in Figure 9. The right column clearly displays a first-order convergence whereas we could prove in Theorem 3.3 only square-root behaviour.</caption><thead><tr><th colspan="3">Fixed β, variable α</th><th colspan="3">Fixed α, variable β</th><th colspan="3">Scaled values β ~ α<sup>-1</sup></th></tr><tr><th>α</th><th>e<sub>L<sup>2</sup></sub></th><th>O<sub>L<sup>2</sup></sub></th><th>β</th><th>e<sub>L<sup>2</sup></sub></th><th>O<sub>L<sup>2</sup></sub></th><th>α</th><th>e<sub>L<sup>2</sup></sub></th><th>O<sub>L<sup>2</sup></sub></th></tr></thead><tbody><tr><td>10<sup>0</sup></td><td>5.29 × 10<sup>-5</sup></td><td>—</td><td>10<sup>-5</sup></td><td>1.27 × 10<sup>-5</sup></td><td>—</td><td>10<sup>0</sup></td><td>1.36 × 10<sup>-4</sup></td><td>—</td></tr><tr><td>10<sup>1</sup></td><td>6.03 × 10<sup>-6</sup></td><td>0.94</td><td>10<sup>-6</sup></td><td>1.49 × 10<sup>-6</sup></td><td>0.93</td><td>10<sup>1</sup></td><td>1.47 × 10<sup>-5</sup></td><td>0.97</td></tr><tr><td>10<sup>2</sup></td><td>6.42 × 10<sup>-7</sup></td><td>0.97</td><td>10<sup>-7</sup></td><td>6.42 × 10<sup>-7</sup></td><td>0.37</td><td>10<sup>3</sup></td><td>1.49 × 10<sup>-6</sup></td><td>1.00</td></tr><tr><td>10<sup>4</sup></td><td>1.31 × 10<sup>-7</sup></td><td>0.69</td><td>10<sup>-8</sup></td><td>6.09 × 10<sup>-7</sup></td><td>0.02</td><td>10<sup>3</sup></td><td>1.31 × 10<sup>-7</sup></td><td>1.06</td></tr><tr><td>10<sup>4</sup></td><td>1.08 × 10<sup>-7</sup></td><td>0.08</td><td>10<sup>-9</sup></td><td>5.98 × 10<sup>-7</sup></td><td>0.01</td><td>10<sup>4</sup></td><td>1.82 × 10<sup>-8</sup></td><td>0.86</td></tr></tbody></table>

### 4.2.6 Undercompressive shock wave(modified KdV-Burgers)

The next OIVP is given for $\gamma > 0$ by

$$ \left\{ \begin{array}{l} u_t + u^2 u_x = \varepsilon u_{xx} - \gamma u_{xxx}, \\ u(x, 0) = u_0(x) = \frac{1}{2} (u_+ + u_- - |u_+ - u_-| \tanh(Ax/\sqrt{\gamma})), \end{array} \right. \qquad (4.13) $$

where

$$ u_{-} = -\sqrt{\frac{2\varepsilon^2}{\gamma}}, \quad u_{+} = -u_{-} - \sqrt{\frac{2\varepsilon^2}{9\gamma}}, \quad A = \frac{u_{-} - u_{+}}{2\sqrt{2}}. $$

For the given values of the solution parameters, this corresponds to a traveling wave solution $u(x, t) = u_0(x - Vt)$, of velocity $V = (u_+^3 - u_-^3)/(u_+ - u_-)$ [6] and consists in a transition from $u_-$ to $u_+$, occurring over a region of space of a length scaling with $\sqrt{\gamma}$. We consider for this test case a sharp transition by setting $\gamma = 10^{-5}$. The computational domain is $I_c = [-0.1, 2]$, discretized over $N = 10000$ cells. We set $\varepsilon = 10^{-2}$. Note that for these values of the parameters, we have approximately $u_- \simeq -4.47$, $u_+ \simeq 3.97$, $A \simeq -2.99$ and $V \simeq 18.02$. The numerical results are gathered in Figure 10. Besides, we show here the influence of the relaxation parameters on the solution, and we consider values of $\alpha$ ranging from 100 to 2500 and values of $\beta$ from $10^{-6}$ down to $10^{-11}$ and for which the results are shown on Figure 11.

We clearly observe that the first-order approximate system is able to display the undercompressive
shock wave in the underlying cubic conservation law. The numerical computation of undercompressive shock waves is difficult because its precise form depends sensitively on the ratio between diffusive and dispersive parameters, see [6]. These can hardly be controlled in numerical schemes unless the shock layer is resolved by the mesh [50]. Our new hyperbolic-parabolic approximation (2.26) allows for stable computations even for underresolved wave structures.

<figure><img src="image_16.png" /><figcaption>Figure 10: Left: Numerical solution for $u^\alpha$ at times $t_0 = 0$, $t_1 = 0.05$, and $t_2 = 0.1$, illustrating the evolution of the undercompressive shock. $x_{1,2} = Vt_{1,2}$ are the theoretical positions of the shock at $t = t_{1,2}$. Right: Comparison between $p^\alpha$ and the finite-difference approximation $u_x^\alpha$ at $t_2 = 0.1$ for several values of $\alpha$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 11: Left: Zoom at the trailing edge of the wave. For $\alpha \to \infty$ we observe convergence to the limit traveling wave. For smaller values of $\alpha$, another rarefaction wave gets visible. Note that the approximate system (2.26) may trigger up to four elementary waves. Right: for fixed value of $\alpha$ and $\beta \to 0$, we observe oscillations and spurious wave speeds for larger values of $\beta$ which are eliminated with $\beta$ vanishing.</figcaption></figure>

The results show that different values of the relaxation parameters exhibit different solution structures, but ultimately show convergence towards the solution of the AIVP. Given that the approximate system is a (4 × 4)-system, the appearance of a rarefaction wave of vanishing amplitude as shown in Figure 11 (Right) in the limit $\alpha \to \infty$ is not surprising. On Figure 11 (Left), the numerical solution shows more drastic changes in its structure as well as the development of spurious oscillations of non-
negligible amplitude for $\beta$ far from the stiff limit. This confirms the fact that both $\alpha$ and $\beta$ have to be scaled carefully in order to correctly capture the original solution's behavior.

### 4.2.7 Dark and bright solitons (Gardner equation)

As final example, we consider as OIVP a Gardner equation given by

$$ \begin{cases} u_t + 6(u - ku^2)u_x + u_{xxx} = 0, \\ u(x, 0) = u_{d,b}(x), \end{cases} \qquad (4.14) $$

where the initial conditions $u_d(x)$ and $u_b(x)$, corresponding to dark or bright solitons, respectively, write as [5]

$$ \begin{aligned} u_d(x) &= u_3 - \frac{u_3 - u_2}{\cosh^2(\theta(x)) - \frac{u_3 - u_2}{u_3 - u_1} \sinh^2(\theta(x))}, \\ u_b(x) &= u_2 + \frac{u_3 - u_2}{\cosh^2(\theta(x)) - \frac{u_3 - u_2}{u_4 - u_2} \sinh^2(\theta(x))}. \end{aligned} $$

Here $u_1 \le u_2 \le u_3 \le u_4$ and $\theta(x) = \frac{1}{2}\sqrt{k(u_3 - u_1)(u_4 - u_2)x}$. The exact solution to (4.14) is then given by $u_{b,d}(x - Vt)$ where the velocity $V$ is given as

$$ V = k(u_1u_2 + u_1u_3 + u_1u_4 + u_2u_3 + u_2u_4 + u_3u_4). \qquad (4.15) $$

These solutions are obtained from a broader class of solutions to the Gardner equations in the limits $u_3 \to u_4$ for $u_d$ and $u_2 \to u_1$ for $u_b$, see [5]. For the numerical simulation of the associated AIVP, we consider a small parameter $\epsilon > 0$, and we take for the dark soliton solution

$$ k = 1, \quad u_1 = 0, \quad u_2 = \epsilon, \quad u_3 = u_4 = 1 - \epsilon, $$

and analogously for the bright soliton

$$ k = 1, \quad u_1 = u_2 = \epsilon, \quad u_3 = 1 - \epsilon, \quad u_4 = 1. $$

This results in table-top-like solitons for $\epsilon \ll 1$. For both parameterizations, we take $\epsilon = 10^{-4}$ and the value of the traveling wave velocity is $V = 1 - \epsilon^2 \simeq 1$. The rest of the parameters for both simulations are as follows. The computational domain $I_c = [-50, 50]$ is discretized over $N = 2000$ cells. We then choose $\alpha = 10^3$ and $\beta = 10^{-6}$. Periodic boundary conditions are used, and we set the final time to $T = 500/V$, corresponding to 10 periods of the solution over the computational domain before returning to its initial state. The numerical results are reported in Figure 12 and underline the validity of our overall approach to approximate diffusive-dispersive equations by extended hyperbolic-parabolic systems.

## 5 Conclusions and future outlook

We have introduced a novel framework for the hyperbolic and hyperbolic-parabolic approximation for a broad class of nonlinear dispersive and diffusive-dispersive equations. The hyperbolic approximation (2.7) is constructed such that a complete set of Riemann invariants can be determined. It is shown to have a Hamiltonian structure that is consistent with the limit dispersive equations. Enforcing the monotonicity on the flux derivative, we have found that the hyperbolic approximate system (2.7) is equipped with an entropy/entropy flux pair that resembles the energy of the limit dispersive equations. Moreover, solutions of the hyperbolic–parabolic approximation (2.26) dissipate this entropy. Based on these findings, we have employed the relative entropy framework to prove the convergence of solutions
<figure><img src="image_18.png" /><figcaption>Figure 12: Left: Numerical solutions to the AIVP corresponding to (4.14) for both bright and dark solitons. Right: Evolution of the total energy error $E^\alpha - E_0^\alpha$ over time. $I_c = [-50, 50]$, $N = 2000$, $\alpha = 10^3$, $\beta = 10^{-6}$ and $T = 500/V$.</figcaption></figure>

of the approximate systems towards smooth solutions of the dispersive and diffusive-dispersive limit equations.

Due to the likewise standard hyperbolic and hyperbolic structure of the approximate systems, it becomes possible to use high-resolution methods from the realm of hyperbolic balance laws for the sake of numerical simulation. For a wide range of diffusive-dispersive equations, we display the effectiveness, accuracy, and robustness of a semi-implicit second-order two-step Lax-Wendroff method. In particular, our focus in the numerical examples was to capture the dynamics of non-classical shock waves, including dispersive and undercompressive shock waves, which are very rare to find in the literature. The numerics suggest convergence of our approximate systems for a significantly larger class of fluxes, independent of the nonlinearity of the flux. The study refers to dispersive equations with nonconvex energy, such as the KdV equation.

As mentioned, the proposed first-order system (2.7) possesses a full set of Riemann invariants. This structure opens the door to possible global existence results for $L^\infty$-initial data via compensated compactness, as well as the development of numerical schemes based on exact Riemann solvers. It would be interesting to analyze whether one can obtain the entropy solution of the hyperbolic system (2.7) from the vanishing diffusion (viscosity) limit of (2.26), i.e., $\varepsilon \to 0$. In this work, we have chosen the scaling between diffusion and dispersion to be 1. For the alternative scaling $\varepsilon^2$, one can also study the diffusion-dispersion limit of the proposed systems as in [21, 22]. We leave these issues for future investigation.

Furthermore, this approach of relaxing higher-order PDEs via energy-consistent lower-order approximations can also be applied to more complex higher-order PDEs, having, for example, nonlinearities in higher-order terms [24,51]. Moreover, this relaxation approach can be readily adapted to coupled multidimensional systems such as Navier-Stokes-Korteweg-type systems, where the energy structure of the underlying mechanisms can be crucial in designing accurate and stable approximate systems [23,52,53].

**Acknowledgements** Funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity, Scales, Randomness (CoScaRa) is gratefully acknowledged.

F.D. acknowledges the support from NextGenerationEU, Azione 247 MUR Young Researchers – SoE line.
**CRediT authorship contribution statement**

**Rahul Barthwal:** Modeling, mathematical analysis, investigation, methodology, writing— original draft, review & editing, **Firas Dhaouadi:** modeling, methodology, numerical simulations, investigation, validation, writing—original draft, review & editing, **Christian Rohde:** investigation, writing—original draft, review & editing, supervision.

**References**

[1] Dafermos CM. 2016 *Hyperbolic conservation laws in continuum physics* vol. 325 *Grundlehren der Mathematischen Wissenschaften*. Berlin: Springer 4th edition.

[2] Whitham GB. 1965 Non-linear dispersive waves. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **283**, 238–261.

[3] Shearer M, Spayd KR, Swanson ER. 2015 Traveling waves for conservation laws with cubic nonlinearity and BBM type dispersion. *Journal of Differential Equations* **259**, 3216–3232.

[4] Bona J, Hong Y. 2022 Numerical study of the generalized Korteweg–de Vries equations with oscillating nonlinearities and boundary conditions. *Water Waves* **4**, 109–137.

[5] Kamchatnov A, Kuo YH, Lin TC, Horng TL, Gou SC, Clift R, El G, Grimshaw RH. 2012 Undular bore theory for the Gardner equation. *Physical Review E—Statistical, Nonlinear, and Soft Matter Physics* **86**, 036605.

[6] El G, Hoefer M, Shearer M. 2017 Dispersive and diffusive-dispersive shock waves for nonconvex conservation laws. *SIAM Review* **59**, 3–61.

[7] Whitham GB. 2011 *Linear and nonlinear waves*. John Wiley & Sons.

[8] Hayes BT, LeFloch PG. 1997 Non-classical shocks and kinetic relations: Scalar conservation laws. *Archive for Rational Mechanics and Analysis* **139**, 1–56.

[9] Hayes B, Shearer M. 1999 Undercompressive shocks and Riemann problems for scalar conservation laws with non-convex fluxes. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics* **129**, 733–754.

[10] LeFloch PG, Shearer M. 2004 Non-classical Riemann solvers with nucleation. *Proceedings of the Royal Society of Edinburgh Section A: Mathematics* **134**, 961–984.

[11] Jacobs D, McKinney B, Shearer M. 1995 Traveling wave solutions of the modified Korteweg-deVries-Burgers equation. *Journal of Differential Equations* **116**, 448–467.

[12] Schulze MR, Shearer M. 1999 Undercompressive shocks for a system of hyperbolic conservation laws with cubic nonlinearity. *Journal of Mathematical Analysis and Applications* **229**, 344–362.

[13] Bertozzi AL, Münch A, Fanton X, Cazabat AM. 1998 Contact line stability and “undercompressive shocks” in driven thin film flow. *Physical Review Letters* **81**, 5169.

[14] Chalons C, Engel P, Rohde C. 2014 A conservative and convergent scheme for undercompressive shock waves. *SIAM Journal on Numerical Analysis* **52**, 554–579.

[15] El G, Hoefer M. 2016 Dispersive shock waves and modulation theory. *Physica D: Nonlinear Phenomena* **333**, 11–65.
[16] Kenig CE, Ponce G, Vega L. 1991 Well-posedness of the initial value problem for the Korteweg-de Vries equation. *Journal of the American Mathematical Society* **4**, 323–347.

[17] Tao T. 2006 *Nonlinear dispersive equations: local and global analysis*. Number 106. American Mathematical Society.

[18] Killip R, Visan M. 2019 KdV is well-posed in $H^{-1}$. *Annals of Mathematics* **190**, 249–305.

[19] Colliander JE, Keel M, Staffilani G, Takaoka H, Tao T. 2003 Sharp global well-posedness for KdV and modified KdV on $\mathbb{R}$ and $\mathbb{T}$. *Journal of the American Mathematical Society* **16**, 705–749.

[20] Farah LG, Linares F, Pastor A. 2011 The supercritical generalized KdV equation: global well-posedness in the energy space and below. *Mathematical Research Letters* **18**, 357–377.

[21] Corli A, Rohde C. 2012 Singular limits for a parabolic–elliptic regularization of scalar conservation laws. *Journal of Differential Equations* **253**, 1399–1421.

[22] Corli A, Rohde C, Schleper V. 2014 Parabolic approximations of diffusive–dispersive equations. *Journal of Mathematical Analysis and Applications* **414**, 773–798.

[23] Dhaouadi F, Dumbser M. 2022 A first order hyperbolic reformulation of the Navier-Stokes-Korteweg system based on the GPR model and an augmented Lagrangian approach. *Journal of Computational Physics* **470**, 111544.

[24] Dhaouadi F, Dumbser M, Gavrilyuk S. 2025 A first-order hyperbolic reformulation of the Cahn-Hilliard equation. *Proceedings A* **481**, 20240606.

[25] Biswas A, Ketcheson DI, Ranocha H, Schütz J. 2025 Traveling-wave solutions and structure-preserving numerical methods for a hyperbolic approximation of the Korteweg-de Vries equation. *Journal of Scientific Computing* **103**, 90.

[26] Favrie N, Gavrilyuk S. 2017 A rapid numerical method for solving Serre–Green–Naghdi equations describing long free surface gravity waves. *Nonlinearity* **30**, 2718.

[27] Gavrilyuk S, Shyue KM. 2022 Hyperbolic approximation of the BBM equation. *Nonlinearity* **35**, 1447.

[28] Gavrilyuk S, Nkonga B, Shyue KM. 2024 The conduit equation: hyperbolic approximation and generalized Riemann problem. *Journal of Computational Physics* **514**, 113232.

[29] Ranocha H, Ketcheson DI. 2025 High-order mass-and energy-conserving methods for the nonlinear Schrödinger equation and its hyperbolization. *arXiv preprint arXiv:2510.14335*.

[30] Duchêne V. 2019 Rigorous justification of the Favrie–Gavrilyuk approximation to the Serre–Green–Naghdi model. *Nonlinearity* **32**, 3772.

[31] Giesselmann J, Ranocha H. 2025 Convergence of hyperbolic approximations to higher-order PDEs for smooth solutions. *arXiv preprint arXiv:2508.04112*.

[32] Giesselmann J, Lattanzio C, Tzavaras AE. 2017 Relative energy for the Korteweg theory and related Hamiltonian flows in gas dynamics. *Archive for Rational Mechanics and Analysis* **223**, 1427–1484.

[33] DiPerna RJ. 1979 Uniqueness of solutions to hyperbolic conservation laws. *Indiana University Mathematics Journal* **28**, 137–188.
[34] Mietka C. 2017 On the well-posedness of a quasi-linear Korteweg-de Vries equation. *Annales Mathématiques Blaise Pascal* **24**, 83–114.

[35] Staffilani G. 1997 On solutions for periodic generalized KdV equations. *IMRN: International Mathematics Research Notices* **1997**.

[36] Guan H, Kuksin S. 2014 The KdV equation under periodic boundary conditions and its perturbations. *Nonlinearity* **27**, R61.

[37] Dhaouadi F, Favrie N, Gavrilyuk S. 2019 Extended Lagrangian approach for the defocusing nonlinear Schrödinger equation. *Studies in Applied Mathematics* **142**, 336–358.

[38] Benjamin TB, Bona JL, Mahony JJ. 1972 Model equations for long waves in nonlinear dispersive systems. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences* **272**, 47–78.

[39] Yong WA. 2004 Entropy and global existence for hyperbolic balance laws. *Archive for Rational Mechanics and Analysis* **172**, 247–266.

[40] Gallego FA, Pazoto AF. 2019 On the well-posedness and asymptotic behaviour of the generalized Korteweg-de Vries-Burgers equation. *Proceedings of the Royal Society of Edinburgh. Section A. Mathematics* **149**, 219–260.

[41] Ikehata R, Takeda H. 2017 Critical exponent for nonlinear wave equations with frictional and viscoelastic damping terms. *Nonlinear Analysis: Theory, Methods & Applications* **148**, 228–253.

[42] Kato T. 1975 The Cauchy problem for quasi-linear symmetric hyperbolic systems. *Archive for Rational Mechanics and Analysis* **58**, 181–205.

[43] Kawashima S. 1984 Systems of a hyperbolic-parabolic composite type, with applications to the equations of magnetohydrodynamics. *Ph.D. Thesis*.

[44] Shizuta Y, Kawashima S. 1985 Systems of equations of hyperbolic-parabolic type with applications to the discrete Boltzmann equation. *Hokkaido Mathematical Journal* **14**, 249–275.

[45] Adogbo JP, Danchin R. 2025 Global existence for multi-dimensional partially diffusive systems. *Journal of Differential Equations* **445**, 51.

[46] Godlewski E, Raviart PA. 2021 *Numerical approximation of hyperbolic systems of conservation laws* vol. 118 *Applied Mathematical Sciences*. Springer-Verlag, New York. Second edition.

[47] Grava T, Klein C. 2007 Numerical solution of the small dispersion limit of Korteweg-de Vries and Whitham equations. *Communications on Pure and Applied Mathematics* **60**, 1623–1664.

[48] Gurevich A, Pitaevskii L. 1973 Nonstationary structure of a collisionless shock wave. *Zhurnal Eksperimentalnoi i Teoreticheskoi Fiziki* **65**, 590–604.

[49] Besse C, Gavrilyuk S, Kazakova M, Noble P. 2022 Perfectly matched layers methods for mixed hyperbolic-dispersive equations. *Water Waves* **4**, 313–343.

[50] Lefloch PG, Rohde C. 2000 High-order schemes, entropy inequalities, and nonclassical shocks. *SIAM Journal on Numerical Analysis* **37**, 2023–2060.

[51] Barthwal R, Dhaouadi F, Rohde C. 2026 A hyperbolic relaxation system for a class of fourth-order nonlinear partial differential equations. *In preparation*.
[52] Huang Q, Rohde C, Yong WA, Zhang R. 2025 A hyperbolic relaxation approximation of the incompressible Navier-Stokes equations with artificial compressibility. *Journal of Differential Equations* **438**, Paper No. 113339, 26.

[53] Keim J, Munz CD, Rohde C. 2023 A relaxation model for the non-isothermal Navier-Stokes-Korteweg equations in confined domains. *Journal of Computational Physics* **474**, Paper No. 111830, 28.