\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 106, pp. 1--16.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/106\hfil Asymptotic waves in Jeffrey media] {Asymptotic dissipative waves in Jeffrey media from the point of view of double-scale method} \author[L. Restuccia, A. Georgescu \hfil EJDE-2017/106\hfilneg] {Liliana Restuccia, Adelina Georgescu} \address{Liliana Restuccia \newline Department of Mathematical and Computer Sciences, Physical Sciences and Earth Sciences, University of Messina, 98166 Messina, Italy} \email{lrestuccia@unime.it} \address{Adelina Georgescu \newline Academy of Romanian Scientists, 050094 Bucharest, Romania} \email{lrestuccia@unime.it} \dedicatory{Communicated by Giovanni Molica Bisci} \thanks{Submitted January 16, 2017. Published April 21, 2017.} \subjclass[2010]{34K35, 34A08, 47045, 65J20} \keywords{Double-scale method; non-linear dissipative waves; \hfill\break\indent asymptotic method; rheological media} \begin{abstract} This article studies previous results on nonlinear dissipative waves in Jeffrey media (viscoanelastic media without memory of order one) done by the first author, from the point of view of double scale method. For these media the equations of motion include second order derivative terms multiplied by a very small parameter. The physical meaning of a new (fast) variable, related to the surfaces across which the solutions or/and some of their derivatives vary steeply, is explained. The three-dimensional case is considered, that contains as a particular case an one-dimensional application worked out in a previous paper. Some known results are revised, other ones are derived and original. The thermodynamic models for Jeffrey media have applications in rheology and in other technological fields of applied sciences. \end{abstract} \maketitle \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} The theoretical interest in nonlinear waves was manifest as early as the years '50 and '60 of the last century and a lot of applications to various branches of physics were worked out \cite{b2, b3, c1, d1, f1, j1, j2, j3, j4, l2, l3}. In the context of rheological media studies on non-linear waves were carried out in \cite{c3,c4,c5}. In previous papers \cite{g4,r1} we sketched out the general use of the double-scale method to nonlinear hyperbolic partial differential equations (PDEs) to study the asymptotic waves and as an application the model governing the motion of inelastic media without shape and bulk memory (Maxwell media) was studied. The mathematical aspects involved in the study of asymptotic waves belong to singular perturbation theory, namely the double-scale method \cite{b1,c6,e1, g1, k1,k6,l1,m1,o1,s1,v1,w1}. The multiple-scale method, and, in particular, the double-scale approach, is appropriate to phenomena which possess qualitatively distinct aspects at various scales. For instance, at some well-determined times or space coordinates, the charac\-teristics of the motion vary steeply, while at larger scale the characteristics are slow and describe another type of motion. In addition, the scales are defined by some small parameters. In \cite{g2, g3} applications of the double scale method were given in the one-dimensional case to study of non-linear asymptotic waves in vi\-sco\-anelastic media with memory and in Jeffrey media, respectively. In this article the three-dimensional case is considered and the double scale method \cite{g1} is applied to investigate non-linear dissipative waves in isotropic vi\-sco\-anelastic media without memory of order one in which a viscous flow phenomenon occurs (Jeffreys media), that were studied in \cite{c4} by the first author (L. R.) in more classical way, following the methodologies developed in \cite{b3} and generalized in \cite{f1}. Only shear phenomena are taken into consideration and the hydrostatic pressure is assumed constant and uniform. For these media the equations of motion include second order derivative terms, multiplied by a very small parameter, that play a very important role because they usually have a balancing effect on the non-linear steepening of waves. In Section 2 the application of double-scale method to non-linear PDES is treated, the various steps in applying this method are introduced and the asymptotic approximations of first and second order are derived. In Section 3 the equations governing the motion of Jeffreys media are introduced and the mecha\-ni\-cal rela\-xation equation for these media is described in the framework of classical irreversible thermodynamics (TIP) with internal variables \cite{c2,k2, k3, k4, k5}. Furthermore, a matrix formulation of these equations is given. In Sections 4 and 5 the propagation into an uniform unperturbed state and the derivation of the asymptotic approximation of first order of the wave front and of the solution are worked out. Some known results are revised, other ones are derived and completely new. \section{Application of double-scale method to nonlinear PDEs} In \cite{c4} it was shown that the motion of viscoanelastic media without memory, in the case where only shear phenomena are taken into consideration and the hydrostatic pressure is constant and uniform, is described by a system of nonlinear partial differential equations (PDEs) having the matrix form \begin{equation} \label{2.1} \begin{gathered} \mathbf{A}^\alpha(\mathbf{U})\mathbf{U}_\alpha + \omega^ {-1} \Big[\mathbf{H}^k \frac{\partial ^2\mathbf{U}}{\partial t\partial x^k} + \mathbf{H}^{ik} \frac{\partial ^2\mathbf{U}}{\partial x^i\partial x^k}\Big] =\mathbf{B}(\mathbf{U}), \\ \alpha = 0,1,2,3 ;\; i,k = 1, 2, 3, \end{gathered} \end{equation} where $ x^0 = t$ (time), $x^1, x^2, x^3$ are the space coordinates, $\mathbf{U}$ is the vector of the unknown functions (which depend on $x^\alpha$), $\mathbf{U}_\alpha = \frac{\partial \mathbf{U}}{\partial x^\alpha}$, $\mathbf{A}^\alpha$, $\mathbf{H}^k$, $\mathbf{H}^{ik}$ are appropriate matrices $9\times9$, and \begin{equation} \label{2.1a} \mathbf{A}^{\alpha}(\mathbf{U})\mathbf{U}_{\alpha} =\mathbf{B}(\mathbf{U}) \end{equation} is the associated system of nonlinear hyperbolic PDEs. The system of PDEs \eqref{2.1} includes terms containing second order derivatives multiplied by a very small parame\-ter. These terms play a very important role because they usually have a balancing effect on the non-linear steepening of waves. In \cite{t1}, using \eqref{2.1}, the pro\-pa\-gation of linear acoustic waves was considered and the velocity and attenuation of the waves were investigated. \par In \cite{c4} the non-linear dissipative waves were worked out (see \cite{b2, b3, c1, d1, f1, j1, j2, j3, j4, l2, l3}) and, in particular, a method developed by Boillat \cite{b3} and generalized by Fusco \cite{f1} was applied to construct asymptotic approximations of order 1 of solutions of the system of equations \eqref{2.1}. \subsection{Asymptotic waves from the point of view of double scale-method} The smooth solutions of systems of type \eqref{2.1} (or type \eqref{2.1a}) that present a steep variation in the normal direction to the associated wavefront are called \textit{asymptotic waves}. Then, there exists a family of hypersurfaces $S$ (defined by the equation $\varphi(x^\alpha)=0$) moving in the Euclidean space $E^{3+1}$ (consisting of points of coordinates $x^\alpha$, $\alpha = 0,1,2,3$, or, equivalently of the time $t=x^0$ and the space coordinates $x^i$, $i=1,2,3$) having equation \begin{equation} \varphi(t,x^i)=\bar{\xi}={\rm const.}, \end{equation} such that the solutions $\mathbf{U}$ or/and some of their derivatives vary steeply across $S$ while along \textit{S } their variation is slow \cite{g4}. From the double scale method point of view this means that around $S$ there exist \textit{ asymptotic internal layers} (see \cite{g1}) such that \textit{the order of magnitude (i.e. the scale)} of the solutions or/and of some of their derivatives inside these layers and far away from them differs very much. In systems of equations of type \eqref{2.1} the coefficient $\omega^ {-1}$ is the small parameter, that is associated with the order of magnitude of the interior layer. Therefore, it is natural to introduce a new independent variable $\xi$, related to the hypersurfaces $S$, \begin{equation} \xi=\omega\bar{\xi}=\omega\varphi(t,x^i), \end{equation} where $\xi=\frac{\varphi(t,x^i)}{\omega^{-1}}$ is asymptotically fixed, i.e. $\xi=Ord (1)$ as $\omega^{-1} \to 0$, and $\omega\gg1$ is a very large parameter, to assume that the solution depends on the old as well as the new variable, i.e. $\mathbf{U}=\mathbf{U}(x^{\alpha},\xi)$, and to consider that $x^\alpha$ and $\xi$ are independent. Taking into account that $\mathbf{U}$ is sufficiently smooth, hence it has sufficiently many bounded derivatives, it follows that, except for the terms containing $\omega$, all other terms are asymptotically fixed and the computation can proceed formally. In this way, if $x^\alpha=x^\alpha(s)$ are the parametric equations of a curve $C$ in $E^{3+1}$, we have $$ \frac{d\mathbf{U}}{ds}=\omega\frac{\partial\mathbf{U}}{\partial\xi} \frac{\partial\varphi}{\partial s} +\frac{\partial \mathbf{U}} {\partial x^\alpha}\frac{dx^\alpha}{ds} $$ (where the dummy index convention is understood). This relation shows that, indeed, along $C$, $\mathbf{U}$ does not vary too much if $C$ belongs to the hypersurface $S$ (in this case $\frac{d\varphi}{ds}= 0$) but has a large variation if $C$ is not situated on $S$. For these reasons $\xi$ is referred to as the fast variable. Recall that the wavefront $\varphi$ is still an unknown function. To determine it, we recall its equation is $\varphi(t, x^1, x^2, x^3) = 0$. This implies that along the wavefront we have $\frac{d\varphi}{dt} = 0$, implying $\frac{\partial \varphi}{\partial t} + \mathbf{v}\cdot\operatorname{grad}\varphi=0$, or equivalently, $ \frac{\frac{\partial \varphi}{\partial t}}{|\operatorname{grad}\varphi|} + \mathbf{v}\cdot \frac{\operatorname{grad}\varphi}{|\operatorname{grad}\varphi|}=0$. Obviously, \begin{equation} \label{n} \frac{\operatorname{grad}\varphi}{|\operatorname{grad}\varphi|}=\mathbf{n}, \end{equation} such that the previous equality reads \begin{equation} \label{mt} \frac{\frac{\partial \varphi}{\partial t}}{|\operatorname{grad}\varphi|}+ \mathbf{v}\cdot \mathbf{n}=0. \end{equation} We introduce the notation \begin{equation} \label{lambda} \lambda=-\frac{\frac{\partial \varphi}{\partial t}}{|\operatorname{grad}\varphi|}, \end{equation} so that \begin{equation} \lambda(\mathbf{U},\mathbf{n})=\mathbf{v} \cdot \mathbf{n}, \end{equation} where $\lambda$ is called the \textit{velocity normal to the progressive wave}, being $\mathbf{n}$ the unit vector normal to the wave front. \par Following the general theory \cite{b3} we introduce the quantity \begin{equation}\label{3.21a} \Psi(\mathbf{U}, \Phi_\alpha) = \varphi_t + |\operatorname{grad}\varphi| \lambda(\mathbf{U}, \mathbf{n}). \end{equation} The characteristic equations for \eqref{3.21a} are \begin{equation}\label{3.27} \frac{d x^\alpha}{d \sigma} = \frac{\partial \Psi}{\partial \Phi_\alpha}, \quad \frac{d \Phi_\alpha}{d \sigma} = - \frac{\partial \Psi}{\partial x^\alpha} \quad (\alpha = 0,1,2,3), \end{equation} where $\sigma$ is the time along the rays. The $i$-th component of the radial velocity $ \mathbf{\Lambda}$ is defined by \begin{equation}\label{3.22} \Lambda_i (\mathbf{U}, \mathbf{n}) \equiv \frac{d x^i}{d \sigma} =\frac{\partial \Psi}{\partial \phi_i} = \lambda n_i + \frac{\partial \lambda}{\partial n_i} - \Big(\mathbf{n} \cdot \frac{\partial \lambda}{\partial \mathbf{n}}\Big) n_i = \lambda n_i + v_i - (n_k v_k) n_i, \end{equation} for $i = 1,2,3$. Hence, \begin{equation}\label{3.22a} \mathbf{\Lambda}(\mathbf{\mathbf{U}}, \mathbf{n}) = \mathbf{v} -(v_n - \lambda) \mathbf{n}. \end{equation} The theory in \cite{b3} enables us to deduce the equation for $\varphi$ by using \eqref{mt}. \subsection{Outline of the various steps in applying the double-scale method} The first step of the double-scale method consists in expressing the derivatives with respect to $x^\alpha$, $\frac{\partial}{\partial x^\alpha}$, in terms of the derivatives with respect to $x^\alpha$ and $\xi$, i.e. \[ \frac{\partial}{\partial x^\alpha} =\frac{\partial}{\partial x^\alpha}+\frac{\partial}{\partial \xi}\frac{\partial \xi}{\partial x^\alpha}= \frac{\partial}{\partial x^\alpha}+\omega\frac{\partial}{\partial \xi}\frac{\partial \varphi}{\partial x^\alpha}, \] so that the derivative $\mathbf{U}_{\alpha} = \frac{\partial \mathbf{U}}{\partial x^\alpha}$ has the form \begin{equation} \label{deriv} \frac{\partial \mathbf{U}}{\partial x^\alpha}\sim \omega^{-1} \Big( \frac{\partial \mathbf{U}^1} {\partial x^\alpha} + \omega\frac{\partial \mathbf{U}^1}{\partial \xi} \frac{\partial \varphi}{\partial x^\alpha}\Big) + \omega^{-1}\frac{\partial \mathbf{U}^2}{\partial \xi} \frac{\partial \varphi}{\partial x^\alpha} + O(\omega^{-2}), \quad \text{as } \omega^{-1} \to 0. \end{equation} In the second step, we look for the solution of the equations as an asymptotic series of powers of the small parameter, say $\epsilon$, namely with respect to the asymptotic sequence $\{1, \varepsilon^{a+1}, \varepsilon^{a+2},\dots ,\}$ or $ \{1, \varepsilon^{\frac{1}{p}}, \varepsilon^{\frac{2}{p}}, \dots , \}$, as $\varepsilon\to 0$. In \cite{c3,c4,c5} it is considered $p = 1$, and $\varepsilon = \omega^{-1}$, such that $\mathbf{U}(x^{\alpha},\xi)$ is written as an asymptotic power series of the small parameter $\omega^{-1}$, i.e. with respect to the asymptotic sequence $1,\omega^{-1}, \omega^{-2},\dots $. Since $\omega^{-1} \to 0$, and the $\mathbf{U}^i$ ($i= 1, 2,\dots$) are functions of $x^{\alpha}$ and $\xi$, it follows that \begin{equation} \label{2.5} \mathbf{U}(x^\alpha , \xi ) \sim \mathbf{U}^0 (x^\alpha , \xi ) + \omega^{-1}\mathbf{U}^1 (x^\alpha , \xi ) + O(\omega^{-2}),\quad \text{as } \omega^{-1}\to 0, \end{equation} where we have assumed that the first approximation $\mathbf{U}^0$ is constant. In \eqref{2.5} $\mathbf{U}^0 (x^\alpha , \xi )$ is a known solution \cite{f1} of \begin{equation} \label{2.1c} \mathbf{A}^{\alpha}(\mathbf{U}^0)\mathbf{U}_{\alpha}(\mathbf{U}^0) =\mathbf{B}(\mathbf{U}^0), \end{equation} where $\mathbf{U}^0$ is taken as the initial unperturbed state (see \eqref{2.1a}). Then, taking into account the form of $\mathbf{A}^\alpha$, $\mathbf H^k$, $\mathbf H^{i k}$ and $\mathbf{B}$, the following asymptotic expansions are deduced: \begin{gather}\label{3.4} \mathbf{A}^\alpha (\mathbf{U}) \sim \mathbf{A}^\alpha (\mathbf{U}^0) + \frac{1}{\omega} \nabla \mathbf{A}^\alpha (\mathbf{U}^0) \mathbf{U}^1 + O \big(\frac{1}{\omega^2}\big), \quad \text{as } \omega^{-1} \to 0,\\ \label{3.5} \mathbf H^k(\mathbf{U}) \sim \mathbf H^k(\mathbf{U}^0) + \frac{1}{\omega}\nabla \mathbf H^k(\mathbf{U}^0) \mathbf{U}^1 + O \big(\frac{1}{\omega^2}\big), \quad \text{as } \omega^{-1} \to 0 \;(k=1,2,3), \\ \label{3.5b} \mathbf H^{i k}(\mathbf{U}) \sim \mathbf H^{i k}(\mathbf{U}^0) + \frac{1}{\omega}\nabla \mathbf H^{i k}(\mathbf{U}^0) \mathbf{U}^1 + O \big(\frac{1}{\omega^2}\big), \quad\text{as } \omega^{-1} \to 0 \; (i, k = 1,2,3), \\ \label{3.6} \mathbf{B}(\mathbf{U}) \sim \mathbf{B}(\mathbf{U}^0) + \frac{1}{\omega}\nabla \mathbf{B}(\mathbf{U}^0)\mathbf{U}^1 + O \big(\frac{1}{\omega^2}\big), \quad \text{as } \omega^{-1} \to 0,\; \end{gather} where $\nabla=\frac{\partial}{\partial\mathbf U}$. The last point of the method consists in introducing the asymptotic expansions \eqref{deriv}--\eqref{3.6} in \eqref{2.1} and matching the obtained series. It follows that \begin{gather}\label{Uuno} (\mathbf{A}^\alpha)_0\Phi_\alpha\frac{\partial\mathbf{U}^1}{\partial\xi} =\mathbf 0\quad(\alpha=0,1,2,3), \\ \label{3.8} \begin{aligned} &(\mathbf{A}^\alpha)_0 \Big( \Phi_\alpha\frac{\partial \mathbf{U}^2}{\partial \xi}\Big) \\ &= - \Big[(\mathbf{A}^\alpha)_0 \frac{\partial \mathbf{U}^1}{\partial x^\alpha} + (\nabla \mathbf{A}^\alpha)_0\mathbf{U}^1 \Big(\Phi_\alpha\frac{\partial \mathbf{U}^1}{\partial \xi}\Big) +(\mathbf{H}^{k})_0 \Phi_0 \Phi_k\frac{\partial^2\mathbf U^1}{\partial\xi^2}\\ &\quad + (\mathbf{H}^{i k})_0 \Phi_i \Phi_k\frac{\partial^2\mathbf U^1}{\partial\xi^2} -(\nabla \mathbf{B})_0\mathbf{U}^1 \Big], \end{aligned} \end{gather} where $\Phi_\alpha=\frac{\partial\varphi}{\partial x^\alpha}$ $(\Phi_k=\frac{\partial\varphi}{\partial x^k},\; k=1,2,3)$ and the symbol $(\dots )_0$ indicates that the quantities are calculated in $\mathbf U^0$. Equation \eqref{Uuno} is linear in $\mathbf{U}^1$, while \eqref{3.8} is affine in $\mathbf{U}^2$. Of course, equations of asymptotic approximations of higher order can be written and they are affine, but their solutions are very difficult. Just to solve the linear equation \eqref{Uuno}, in \cite{c4} the methods developed in \cite{b3,f1} were applied obtaining many results that we present in Sections 4, 5. \section{Equations governing the motion of Jeffreys media} In \cite{k2}, in the framework of the classical irreversible thermodynamics with internal variables, Kluitenberg developed a theory for mechanical relaxation phenomena in rheological media, assuming that, when several microscopic phenomena give rise to inelastic deformation, the tensor of the total strain $\varepsilon_{\alpha\beta}$ can be split in two parts: $\varepsilon_{\alpha\beta}=\varepsilon_{\alpha\beta}^{el} + \varepsilon_{\alpha\beta}^{in}$, where the tensors $\varepsilon_{\alpha\beta}^{el}$ and $\varepsilon_{\alpha\beta}^{in}$ describe the ela\-stic and inelastic strains, respectively. It is supposed that the inelastic deformation $\varepsilon_{\alpha\beta}^{in}$ is due to the defects of the lattice (slip, dislocations,\dots ) and to the influence of microscopic stress fields, surrounding imperfections in the medium, and that it can be split in $n$ contributions $\varepsilon_{\alpha\beta}^{(k)}$ $(k=1,2,\dots ,n)$: $\varepsilon_{\alpha\beta}^{in}=\sum_{k=1}^n\varepsilon_{\alpha\beta}^{(k)}$ (being $n$ arbitrary), that are introduced as \textit{internal variables} in the thermodynamical state vector. Furthermore, in the theory it is assumed that the gradient of the displacement field is small and this implies that the deformations and the rotations are small. In this case the strain tensor $\varepsilon_{ik}$ is defined by $\varepsilon_{ik} =\frac{1}{2}\left(\frac{\partial u_i}{\partial x^k} + \frac{\partial u_k}{\partial x^i}\right)$ ($i,k = 1,2,\dots ,n$), where $u_i$ is the $i$-th component of the displacement field $\mathbf{U}$ and $x^i$ is the $i$-th component of the position vector $\mathbf{x}$ in Eulerian coordinates in a Cartesian reference frame. \par In \cite{k2} for shear phenomena in isotropic viscoanelastic media without memo\-ry, where the hydrostatic pressure is assumed constant and uniform, the follo\-wing mechanical relaxation equation between the deviators $\tilde{\tau}_{ik}$ of the mechanical stress tensor and $\tilde{\epsilon}_{ik}$ of the strain tensor was given by \begin{equation} R^{(\tau)}_{(d)0}\tilde{\tau}_{ik} + \sum_{m = 1}^{n - 1} R^{(\tau)}_{(d)m}\frac{d^m}{dt^m}\tilde{\tau}_{ik} + \frac{d^n}{dt^n}\tilde{\tau}_{ik} = R^{(\epsilon)}_{(d)0}\tilde{\epsilon}_{ik} + \sum_{m = 1}^{n + 1} R^{(\epsilon)}_{(d)m}\frac{d^m}{dt^m}\tilde{\epsilon}_{ik} \end{equation} for $i,k = 1,2,3$, being the quantities $R^{(\tau)}_{(d)m}$ ($m = 0,1, \dots , n - 1$), $R^{(\epsilon)}_{(d)m}$ ($m = 0,1, \dots , n + 1$), algebraic functions of the coefficients occurring in the phenomenological equations and in the equations of state. Mechanical relaxation equations were derived from this more general mentioned above relation for Maxwell, Jeffreys, Burgers media and other special rheological materials (see \cite{c2,k2, k3, k4, k5}). Assuming that only one microscopic phenomenon gives rise to ine\-lastic strain, in the isotropic case and for shear phenomena, when the hydrostatic pressure is assumed constant and uniform, the relaxation equation describing the mechanical behaviour of viscoanelastic media without memory (Jeffreys media) of order one $(n = 1)$ (i.e. when only one internal variable of mechanical origin is taken into consideration), can be written in the following form \begin{equation}\label{2.4} R^{(\tau)}_{(d)0}\tilde{P}_{ik} + \frac{d}{d t}\tilde{P}_{ik} + R^{(\epsilon)}_{(d)1} \frac{d}{d t}\tilde{\epsilon}_{ik} + R^{(\epsilon)}_{(d)2} \frac{d^2}{d t^2}\tilde{\epsilon}_{ik}= 0, \end{equation} where $\tilde{P}_{ik}$ and $\tilde{\epsilon}_{ik}$ are the deviators of the mechanical pressure tensor $P_{ik}$ and of the strain tensor $\epsilon_{ik}$, respectively. We define $P_{ik}$ in terms of the symmetric Cauchy stress tensor $P_{ik} = -\tau_{ik}$ ($i,k = 1, 2, 3$) and the following quantities \begin{gather*} \tilde{P}_{ik} = P_{ik} - \frac{1}{3} P_{ss} \delta_{ik},\quad P =\frac{1}{3} P_{ss},\quad P_{ss} = tr P,\\ P_{ik} = \tilde{P}_{ik} + P \delta_{ik},\quad \tilde{P}_{ss}= 0, \end{gather*} where the hydrostatic pressure $P$ is the scalar part of the tensor $P_{ik}$. In \eqref{2.4} the coefficients satisfy the relations \begin{gather}\label{2.4a} R^{(\tau)}_{(d)0} = a^{(0,0)}\eta_{s}^{(1,1)} \geq 0, \\ \label{2.4b} R^{(\epsilon)}_{(d)1} = a^{(0,0)}\big[\left(1+\eta_{s}^{(0,1)}\right)^2 +\eta_{s}^{(0,0)}\eta_{s}^{(1,1)}\big]\geq 0, \\ \label{2.4c} R^{(\epsilon)}_{(d)2} = \eta_{s}^{(0,0)} \geq 0, \end{gather} where $a^{(0,0)}$ is a scalar constant which occurs in the equations of state, while the coefficients $\eta_s^{(0,0)}$, $\eta_s^{(0,1)}$ and $\eta_s^{(1,1)}$ are called \textit{phenomenological coefficients} and represent fluidities. The balance equations for the mass density and momentum in the case of Jeffreys media read \cite{c4} \begin{gather}\label{4.13b} \frac{\partial \rho}{\partial t} + \frac{\partial }{\partial x_i}(\rho v_i) = 0 \quad (i = 1,2,3), \\ \label{4.13} \rho \big(\frac{\partial }{\partial t} v_i + v_k \frac{\partial }{\partial x^k} v_i \big) + \frac{\partial }{\partial x^k} \tilde{P}_{ik} = 0, \end{gather} where $v_i = d u_i/dt$ is the $i$-th component of the velocity field and the force per unit mass is neglected. Thus, \eqref{2.4} becomes \begin{equation}\label{2.10} \begin{aligned} &R^{(\tau)}_{(d)0} \tilde{P}_{ik} + \frac{\partial}{\partial t}\tilde{P}_{ik} + v_p \frac{\partial}{\partial x^p}\tilde{P}_{ik} \\ &+ R^{(\epsilon)}_{(d)1}\Big[\frac{1}{2}\Big(\frac{\partial }{\partial x^k} v_i + \frac{\partial }{\partial x^i}v_k\Big) - \frac{1}{3}\frac{\partial}{\partial x^p} v_p \delta_{ik}\Big] \\ & + R^{(\epsilon)}_{(d)2}\Big[\frac{1}{2} \Big(\frac{\partial^2 }{\partial t \partial x^k} v_i + v_p \frac{\partial^2 }{\partial x^p \partial x^k}v_i + \frac{\partial^2 }{\partial t \partial x^i} v_k + v_p \frac{\partial^2 }{\partial x^p \partial x^i} v_k\Big) \Big] \\ &- R^{(\epsilon)}_{(d)2}\Big[\frac{1}{3} \Big( \frac{\partial^2 }{\partial t \partial x^p} v_p + v_q \frac{\partial^2 }{\partial x^q \partial x^p} v_p\Big) \delta_{ik} \Big] = 0 \quad (p, q=1, 2, 3), \end{aligned} \end{equation} where the relations $\frac{d \varepsilon_{ik}}{dt} = \frac{1}{2} \left(\frac{\partial v_i }{\partial x^k} + \frac{\partial v_k }{\partial x^i}\right)$ and $\frac{d}{dt} = \frac{\partial}{\partial t} + v_p\frac{\partial}{\partial x^p}$ are used. \subsection{Matrix form of equation system governing the motion of Jeffreys media} It is easy to see that the system of equations \eqref{4.13b}--\eqref{2.10}, together with \eqref{2.4a}--\eqref{2.4c}, takes the matrix form \eqref{2.1}, where $\mathbf{A}^0(\mathbf{U}) = \mathbf{I}$ is the identity matrix, $R^{(\epsilon)}_{(d)2} = \omega^{-1}R'^{(\epsilon)}_{(d)2}$, the matrices $\mathbf{A}^i$, $\mathbf{H}^i$ and $\mathbf{H}^{ik}$ $(i,k = 1,2,3)$ are appropriate $9\times 9$ matrices, given in the Appendix, and \begin{gather}\label{2.14} \mathbf{U} = (\rho, v_1, v_2, v_3, \tilde P_{11}, \tilde P_{12}, \tilde P_{13}, \tilde P_{22}, \tilde P_{23})^T, \\ \label{2.15} \mathbf{B} = (0, 0, 0, 0,\tilde P_{11}^*, \tilde P_{12}^*, \tilde P_{13}^*, \tilde P_{22}^*, \tilde P_{23}^*)^T, \\ \label{2.17} \tilde P_{ik}^* = - {R}_{(d)0}^{(\tau)} \tilde{P}_{ik} = - a^{(0,0)}\eta_{s}^{(1,1)} \tilde{P}_{ik}\quad (i,k = 1,2,3). \end{gather} The symbol $(\dots )^T$ means that $\mathbf{U}$ and $\mathbf{B}$ are column vectors of 9 components. \section{Propagation into a uniform unperturbed state} Let us consider an uniform unperturbed state in which $\mathbf{U}^0$, solution of \eqref{2.1c}, is \begin{equation}\label{3.1} \mathbf{U}^0 = (\rho^0, 0, 0, 0, 0, 0, 0, 0, 0)\quad (\rho^0 = {\rm const.}). \end{equation} If the quantities \eqref{n} and \eqref{lambda} are introduced in the expression \eqref{Uuno} we obtain \begin{equation}\label{3.1a} (A_{0n} - \lambda I)\frac{\partial \mathbf{U}^1}{\partial \xi} = 0, \end{equation} where $(\mathbf{A}_n)_0=A_{0n} $ and $\mathbf{A}_n (\mathbf{U}) = \mathbf{A}^i n_i$ is an appropriate $9\times 9$ matrix. Using expressions \eqref{4.16}--\eqref{4.18} in equation \eqref{3.1a}, $\mathbf{A}_n (U)$ assumes the form \begin{equation}\label{3.13} \begin{aligned} &\mathbf{A}_n (\mathbf{U}) = \mathbf{A}^i n_i \\ &= \begin{pmatrix} v_n & \rho n_1 & \rho n_2 & \rho n_3 & 0 & 0 & 0 & 0 & 0 \\ 0 & v_n & 0 & 0 & \frac{n_1}{\rho} & \frac{n_2}{\rho} & \frac{n_3}{\rho} & 0 & 0 \\ 0 & 0 & v_n & 0 & 0 & \frac{n_1}{\rho} & 0 & \frac{n_2}{\rho} & \frac{n_3}{\rho} \\ 0 & 0 & 0 & v_n & - \frac{n_3}{\rho} & 0 & \frac{n_1}{\rho} & - \frac{n_3}{\rho} & \frac{n_2}{\rho} \\ 0 & \frac{2}{3} R^{(\epsilon)}_{(d)1}n_1 & -\frac{1}{3} R^{(\epsilon)}_{(d)1}n_2 & -\frac{1}{3} R^{(\epsilon)}_{(d)1}n_3 & v_n & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_2 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_1 & 0 & 0 & v_n & 0 & 0 & 0 \\ 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_3 & 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_1 & 0 & 0 & v_n & 0 & 0 \\ 0 & - \frac{1}{3} R^{(\epsilon)}_{(d)1}n_1 & \frac{2}{3} R^{(\epsilon)}_{(d)1}n_2 & - \frac{1}{3} R^{(\epsilon)}_{(d)1}n_3 & 0 & 0 & 0 & v_n & 0 \\ 0 & 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_3 & \frac{1}{2} R^{(\epsilon)}_{(d)1}n_2 & 0 & 0 & 0 & 0 & v_n \\ \end{pmatrix} \end{aligned} \end{equation} In \cite{c4} it was found that the eigenvalues are real and the eigenvectors of the matrix $\mathbf{A}_n$ are linearly independent, so that the system of PDEs \eqref{2.1a} is hyperbolic (see \cite{g4} for the definition of hyperbolicity). The eigenvalues of $\mathbf{A}_n (\mathbf{U})$ are \begin{equation}\label{3.15} \lambda_1 = \mathbf{v }\cdot \mathbf{n }= v_n,\quad \lambda_2^{(\pm)} = v_n \pm \Big(\frac{R^{(\epsilon)}_{(d)1}}{2\rho}\Big)^{1/2}, \quad \lambda_3^{(\pm)} = v_n \pm \Big(\frac{2 R^{(\epsilon)}_{(d)1}}{3\rho}\Big)^{1/2}, \end{equation} where $\lambda_1$ has multiplicity 3, both $\lambda_2^{(+)}$ and $\lambda_2^{(-)}$ have multiplicity 2 and $\lambda_3^{(\pm)}$ are simple. Furthermore, by using \eqref{3.1a}, \eqref{3.15} and \eqref{3.13}, the eigenvectors $\mathbf{r}_1,\; \mathbf{r}_2^{(\pm)}\mathbf{r}_3^{(\pm)}$ corresponding to $\lambda_1$, $\lambda_2^{(\pm)}\lambda_3^{(\pm)}$ are \begin{gather} \label{r1} \begin{aligned} \mathbf{r}_1 &= \Big(1,0,0,0,\frac{-(n_2 + n_3)}{n_1}, 1,1, \frac{ n_1^2(n_3-n_2) + n_3^2(n_2 + n_3) }{ n_1(n_2^2 + n_3^2)}, \\ &\quad \frac{-(n_2 + n_3)( n_1^2 + n_2 n_3) }{ n_1(n_2^2 + n_3^2)} \Big), \end{aligned} \\ \label{r2} \begin{aligned} \mathbf{r}_2^{(\pm)} &= \Big(0, 1, \frac{ n_2^2 - 1}{n_1 n_2},\frac{ n_3}{n_1 }, \pm\sqrt{2\rho R^{(\epsilon)}_{(d)1}}n_1, \pm \sqrt{2\rho R^{(\epsilon)}_{(d)1}} \frac{ 2n_2^2 - 1}{2n_2}, \\ &\quad \pm\sqrt{2\rho R^{(\epsilon)}_{(d)1}}n_3, \pm \sqrt{2\rho R^{(\epsilon)}_{(d)1}} \frac{n_2^2 - 1}{n_1}, \pm \sqrt{2\rho R^{(\epsilon)}_{(d)1}} \frac{n_3(2n_2^2 - 1)}{2n_1n_2} \Big), \end{aligned} \\ \label{r3} \begin{aligned} \mathbf{r}_3^{(\pm)} &= \Big(\frac{3\rho R^{(\epsilon)}_{(d)1}}{2}\Big)^{1/2} \Big(\pm\frac{1}{R^{(\epsilon)}_{(d)1}}\rho, \pm\Big(\frac{2}{3\rho R^{(\epsilon)}_{(d)1}}\Big)^{1/2} \mathbf{n}, \\ &\quad \pm\frac{3 n_1^2 - 1}{3}, \pm n_1 n_2, \pm n_1 n_3, \pm \frac{3 n_2^2 - 1}{3}, \pm n_2 n_3 \Big). \end{aligned} \end{gather} In \cite{c4} it was seen that the discontinuity waves which are propagating with the velocities given by $\lambda_1$ and $\lambda_2^{(\pm)}$ satisfy the Lax-Boillat exceptionality condition \cite{b2} because $ \nabla\lambda_1\cdot\textbf{ r}_1 = 0$ and $ \nabla\lambda_2^{(\pm)}\cdot\textbf{ r}_2^{(\pm)} = 0$ (with $\nabla=\frac{\partial}{\partial\mathbf U}$), while the discontinuity waves whose velocities of propagation are $\lambda_3^{(\pm)}$ do not possess this property being $\nabla \lambda_3^{(+)}\cdot\mathbf{r}_3^{(+)}=\frac{1}{2}$. Thus, only for $\lambda_3^{(\pm)}$ our results are valid. Then, we fix our attention on $\lambda =\lambda_3^{(+)}$ which corresponds to a progressive fast longitudinal wave traveling to the right. The left eigenvector $\mathbf{l}_3^{(+)}$ corresponding to $\lambda_3^{(+)}$ is \begin{equation}\label{3.18} \mathbf{l}_3^{(+)} =\Big(\frac{3\rho}{2R^{(\epsilon)}_{(d)1}}\Big)^{1/2} \Big[ 0, \Big(\frac{R^{(\epsilon)}_{(d)1}}{3\rho}\Big)^{1/2} \frac{\mathbf{n}}{n_3}, \frac{n_1^2 - n_3^2}{\rho n_3}, \frac{2n_1 n_2}{\rho n_3}, \frac{2 n_1}{\rho}, \frac{n_2^2 - n_3^2}{\rho n_3}, \frac{2 n_2}{\rho}\Big] . \end{equation} Furthermore, the left eigenvector $\mathbf{l}_3^{(+)}$ and the right eigenvector $\mathbf{r}_3^{(+)}$ satisfy the relation \begin{equation}\label{lr} \mathbf{l}_3^{(+)} \cdot \mathbf{r}_3^{(+)} = \frac{2}{n_3}. \end{equation} From \eqref{3.1a} we have \begin{equation}\label{3.20} \mathbf{U}^1(x^\alpha, \xi) = u(x^\alpha, \xi) \mathbf{r}_3^{(+)} (\mathbf{U}^0, \mathbf{n}) + \mathbf{v}^1(x^\alpha), \end{equation} where $u$ is a scalar function to be determined and $\mathbf{v}^1$ is an arbitrary function of integration which can be taken as zero, without loss of generality (see \cite{b3}). It may be observed that in \eqref{3.20} $u$ gives rise to the phenomenon of the distortion of the signals and this term governs the first-order perturbation obeying a non-linear partial differential equation (see Section 4). We conclude this section by showing how the wave front $\varphi(t, x^1, x^2, x^3) = 0$ can be determined (see \cite{c4}). From equation \eqref{3.22a}, for $\lambda = \lambda_3^{(+)}$ we obtain \begin{equation}\label{3.24} \mathbf{ \Lambda}(\mathbf{U}, \mathbf{n}) = \mathbf{v} +\Big(\frac{2R^{(\epsilon)}_{(d)1}}{3\rho}\Big)^{1/2} \mathbf{n}. \end{equation} Since we are considering the propagation into a uniform unperturbed state it is known \cite{b3} that the wave front $\varphi$ satisfies the partial differential equation \begin{equation}\label{3.25} \Psi(\mathbf U^0, \Phi_\alpha) = \varphi_t + |\operatorname{grad} \varphi|\lambda_3^{(+)}(\mathbf U^0, \mathbf n) = \Psi^0 = 0, \end{equation} and so \begin{equation}\label{3.26} \Lambda_i(\mathbf{U}^0, \mathbf{n}) = \frac{\partial \Psi^0}{\partial \Phi_i} \quad (i = 1,2,3). \end{equation} The characteristic equations for \eqref{3.25} are \begin{equation}\label{3.27b} \frac{d x^\alpha}{d \sigma} = \frac{\partial \Psi^0}{\partial \Phi_\alpha}, \quad \frac{d \Phi_\alpha}{d \sigma} = - \frac{\partial \Psi^0}{\partial x^\alpha} \quad (\alpha = 0,1,2,3), \end{equation} where $\sigma$ is the time along the rays. From \eqref{3.24} and \eqref{3.26} the expression of the radial velocity along the rays is \begin{equation}\label{3.29} \Lambda_i^0(\mathbf{U}^0, \mathbf{n}^0) = \Big(\frac{2 R^{(\epsilon)}_{(d)1}}{3\rho^0}\Big)^{1/2} {n}_i^0, \end{equation} where $\mathbf{n}^0$ is a constant value of $\mathbf{n}$. By integration of \eqref{3.27} one obtains \begin{equation}\label{3.30} x^0 = t = \sigma, \quad x^i = (x^i)^0 + \Lambda_i^0(\mathbf{U}^0, \mathbf{n}^0)t, \;\; \end{equation} with $(x^i)^0 = (x^i)_{t = 0}$ ($i = 1,2,3$). If we denote by $\varphi^0$ the given initial surface, we have $$ (\varphi)_{t = 0} = \varphi^0 \left[(x^i)^0\right] $$ and $\mathbf{n}^0$ represents the normal vector at the point $(x^i)^0$ defined by $$ \mathbf{n}^0 = \Big(\frac{\operatorname{grad}\varphi}{|\operatorname{grad}\varphi|} \Big)_{t = 0} = \frac{grad^0 \varphi^0}{|grad^0 \varphi^0|}, \quad \text{where } (grad^0)_i \equiv \frac{\partial }{\partial (x^i)^0} \quad (i = 1,2,3). $$ Then, $\mathbf{x} = \mathbf{x}|_{t=0} + \mathbf{\Lambda}^0 t$ and since the Jacobian $J$ of the transformation $\mathbf{x} \to \mathbf{x}|_{t=0}$ is nonvanishing, i.e. $$ J = \det|\delta_{ik} + \frac{ \partial \Lambda^0_k}{ \partial (x^i)^0}t |\neq 0 \quad (i,k = 1,2,3), $$ $(x^i)^0$ can be deduced from equations \eqref{3.30}$_2$ and $ \varphi$ in the first approximation takes the following form \begin{equation}\label{3.36} \varphi(t, x^i) = \varphi^0(x^i - \Lambda_i^0 t) = \varphi^0\Big(x^i -\Big(\frac{2 R^{(\epsilon)}_{(d)1}}{3\rho^0}\Big)^{1/2} {n}_i^0t\Big). \end{equation} \section{First approximation of wavefront and of $\mathbf{U}$} Using \eqref{3.8} and \eqref{3.20} (see \cite{b3} and \cite{c4}), the following equation for $u(x^\alpha, \xi)$ can be obtained: \begin{equation}\label{4.1} \frac{\partial u}{\partial \sigma} + \left(\nabla \Psi \cdot \mathbf{r}_3^{(+)}\right)_0 u \frac{\partial u}{\partial \xi} + \frac{1}{\sqrt{J}}\frac{\partial \sqrt{J}}{\partial \sigma} u +\mu^0\frac{\partial^2 u}{\partial \xi^2} = \nu^0 u, \end{equation} where \begin{gather}\label{4.3} (\nabla \Psi \cdot \mathbf{r}_3^{(+)})_0= (|\operatorname{grad}\varphi|)_0 \big(\nabla \lambda_3^{(+)} \cdot \mathbf{r}_3^{(+)}\big)_0, \\ \mu^0=\frac{\big[\mathbf l_3^{(+)}\cdot\big(\mathbf H^{k} \frac{\partial \varphi}{\partial t} \frac{\partial \varphi}{\partial x^k} + \mathbf H^{i k} \frac{\partial \varphi}{\partial x^i} \frac{\partial \varphi}{\partial x^k}\big)\big(\mathbf r_3^{(+)}\big)\big]_0} {\big(\mathbf l_3^{(+)}\cdot\mathbf r_3^{(+)}\big)_0}, \\ \label{4.4} \nu^0 = \frac{\big(\mathbf{l}_3^{(+)}\cdot \nabla \mathbf{B } \mathbf{r}_3^{(+)}\big)_0} {\big(\mathbf{l}_3^{(+)} \cdot \mathbf{r}_3^{(+)}\big)_0}. \end{gather} Straightforward computations give \begin{equation}\label{4.6aa} \big(\nabla \Psi \cdot \mathbf{r}_3^{(+)}\big)_0 = \frac{1}{2}(|\operatorname{grad}\varphi|)_0, \end{equation} where $$ \nabla \lambda_3^{(+)} =\frac{\partial\lambda_3^{(+)}}{\partial \mathbf{U}}\equiv \Big[ -\frac{1}{2\rho} \Big(\frac{2 R^{(\epsilon)}_{(d)1} }{3 \rho}\Big)^{1/2}, n_1, n_2, n_3, 0, 0, 0, 0, 0, 0\Big] $$ and \begin{gather}\label{4.6a} \nabla \lambda_3^{(+)}\cdot\mathbf{r}_3^{(+)}=\frac{1}{2}, \\ \label{4.6} \big(\mathbf{l}_3^{(+)}\cdot \nabla \mathbf{B }\; \mathbf{r}_3^{(+)}\big)_0 = -\frac{R^{(\tau)}_{(d)0}}{n_3},\quad \nu^0 = - \frac{R^{(\tau)}_{(d)0}}{2} ,\quad \mu^0 = \frac{\big(\frac{\partial \varphi}{\partial t} |\operatorname{grad}\varphi|\big)_0 R'^{(\epsilon)}_{(d)2}} {\sqrt{6\rho^0 R^{(\epsilon)}_{(d)1}}}, \end{gather} where we used equations \eqref{lr} and \eqref{0,1}--\eqref{3,3}. By using the transformation of variables (see \cite{f1}) \begin{equation}\label{4.7} u = \frac{v}{\sqrt{J}} e^{w}, \quad \kappa = \int_0^\sigma \frac{e^w}{\sqrt{J}} \left( \nabla \Psi \cdot \mathbf{r}_3^{(+)}\right)_0 d\sigma, \quad \text{with } w = \int_0^\sigma \nu^0 d\sigma, \end{equation} equation \eqref{4.1} can be reduced to an equation of the type \begin{equation}\label{4.9} \frac{\partial v}{\partial \kappa} + v \frac{\partial v}{\partial \xi}+\hat{\mu}^0\frac{\partial^2 v}{\partial \xi^2} = 0,\quad \text{with } \hat{\mu}^0 = \frac{\mu^0\sqrt{J}e^{-w}}{\big( \nabla \Psi \cdot \mathbf{r}_3^{(+)} \big)_0}, \end{equation} which is similar to Burger's equation and is valid along the characteristic rays. Using the obtained results \eqref{4.6aa}--\eqref{4.6}, $\kappa$, $w$ and $\hat{\mu}^0$ are given by $$ \kappa = \int_0^\sigma \frac{1}{2}(|\operatorname{grad}\varphi|)_0 \frac{e^{w}}{\sqrt{J}} d\sigma, \quad w = \nu^0\sigma\quad \text{and}\quad \hat{\mu}^0 = \frac{2\mu^0\sqrt{J}e^{ - \nu^0 \sigma}} {(|\operatorname{grad}\varphi|)^0}. $$ Equation \eqref{4.9}$_1$ can be reduced to the semilinear heat equation \cite{h1} \begin{equation} \frac{\partial h}{\partial \kappa}=\hat\mu^0\frac{\partial^2 h}{\partial\xi^2}-h \log\frac{h}{\hat\mu^0}\frac{d\hat\mu^0}{d\kappa}, \end{equation} for which the solution is known, using the following Hopf transformation \cite{h1} \begin{equation} v(\xi,\kappa)=\hat\mu_0\frac{\partial}{\partial\xi}\log h(\xi,\kappa). \end{equation} \subsection*{Conclusions} In this article we presented a system of non-linear hyperbolic PDEs describing isotropic visco-inelastic media without memory of order one (Jeffrey media). Because a theory has an added value if we test it from the mathematical point of view, we manage to find possible solutions of the obtained thermodynamic model. But since the achievement of a closed-form solution of nonlinear PDEs is rare we look for the solution in the form of an asymptotic sequence of powers of some small parameter, which is related to the thickness of internal layers, across which the solutions or/and some of their derivatives varies steeply. We describe the various steps in applying the double scale method and we study the propagation of a particular solution into a uniform unperturbed state, obtaining the first approximation equation of the wave front of the solution. The three-dimensional case is treated. The thermodynamic models for Jeffrey media may have relevance in many fundamental technological sectors and in particular in rheology. The asymptotic waves for the media under consideration were studied by the first author (L. R.) in a more classical way in previous paper. \section{Appendix} The system of equations \eqref{4.13b}--\eqref{2.10} takes the matrix form \eqref{2.1}, where the matrices $\mathbf{A}^i$ ($i = 1,2,3$) have the form \begin{gather}\label{4.16} \mathbf{A}^1 = \begin{pmatrix} v_1 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & v_1 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 & 0 \\ 0 & 0 & v_1 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 \\ 0 & 0 & 0 & v_1 & 0 & 0 & \frac{1}{\rho} & 0 & 0 \\ 0 & \frac{2}{3}R^{(\epsilon)}_{(d)1} & 0 & 0 & v_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R^{(\epsilon)}_{(d)1}& 0 & 0 & v_1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R^{(\epsilon)}_{(d)1}& 0 & 0 & v_1 & 0 & 0 \\ 0 & - \frac{1}{3}R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & 0 & 0 & v_1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & v_1 \end{pmatrix}, \\ \label{4.17} \mathbf{A}^2 = \begin{pmatrix} v_2 & 0 & \rho & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & v_2 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 & 0 \\ 0 & 0 & v_2 & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 \\ 0 & 0 & 0 & v_2 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \\ 0 & 0 & -\frac{1}{3}R^{(\epsilon)}_{(d)1} & 0 & v_2 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & v_2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0& 0 & 0 & v_2 & 0 & 0 \\ 0 & 0 & \frac{2}{3}R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & 0 & v_2 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & 0 & v_2 \end{pmatrix}, \\ \label{4.18} \mathbf{A}^3 = \begin{pmatrix} v_3 & 0 & 0 & \rho & 0 & 0 & 0 & 0 & 0 \\ 0 & v_3 & 0 & 0 & 0 & 0 & \frac{1}{\rho} & 0 & 0 \\ 0 & 0 & v_3 & 0 & 0 & 0 & 0 & 0 & \frac{1}{\rho} \\ 0 & 0 & 0 & v_3 & - \frac{1}{\rho} & 0 & 0 & - \frac{1}{\rho} & 0 \\ 0 & 0 & 0 & -\frac{1}{3}R^{(\epsilon)}_{(d)1} & v_3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & v_3 & 0 & 0 & 0 \\ 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & 0 & v_3 & 0 & 0 \\ 0 & 0 & 0 & - \frac{1}{3} R^{(\epsilon)}_{(d)1} & 0 & 0 & 0 & v_3 & 0 \\ 0 & 0 & \frac{1}{2} R^{(\epsilon)}_{(d)1}& 0 & 0 & 0 & 0 & 0 & v_3 \end{pmatrix}. \end{gather} Furthermore, in the system of equations \eqref{2.1}, the matrices $\mathbf{H}^{i}$ and $\mathbf{H}^{ik}$ ($i,k = 1,2,3$) have the following form (with $R^{(\epsilon)}_{(d)2} = \omega^{-1} R'^{(\epsilon)}_{(d)2})$: \begin{gather}\label{0,1} \mathbf{H}^{1}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 \\ 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix},\\ \label{0,2} \mathbf{H}^{ 2}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{0,3} \mathbf{H}^{ 3}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{11} \mathbf{H}^{1 1}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 \\ 0 &-\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{22} \mathbf{H}^{1 2}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 &\frac{2}{3} R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{1,3} \mathbf{H}^{1 3}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_1& 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_1 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{2,1} \mathbf{H}^{2 1}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 \\ 0 &-\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{2,2} \mathbf{H}^{2 2}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 &-\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 &\frac{2}{3} R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{2'3'} \mathbf{H}^{2 3}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} v_2& 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & - \frac{1}{3}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_2 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{3,1} \mathbf{H}^{3 1}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{3,2} \mathbf{H}^{3 2}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 &-\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{2}{3}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \\ \label{3,3} \mathbf{H}^{3 3}= \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{1}{3}R'^{(\epsilon)}_{(d)2} v_3& 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2}R'^{(\epsilon)}_{(d)2} v_3 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}. \end{gather} \subsection*{Dedication and acknowledgements} L. R. dedicates this paper to the memory of Prof. Adelina Georgescu, great friend, eminent scientist and invaluable teacher. This work was supported by a grant of Messina University, by Accademia Peloritana dei Pericolanti and by Bonino-Pulejo Foundation. It was begun during a stay of Prof. A. Georgescu at Department of Mathematics at Messina University, in occasion of a series of her lectures on ``Asymptotic methods with applications to waves and shocks''. L. R. is deeply grateful for her precious advices and enlightening discussions regarding physical concepts and mathematical tools. \begin{thebibliography}{99} \bibitem{b1} N. N. Bogoliubov, Yu. A. Mitropolsk$\ddot{u}$, \emph{Asymptotic methods in the theory of nonlinear oscillations}, Gordon and Breach, New York, 1961. \bibitem{b2} G. Boillat; \emph{La propagation des ondes}, Gauthier-Villars, Paris, 1965. \bibitem{b3} G. Boillat; Ondes asymptotiques nonlin\'{e}aires, \emph{Annali di Matematica Pura ed Applicata}, \textbf{91} (IV) (1976), 31-44. \bibitem{c1} Y. Choquet-Bruhat; Ondes asymptotiques et approch\'{e}es pour syst\`emes d'\'{e}quations aux d\'{e}rive\'{e}s partielles nonlin\'{e}aires, \textit{J. Math. Pures et Appl.}, 1968. \bibitem{c2} V. Ciancio, G. A. Kluitenberg; On linear dynamical equations of state for isotropic media - II - Some cases of special interest, \emph{Physica}, \textbf{99 A} (1979), 592-600. \bibitem{c3} V. Ciancio, L. Restuccia; Asymptotic waves in anelastic media without memory (Maxwell media), \emph{Physica}, \textbf{131 A} (1985), 251-262. \bibitem{c4} V. Ciancio, L. Restuccia; Nonlinear dissipative waves in viscoanelastic media, \emph{Physica}, \textbf{132 A} (1985), 606-616. \bibitem{c5} V. Ciancio, L. Restuccia; The generalized Burgers equation in viscoanelastic media with memory, \emph{Physica}, \textbf{142 A}(1987), 309-320. \bibitem{c6} J. D. Cole; \emph{Perturbation methods in applied mathematics}, Blaisdell, Waltham, MA, 1968. \bibitem{d1} A. Donato, A. M. Greco; \emph{ Metodi qualitativi per onde non lineari}, Quaderni del C. N. R., Gruppo Nazionale di Fisica Matematica, 11-th Scuola Estiva di Fisica Matematica, Ravello, 1986, 8-20 September (1986). \bibitem{e1} W. Eckhaus; \emph{Asymptotic analysis of singular perturbation}, North- Holland, Amsterdam, 1979. \bibitem{f1} D. Fusco; Onde non lineari dispersive e dissipative, \emph{Bollettino U.M.I.}, \textbf{5} (16-A) (1979), 450-458. \bibitem{g1} A. Georgescu; \emph{Asymptotic treatment of differential equations}, Chapman and Hall, London, 1995. \bibitem{g2} A. Georgescu, L. Restuccia; Asymptotic waves as solutions of nonlinear PDEs deduced by double-scale method in viscoanelastic media with memory, \emph{Romai Journal}, \textbf{6}(2) (2010), 139-155. \bibitem{g3} A. Georgescu, L. Restuccia; An application of double-scale method to the study of non-linear dissipative waves in Jeffreys media, \emph{Annals of the Academy of Romanian Scientists, Series on Mathematics and its Applications}, \textbf{3}(1) (2011), 116-134. \bibitem{g4} A. Georgescu, L. Restuccia; Asymptotic waves from the point of view of double-scale, \emph{Atti Accademia Peloritana dei Pericolanti} \textbf{LXXXIV} (2006), DOI:10.1478/C1A0601005. \bibitem{h1} E. Hopf; The partial differential equation $u_t+uu_x=\mu u_{xx}$, \emph{Comm. Pure Appl. Math.}, \textbf{3} (1950), 201-230. \bibitem{j1} A. Jeffrey; The development of jump discontinuities in nonlinear hyperbolic systems of equations in two independent variables; \emph{Arch. Rational Mech. Anal.}, \textbf{14} (1963), 27-37. \bibitem{j2} A. Jeffrey; The propagation of weak discontinuities in quasilinear symmetric hyperbolic systems, \emph{Zamp} \textbf{14}: (1963), 301-314. \bibitem{j3} A. Jeffrey; \emph{Quasilinear hyperbolic systems and waves}, Pitman, London, 1976. \bibitem{j4} A. Jeffrey, T. Taniuti; \emph{Nonlinear wave propagation}, Academic Press, Inc., New York, 1964. \bibitem{k1} J. Kevorkian; \emph{M\'{e}thodes des \'{e}chelles multiple}, S\'{e}minaire de l'\'{E}cole Nationale Sup\'{e}rieure de Techniques Avanc\'{e}es, 1972. \bibitem{k2} G. A. Kluitenberg; A thermodynamic derivation of the stress-strain relations for Burgers media and related substancesm, \emph{Physica}, \textbf{38} (1968), 513-548. \bibitem{k3} G. A. Kluitenberg; On the thermodynamics of viscosity and plasticity, \emph{Physica}, \textbf{29} (1963), 633-652. \bibitem{k4} G. A. Kluitenberg; On heat dissipation due to irreversible mechanical phenomena in continuous media, \emph{Physica}, \textbf{35} (1967), 117-192. \bibitem{k5} G. A. Kluitenberg, V. Ciancio; On linear dynamical equations of state for isotropic media - I - General formalism, \emph{Physica}, \textbf{93 A} (1978), 273-286. \bibitem{k6} N. M. Krylov, N. N. Bogoliubov, \emph{Introduction to nonlinear mechanics}, Izd. AN USS, 1937 (Russian) \bibitem{l1} P. A. Lagerstrom, R. G. Casten, Basic concepts underlying singular perturbation technique, \emph{SIAM Rev.}, \textbf{14} (1) (1972), 63-120. \bibitem{l2} P. D. Lax; \emph{Contributions to the theory of partial differential equations}, Princeton University Press, Princeton, 1954. \bibitem{l3} P. D. Lax; Nonlinear hyperbolic equations, \emph{Comm. Pure Appl. Math.}, \textbf{6} (1983), 231-258. \bibitem{m1} Yu. A. Mitropolsk\"u; \emph{Probl$\grave{e}$mes de la th\'{e}orie asymptotique des oscillations non-stationnaires}, Gauthier-Villars, Paris, 1966. \bibitem{o1} R. E. Jr. O'Malley, Topics in singular perturbations, \emph{Adv. Math.}, \textbf{2} (4), (1968) (1968), 365-470. \bibitem{r1} L. Restuccia, A. Georgescu; Determination of asymptotic waves in Maxwell media by double-scale method, \emph{Technische Mechanik}, \textbf{28} (2) (2008), 140-151. \bibitem{r2} L. Restuccia, G. A. Kluitenberg; On the heat dissipation function for irreversible mechanical phenomena in anisotropic media, \emph{Rendiconti del Seminario Matematico di Messina}, \textbf{7} (II , Tomo XXII): (2000), 169-187. \bibitem{s1} D. R. Smith, The multivariable method in singular perturbation analysis, \emph{SIAM Rev.}, \textbf{17} (1975), 221-273. \bibitem{t1} E. Turrisi, V. Ciancio, G. A. Kluitenberg; On the propagation of linear transverse acoustic waves in isotropic media with mechanical relaxation phenomena due to viscosity and a tensorial internal variable II. Some cases of special interest (Poynting-Thomson, Jeffreys, Maxwell, Kelvin-Voigt, Hooke and Newton media), \textit{ Physica}, \textbf{116 A} (1982), 594. \bibitem{v1} G. Veronis; A note on the method of multiple scales, \emph{A. Appl. Math.} \textbf{38} (1980), 363-368. \bibitem{w1} D. J. Wolkind, Singular perturbation techniques, \emph{SIAM Rev.}, \textbf{19} (3) (1977), 502-516. \end{thebibliography} \end{document}