Dispersive homogenized models and coefficient formulas for waves in general periodic media

Abstract

We analyze a homogenization limit for the linear wave equation of second order. The spatial operator is assumed to be of divergence form with an oscillatory coefficient matrix $a^{ε}$ that is periodic with characteristic length scale ε; no spatial symmetry properties are imposed. Classical homogenization theory allows to describe solutions $u^{ε}$ well by a non-dispersive wave equation on fixed time intervals $(0, T)$ . Instead, when larger time intervals are considered, dispersive effects are observed. In this contribution we present a well-posed weakly dispersive equation with homogeneous coefficients such that its solutions $w^{ε}$ describe $u^{ε}$ well on time intervals $(0, T ε^{- 2})$ . More precisely, we provide a norm and uniform error estimates of the form $∥ u^{ε} (t) - w^{ε} (t) ∥ ⩽ C ε$ for $t \in (0, T ε^{- 2})$ . They are accompanied by computable formulas for all coefficients in the effective models. We additionally provide an ε-independent equation of third order that describes dispersion along rays and we present numerical examples.

Keywords

wave equation large time homogenization dispersive model Bloch analysis

1. Introduction

Waves in heterogeneous media exhibit dispersion. This fact is well known in physics, it can be observed in microscopically discrete systems [19,21], but it can also be observed for waves that obey microscopically the classical, non-dispersive wave equation with variable coefficients. Our aim in this contribution is to cast the effect in mathematical terms, to present a well-posed, dispersive effective wave equation, and to provide computable formulas for the (homogeneous) coefficients in the effective equation.

Our analysis concerns solutions $u^{ε} : R^{n} \times (0, \infty) \to R$ , $n \in {1, 2, 3}$ , of the linear wave equation in periodic media, $\partial_{t}^{2} u^{ε} (x, t) = \nabla \cdot (a^{ε} (x) \nabla u^{ε} (x, t)) .$ (1.1) The medium is characterized by a positive, symmetric coefficient matrix field $a^{ε} : R^{n} \to R^{n \times n}$ . We are interested in periodic media with a small periodicity length-scale $ε > 0$ , and assume that $a^{ε} (x) = a_{Y} (x / ε)$ where $a_{Y} : R^{n} \to R^{n \times n}$ is periodic with the periodicity of the unit cell $Y = {(- π, π)}^{n}$ . Except for positivity, matrix symmetry and periodicity, no assumptions on $a_{Y} (\cdot)$ are made (in contrast to our earlier paper [12], where certain spatial symmetries are exploited). Our interest is to describe the solutions $u^{ε}$ for large times, $t \sim ε^{- 2}$ . For classical homogenization results (derivation of effective equations on fixed time intervals) we refer to [6,17] and mention here that, due to energy conservation, even classical homogenization results for the wave equation are much more involved than corresponding results e.g. for the heat equation (the “intermediate case”, the wave equation with damping is considered in [20]). To simplify the exposition, we work here with smooth coefficients $a^{ε}$ . Working with less restrictive assumptions does not lead to major difficulties (in contrast e.g. to observability results, where the regularity of the coefficient seems to be crucial [7]).

In order to have a well-defined object $u^{ε}$ , we must complement the wave equation with an initial condition. For convenience, we restrict our analysis to a vanishing initial velocity, i.e. to initial data $u^{ε} (x, 0) = f (x), \partial_{t} u^{ε} (x, 0) = 0 .$ (1.2) In our mathematical results, we will assume smoothness of f. More precisely, we assume that $f \in L^{2} (R^{n}) \cap L^{1} (R^{n})$ has the Fourier representation $f (x) = \frac{1}{{(2 π)}^{n / 2}} \int_{R^{n}} F_{0} (k) e^{+ i k \cdot x} d k,$ (1.3) where $F_{0} : R^{n} \to C$ has compact support $K \subset R^{n}$ .

We note that our assumptions imply the smoothness $f \in C^{\infty} (R^{n})$ . Less regular initial data can be treated with the help of our results, exploiting the linearity of the equations: Decomposing initial data with bounded energy into two parts, our results can be applied to the smooth part, while the other part generates an error that is, for all times, small in energy norm.

1.1. Known results on dispersive models

The contribution [16] started a series of articles [13–16] which is concerned with the derivation of dispersive models for the wave equation. The authors perform asymptotic (two-scale) expansions of $u^{ε}$ in ε and obtain with their formal calculations a fourth order equation of the form $\partial_{t}^{2} U^{ε} = {A D}^{2} U^{ε} - ε^{2} {C D}^{4} U^{ε},$ (1.4) where A and C are homogeneous coefficients and D denotes spatial derivatives. They call this equation “bad Boussinesq equation”, a well-chosen name, considering the fact that the equation is ill-posed (in the homogenization process, a positive matrix A and a non-positive tensor C appear). We note that in the earlier article [23] this equation also appears (with a sign typo) as a result of a Bloch analysis, but it is not further analyzed in [23]. In [13–16] various approaches for a further (analytical and numerical) exploitation of Eq. (1.4) are investigated: regularizations, non-local approximations and multiple time scales.

The first rigorous result that establishes a dispersive model for the wave Eq. (1.1) appeared in [18]. In that work, which is concerned with the one-dimensional case, the well-posed dispersive Eq. (1.10) is formulated and an error estimate similar to (2.1) is derived. The method of proof is very different from our approach here (which is as in [12]): Adaption operators are constructed and used to adapt smooth solutions of the homogeneous dispersive system to the periodic medium. After the adaption, direct energy procedures can be applied.

Another mathematical derivation of dispersive limits is performed in [4,5]. The wave equation is scaled as in our setting (time scales of order $ε^{- 1}$ are investigated), but the initial data are assumed to be oscillatory at scale ε and are described by Bloch wave packets. In this setting, the effective diffraction can be described by a Schrödinger equation for the envelope function. Another scaling of the system is analyzed in [2], where large potentials instead of large time spans are considered.

1.1.1. Bloch analysis

The central tool in a Bloch analysis is the Bloch expansion of an arbitrary function (in our case the solution $u^{ε}$ ). While in a Fourier expansion one uses the dual variable $k \in R^{n}$ , the Bloch expansion uses two dual variables, k and m. Since $m \in N_{0}$ is an additional parameter, the other parameter varies only in a restricted domain, the Brillouin zone, $k \in Z : = {(- 1 / 2, 1 / 2)}^{n}$ . The basis functions $e^{i k \cdot x}$ of the Fourier analysis are replaced by solutions $ψ_{m} (\cdot, k)$ of the Bloch eigenvalue problem $- (\nabla_{y} + i k) \cdot (a_{Y} (y) (\nabla_{y} + i k) ψ_{m} (y, k)) = μ_{m} (k) ψ_{m} (y, k) .$ (1.5) Here $ψ_{m} (\cdot, k) : Y \to C$ is a periodic function, $ψ_{m} (\cdot, k) \in H_{per}^{1} (Y)$ , $0 ⩽ μ_{0} (k) ⩽ μ_{1} (k) ⩽ \dots$ are the ordered, real eigenvalues. Equation (1.5) is constructed in such a way that the function $w_{m} (y, k) : = ψ_{m} (y, k) e^{i k \cdot y}$ is an eigenfunction to the elliptic operator with periodic coefficients $a_{Y}$ .

Bloch wave homogenization theory establishes that the effective behavior of $u^{ε}$ in the limit $ε \to 0$ is characterized solely by the behavior of the smallest eigenvalue $μ_{0} (k)$ in a neighborhood of $k = 0 \in Z$ . For such results in classical homogenization settings, we refer to [3,8–11]. In Fig. 1 we plot the Bloch wave functions $ψ_{0} (y, k)$ for $k = (1 / 2, 0)$ and $k = (1 / 2, 1 / 2)$ . For an illustration of the eigenvalue structure see Fig. 2(a).

Fig. 1.

The Bloch wave $ψ_{0} (y, k)$ at $k = (1 / 2, 0)$ in (a) and (b) and at $k = (1 / 2, 1 / 2)$ in (c) and (d) corresponding to $μ_{0}$ for $a_{Y}$ from (4.8). (a) $Re (ψ_{0} (y, k))$ . (b) $Im (ψ_{0} (y, k))$ . (c) $Re (ψ_{0} (y, k))$ . (d) $Im (ψ_{0} (y, k))$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

Fig. 2.

(a) The first three eigenvalues $μ_{0} (k)$ , $μ_{1} (k)$ , $μ_{2} (k)$ of (1.5) for $a_{Y}$ from (4.8). (b) Experimental convergence in ε of the error $e_{ε} : = {∥ (u^{ε} - w^{ε}) (t = ε^{- 2} / 2) ∥}_{L^{2} (Ω)}$ . A logarithmic scale is used on both axes. The line with slope 1.1 was obtained by a linear interpolation of the error in the logarithmic scale for the 5 smallest ε-values. The experimental rate is thus $e_{ε} \sim ε^{1.006}$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

The Bloch wave homogenization method was used for an analysis of higher order effects of the heterogeneity of the medium in the influential article [23]. That article does not formulate a well-posed dispersive effective equation (hence, in particular, it does not provide an error estimate), but it gives a lot of insight into the dispersive limit: Even the effective long time behavior of $u^{ε}$ is characterized by $μ_{0} (k)$ and its behavior near $k = 0$ . Expanding $μ_{0} (k)$ in a Taylor series around $k = 0$ , we may write $μ_{0} (k) = \sum A_{l m} k_{l} k_{m} + \sum C_{l m n q} k_{l} k_{m} k_{n} k_{q} + O ({| k |}^{6}),$ (1.6) where odd derivatives vanish due to the symmetry $μ_{0} (k) = μ_{0} (- k)$ , sums are over repeated indices. The matrix A and the tensor C provide the coefficients in the formal Eq. (1.4), where ${C D}^{4}$ is the spatial fourth order operator ${C D}^{4} = \sum C_{i j k l} \partial_{i} \partial_{j} \partial_{k} \partial_{l} .$ (1.7) While A is positive definite and symmetric, C turns out to be negative semi-definite, a fact that is shown in [10]. As a consequence, the differential operator ${A D}^{2}$ is negative and the operator $- ε^{2} {C D}^{4}$ is non-negative. For this reason, Eq. (1.4) is ill-posed.

Let us be more precise about the arguments in the Bloch analysis. The key idea is to perform a Bloch expansion of the solution $u^{ε}$ . We use $L^{2} (Y)$ -normalized Bloch eigenfunctions $ψ_{m}$ solving (1.5) and define scaled functions $ψ_{m}^{ε} (x, k) : = ψ_{m} (x / ε, ε k)$ . The initial values f are expanded in terms of these basis functions with Bloch coefficients ${\hat{f}}_{m}^{ε} (k)$ . The wave equation demands that second time derivatives have the same effect on $u^{ε}$ as the application of the elliptic operator with periodic coefficients. Formally, this leads to the following expression of the solution: $u^{ε} (x, t) = \sum_{m = 0}^{\infty} \int_{Z / ε} {\hat{f}}_{m}^{ε} (k) ψ_{m} (x / ε, ε k) e^{i k \cdot x} Re (e^{i t \sqrt{μ_{m} (ε k)} / ε}) d k .$ (1.8) In fact, the expansion (1.8) can be rigorously justified, see Lemma 2.1 of [12]. In the next steps, this formula is simplified for small $ε > 0$ : One realizes, to leading order in ε, that only $m = 0$ has to be considered, that ${\hat{f}}_{m}^{ε}$ can be replaced by the Fourier transform $F_{0}$ of the initial values, and that $ψ_{0}$ can be replaced by the constant ${(2 π)}^{- n / 2}$ . Expanding finally $μ_{0} (ε k)$ in ε, one finds the following expression, which can be used to define an approximate solution $v^{ε}$ , $\begin{array}{rclr} v^{ε} (x, t) & : = & {(2 π)}^{- n / 2} \frac{1}{2} \sum_{\pm} \int_{K} F_{0} (k) e^{i k \cdot x} exp (\pm i t \sqrt{\sum A_{l m} k_{l} k_{m}}) \\ \times exp (\pm \frac{i ε^{2}}{2} t \frac{\sum C_{l m n q} k_{l} k_{m} k_{n} k_{q}}{\sqrt{\sum A_{l m} k_{l} k_{m}}}) d k . & (1.9) \end{array}$ Equation (1.4) with tensors A and C is constructed in such a way that (formally) the function $v^{ε}$ is a solution up to errors of order $ε^{4}$ .

1.1.2. Rigorous approximation results

A rigorous mathematical analysis can be performed when Eq. (1.4) is transformed into a well-posed equation, using the replacement ${A D}^{2} U^{ε} \approx \partial_{t}^{2} U^{ε}$ to re-write the operator ${C D}^{4}$ . The first utilization of this trick for a rigorous result seems to be in the treatment [18] in the one-dimensional case. More recently, we were able to exploit the same trick in arbitrary dimension in [12] (under quite strong spatial symmetry assumptions on the coefficient field $a_{Y} (\cdot)$ ).

In that contribution, we use a rigorous Bloch wave analysis to show that $v^{ε}$ of (1.9) approximates $u^{ε}$ in appropriate norms. In a second step, we show with energy methods that $v^{ε}$ is close to the solution $w^{ε}$ of the well-posed, weakly dispersive equation $\partial_{t}^{2} w^{ε} = {A D}^{2} w^{ε} + ε^{2} {E D}^{2} \partial_{t}^{2} w^{ε} - ε^{2} {F D}^{4} w^{ε} .$ (1.10) The positive semi-definite and symmetric tensors E and F are constructed in such a way that $- {C D}^{4} = {E D}^{2} {A D}^{2} - {F D}^{4} .$ (1.11) Together, the two estimates provide an estimate for $u^{ε} - w^{ε}$ . This shows that the weakly dispersive Eq. (1.10) is a valid replacement for the original Eq. (1.1) on large time intervals.

1.2. New results

In the article at hand, we obtain the long-time homogenization result for very general coefficient fields $a_{Y} (\cdot)$ . In particular, we show that the well-posed Eq. (1.10) provides the effective description of solutions for large times, characterizing dispersion in arbitrary dimension and without spatial symmetry assumptions. Furthermore, we provide explicit formulas for the effective coefficients.

Decomposition lemma and approximation result

In order to show the approximation result, we can rely on the Bloch wave analysis of [12]. The only new ingredient is a considerably developed decomposition lemma: Lemma 2.5 yields that, without any structural assumptions on C, the differential operator ${C D}^{4}$ can be written as in (1.11) for appropriate semi-definite and symmetric tensors E and F (using the given positive, symmetric matrix A). In particular, the lemma allows to decompose the operator ${C D}^{4}$ also when the coefficients $a_{Y} (\cdot)$ have no spatial symmetries.

Once the decomposition lemma is established, we can apply results of [12]. We obtain that (1.10) is well posed and that solutions $w^{ε}$ approximate the solutions $u^{ε}$ . More precisely, from the analysis of [12], we obtain an error estimate of the form $∥ u^{ε} (t) - w^{ε} (t) ∥ ⩽ C_{0} ε$ uniformly in $t \in (0, T ε^{- 2})$ in general periodic media. This approximation result is stated and proved in Section 2.

An ε-independent third order dispersive equation

Our aim in Section 3 is to provide a simplified model in which no ε-dependence occurs. The approximation result of Theorem 2.2 allows to analyze, instead of the solution $u^{ε}$ of the original problem, the solution $w^{ε}$ of the weakly dispersive Eq. (1.10) . In Section 3 we analyze $w^{ε}$ in two dimensions in polar coordinates. On every ray through the origin, for an appropriate scaling of the solution, we can perform the limit $ε \to 0$ . The result is an equation that determines the shape of pulses, namely a linear third order equation (a linearized KdV-equation). All coefficients in this equation are computable from the coefficient field $a_{Y}$ .

Algorithms to compute effective quantities

In Section 4 we present an algorithm that provides the homogeneous coefficients in all effective equations, i.e. A and C (in a more direct form than as derivatives of the Bloch eigenvalue), the coefficients E and F, and the coefficients of the linearized KdV-equation.

The numerical results of Section 4 compare solutions to the original problem with solutions to the effective problems. We find a remarkable qualitative and quantitative agreement also for moderate ε, and the correct experimental convergence rates for $ε \to 0$ .

2. A weakly dispersive effective equation

In Section 2.1 we formulate our main result. The main part of its proof can be obtained by applying results of [12]; this is the subject of Section 2.2. The new ingredient is the decomposition lemma, which is shown in Section 2.3.

2.1. Main approximation result

We emphasize that the set $Y \subset R^{n}$ , the reciprocal cell $Z : = {(- 1 / 2, 1 / 2)}^{n} \subset R^{n}$ , and the support $K \subset R^{n}$ are fixed data of the problem. Given is also the coefficient field $a_{Y}$ that determines $a^{ε} (x) = a_{Y} (x / ε)$ .

Assumption 2.1.
The coefficient field $a_{Y} : R^{n} \to R^{n \times n}$ is Y-periodic for the cube $Y : = {(- π, π)}^{n} \subset R^{n}$ and has the regularity $a_{Y} \in C^{1} (R^{n}, R^{n \times n})$ . The matrix $a_{Y} (y)$ is symmetric for every point $y \in R^{n}$ , i.e. ${(a_{Y} (y))}_{i j} = {(a_{Y} (y))}_{j i}$ for all $i, j \in {1, \dots, n}$ . The field is positive definite: for some $γ > 0$ there holds $\sum_{i, j = 1}^{n} {(a_{Y} (y))}_{i j} ξ_{i} ξ_{j} ⩾ γ {| ξ |}^{2}$ for every $y \in R^{n}$ and every $ξ \in R^{n}$ .

Our main approximation result is stated in the following theorem. The result is very similar to the main theorem in [12]. The difference is that we do not assume any spatial symmetry of $a_{Y} (\cdot)$ (such as a reflection symmetry in each coordinate direction or symmetry with respect to exchanging coordinate axes). We note that dimensions $n > 3$ can be treated by assuming higher regularity properties.
Theorem 2.2. (Approximation)

Let $ε = ε_{l} \to 0$ be a sequence of positive numbers and $n \in {1, 2, 3}$ be the dimension. Let the medium $a_{Y} : R^{n} \to R^{n \times n}$ satisfy Assumption 2.1 and let the initial data $f : R^{n} \to R$ be as in (1.3). We use the coefficient matrices A and Cdefined in (1.6) . Let E and Fbe positive semi-definite such that (1.11) holds (the existence is established in Lemma 2.5). Then the following holds:

Well-posedness. Equation (1.10) with initial condition (1.2) has a unique solution $w^{ε}$ for all positive times (see Theorem 2.4 for function spaces).

Error estimate. Let $w^{ε}$ be the solution of (1.10), and let $u^{ε}$ be the solution of (1.1) for the same initial condition (1.2). Then, with a constant $C_{0} = C_{0} (a_{Y}, T_{0}, f)$ , there holds the error estimate $sup_{t \in [0, T_{0} ε^{- 2}]} {∥ u^{ε} (\cdot, t) - w^{ε} (\cdot, t) ∥}_{L^{2} (R^{n}) + L^{\infty} (R^{n})} ⩽ C_{0} ε .$ (2.1)

Here we use, for two Banach spaces with norms ${∥ \cdot ∥}_{X}$ and ${∥ \cdot ∥}_{Y}$ , the weaker norm ${∥ u ∥}_{X + Y} : = inf {{∥ u_{1} ∥}_{X} + {∥ u_{2} ∥}_{Y} : u = u_{1} + u_{2}}$ . Such a norm appears in our main Theorem 2.2, since different contributions to the error $u^{ε} - w^{ε}$ are measured in different norms.

2.2. Proof of Theorem 2.2

The following corollary is the central result of the Bloch analysis. It is derived with mathematical rigor in [12]; it provides a comparison between the solution $u^{ε}$ of the heterogeneous wave equation with the explicitly defined function $v^{ε}$ .

Theorem 2.3. (Corollary 2.5 of [12])

Let Assumption 2.1 be satisfied. Let $u^{ε}$ be the solution of (1.1) and let $v^{ε}$ be defined by (1.9) . Then $sup_{t \in [0, T_{0} ε^{- 2}]} {∥ u^{ε} (\cdot, t) - v^{ε} (\cdot, t) ∥}_{L^{2} (R^{n}) + L^{\infty} (R^{n})} ⩽ C_{0} ε .$ (2.2)

The following theorem is based on energy methods. It provides the comparison between the solution $w^{ε}$ of the weakly dispersive (homogeneous) equation and the explicit function $v^{ε}$ .

Theorem 2.4. (Theorem 3.3 of [12])

Let A, C, E, F be tensors with the properties: $A \in R^{n \times n}$ is symmetric and positive definite, $\sum_{i j} A_{i j} ξ_{i} ξ_{j} ⩾ γ {| ξ |}^{2}$ for some $γ > 0$ , $E \in R^{n \times n}$ and $F \in R^{n \times n \times n \times n}$ are positive semi-definite and symmetric, $C \in R^{n \times n \times n \times n}$ allows the decomposition (1.11). Then the following holds.

( 1 ) Well-posedness. For the initial datum $f \in H^{2} (R^{n})$ , the equation $\begin{matrix} \partial_{t}^{2} w^{ε} - {A D}^{2} w^{ε} - ε^{2} \partial_{t}^{2} {E D}^{2} w^{ε} + ε^{2} {F D}^{4} w^{ε} = 0, \\ w^{ε} (\cdot, 0) = f, \partial_{t} w^{ε} (\cdot, 0) = 0 \end{matrix}$ (2.3) has a unique solution $w^{ε} \in L^{\infty} (0, T_{0} ε^{- 2}; H^{2} (R^{n})) \cap W^{1, \infty} (0, T_{0} ε^{- 2}; H^{1} (R^{n}))$ .

( 2 ) Approximation. Let $v^{ε}$ be defined by (1.9) where $F_{0}$ and fare related by (1.3) . Let $w^{ε}$ be a solution of (2.3) . Then there holds $sup_{t \in [0, T_{0} ε^{- 2}]} {∥ \partial_{t} (v^{ε} - w^{ε}) (\cdot, t) ∥}_{L^{2} (R^{n})} + sup_{t \in [0, T_{0} ε^{- 2}]} {∥ \nabla (v^{ε} - w^{ε}) (\cdot, t) ∥}_{L^{2} (R^{n})} ⩽ C_{0} ε^{2}$ (2.4) with a constant $C_{0}$ that is independent of ε.

Theorem 2.2 is a consequence of the above results and a new decomposition lemma that is shown in the next subsection.

The estimate (2.2) provides that $∥ u^{ε} - v^{ε} ∥$ is of order ε. The norms coincide with the ones in the claim (2.1). It therefore remains to estimate the difference $∥ v^{ε} - w^{ε} ∥$ .

We define A and C through (1.6) and note that A is positive definite and symmetric. The decomposition result of Lemma 2.5 allows to construct E and F such that (1.11) is satisfied. This means that Theorem 2.4 can be applied. It provides the well-posedness claim and the estimate (2.4), which shows that norms of derivatives of $v^{ε} - w^{ε}$ are of order $ε^{2}$ . The norms can be transformed with an interpolation lemma (see Lemma 3.4 in [12]); the result is an estimate for $sup_{t \in [0, T_{0} ε^{- 2}]} {∥ v^{ε} (\cdot, t) - w^{ε} (\cdot, t) ∥}_{L^{2} + L^{\infty}}$ of order ε. We therefore obtain (2.1).

2.3. Decomposition lemma

Our aim is to construct coefficient tensors $E \in R^{n \times n}$ and $F \in R^{n \times n \times n \times n}$ such that the differential operator ${C D}^{4}$ can be re-written as in (1.11) . We impose that E and F are positive semi-definite and symmetric, i.e. $\sum_{i, j, k, l = 1}^{n} F_{i j k l} ξ_{i j} ξ_{k l} ⩾ 0 for every ξ \in R^{n \times n}, F_{i j k l} = F_{k l i j},$ (2.5) and similarly $\sum_{i, j = 1}^{n} E_{i j} η_{i} η_{j} ⩾ 0$ for every $η \in R^{n}$ and $E_{i j} = E_{j i}$ . The symmetry relations must hold for all indices $i, j, k, l \in {1, \dots, n}$ . In this section, $n \in {1, 2, 3, 4, \dots}$ can also be larger than 3.

Lemma 2.5. (Decomposability)

Let $A \in R^{n \times n}$ be a symmetric and positive definite matrix and let $C \in R^{n \times n \times n \times n}$ be arbitrary. Then there exist symmetric and positive semi-definite tensors $E \in R^{n \times n}$ and $F \in R^{n \times n \times n \times n}$ such that the differential operator ${C D}^{4}$ can be written as in (1.11) .

We remark that the tensors E and F in the above lemma are not uniquely determined. On the contrary, having found one pair of tensors, it is easy to construct a large family of other admissible tensors as follows: Let E and F be admissible and let $\hat{E} \in R^{n \times n}$ be an arbitrary symmetric, positive semidefinite tensor. We construct new tensors $\tilde{E}$ , $\tilde{F}$ by setting $\tilde{E} : = E + \hat{E}$ and $\tilde{F} : = F + \hat{E} \otimes A$ . A simple calculation shows that $\tilde{E}$ , $\tilde{F}$ are admissible in the sense of Lemma 2.5.

Concerning the non-uniqueness of E and F we further remark that the corresponding different weakly dispersive equations differ only up to terms of order $O (ε^{4})$ . Such errors are not relevant on time scales of order $O (ε^{- 2})$ , hence our approximation theorem holds for every admissible choice of E and F. On the other hand, we may ask the following: Is there a natural choice for E and F? This might be the case in the sense of a variational principle or in the sense of error minimization.

The proof of the decomposition Lemma 2.5 consists of two steps. In the first (and essential) step we show that the decomposability result holds for diagonal matrices $A \in R^{n \times n}$ and with rotated derivative operators. In the second step, general matrices A are treated by diagonalization.

Lemma 2.6. (Decomposability for diagonal matrices A)

Let Abe a positive definite diagonal matrix, $A = diag (a_{1}, a_{2}, \dots, a_{n})$ . Let $S \in SO (n)$ be an orthogonal matrix and let $C \in R^{n \times n \times n \times n}$ be arbitrary. Then there exist symmetric and positive semi-definite tensors $E \in R^{n \times n}$ and $F \in R^{n \times n \times n \times n}$ such that $- C {\tilde{D}}^{4} = E {\tilde{D}}^{2} A {\tilde{D}}^{2} - F {\tilde{D}}^{4} .$ (2.6) Here, the operator $\tilde{D} : = S D$ denotes a rotated derivative.

Proof.

Step 1: Reduction to matrices C with only one non-trivial entry. We note that the relation (2.6) is additive in the following sense: Let $C^{(1)}$ and $C^{(2)}$ be two matrices, and let (2.6) be satisfied for $C^{(m)}$ with tensors $E^{(m)}$ and $F^{(m)}$ , $m \in {1, 2}$ . Then (2.6) holds for $C^{(1)} + C^{(2)}$ with the two tensors $E^{(1)} + E^{(2)}$ and $F^{(1)} + F^{(2)}$ . We exploit here that the sum of symmetric, semi-definite tensors is again symmetric and semi-definite.

This observation implies that it is sufficient to consider a tensor C that has only one non-trivial entry. We denote a tensor with only one entry 1 in the canonical way as $e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ . An arbitrary tensor C can be written as a sum, $C = \sum_{α β γ δ} C_{α β γ δ} e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ . After constructing tensors $E^{(α, β, γ, δ)}$ and $F^{(α, β, γ, δ)}$ according to the one-entry tensor $C_{α β γ δ} e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ , we find E and F according to C by a summation, $E = \sum_{α β γ δ} E^{(α, β, γ, δ)}$ and $F = \sum_{α β γ δ} F^{(α, β, γ, δ)}$ .

In the following construction, we restrict ourselves to a fixed choice of indices, $(α, β, γ, δ) \in {1, \dots, n}^{4}$ . For a number $c \in R$ , we can consider a tensor C of the form $C = c e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ .

Our aim is to re-write the differential operator $- C {\tilde{D}}^{4} = c {\tilde{D}}_{α} {\tilde{D}}_{β} {\tilde{D}}_{γ} {\tilde{D}}_{δ}$ . We use here ${\tilde{D}}_{i} : = \sum_{j = 1}^{n} S_{i j} \partial_{x_{j}}$ for the ith component of the rotated gradient. In the following, ${a}_{+} : = max {0, a}$ denotes the positive part of a number $a \in R$ .

Step 2: Construction of E and Ffor $C = c e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ , where at least two indices coincide. Case 1.

The indices α, β, γ, δ contain two pairs, i.e. $(α, β, γ, δ) = (i, i, j, j)$ or $(α, β, γ, δ) = (i, j, i, j)$ or $(α, β, γ, δ) = (i, j, j, i)$ for $i, j \in {1, \dots, n}$ (four equal indices are allowed). We restrict ourselves to $(α, β, γ, δ) = (i, i, j, j)$ , the permutations define the same operator $C {\tilde{D}}^{4}$ . We define the tensors $E = E^{(α, β, γ, δ)}$ and $F = F^{(α, β, γ, δ)}$ through $E_{i i} : = \frac{{- c}_{+}}{a_{j}}, F_{i j i j} : = {c}_{+}, F_{i m i m} : = \frac{{- c}_{+}}{a_{j}} a_{m}$ for all $m \in {1, \dots, n}$ with $m \neq j$ . All other entries of E and F are set to zero.

Properties of E and F. By definition E and F are symmetric and positive semi-definite. A direct calculation yields the decomposition property: $\begin{array}{rcl} E {\tilde{D}}^{2} A {\tilde{D}}^{2} - F {\tilde{D}}^{4} \\ = (\frac{{- c}_{+}}{a_{j}} {\tilde{D}}_{i}^{2}) (\sum_{m} a_{m} {\tilde{D}}_{m}^{2}) - {c}_{+} {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} - \sum_{m \neq j} \frac{{- c}_{+}}{a_{j}} a_{m} {\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} \\ = {- c}_{+} {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} + \sum_{m \neq j} \frac{{- c}_{+}}{a_{j}} a_{m} {\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} - {c}_{+} {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} - \sum_{m \neq j} \frac{{- c}_{+}}{a_{j}} a_{m} {\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} \\ = ({- c}_{+} - {c}_{+}) {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} = - c {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} = - C {\tilde{D}}^{4} . \end{array}$ Case 2.

The indices α, β, γ, δ contain three identical entries, i.e. $(α, β, γ, δ) = (i, i, i, j)$ or $(α, β, γ, δ) = (i, i, j, i)$ or $(α, β, γ, δ) = (i, j, i, i)$ or $(α, β, γ, δ) = (j, i, i, i)$ for $i, j \in {1, \dots, n}$ with $i \neq j$ . We restrict ourselves to $(α, β, γ, δ) = (i, i, i, j)$ since the other cases define the same differential operator. We define the tensors $E = E^{(α, β, γ, δ)}$ and $F = F^{(α, β, γ, δ)}$ through $\begin{array}{rcl} E_{i j} : = E_{j i} : = - \frac{c}{2 a_{i}} = : \tilde{c}, E_{i i} : = E_{j j} : = | \tilde{c} |, \\ F_{i m i m} : = F_{j m j m} : = | \tilde{c} | a_{m} for all m \in {1, \dots, n}, \\ F_{i m j m} : = F_{j m i m} : = \tilde{c} a_{m} for all m \in {1, \dots, n} with m \neq i . \end{array}$ All other entries of E and F are set to zero.

Properties of E and F. By definition, E and F are symmetric. Concerning the positive semi-definiteness of E and F we calculate, for arbitrary $ξ \in R^{n}$ and $ζ \in R^{n \times n}$ , $\begin{array}{rcl} \sum_{l, m} E_{l m} ξ_{l} ξ_{m} = | \tilde{c} | ξ_{i}^{2} + | \tilde{c} | ξ_{j}^{2} + 2 \tilde{c} ξ_{i} ξ_{j} ⩾ | \tilde{c} | ξ_{i}^{2} + | \tilde{c} | ξ_{j}^{2} - (| \tilde{c} | ξ_{i}^{2} + | \tilde{c} | ξ_{j}^{2}) = 0, \\ \begin{matrix} \sum_{l, m, p, q} F_{l m p q} ζ_{l m} ζ_{p q} & = \sum_{m} | \tilde{c} | a_{m} ζ_{i m}^{2} + \sum_{m} | \tilde{c} | a_{m} ζ_{j m}^{2} + \sum_{m \neq i} 2 \tilde{c} a_{m} ζ_{i m} ζ_{j m} \\ ⩾ \sum_{m} | \tilde{c} | a_{m} ζ_{i m}^{2} + \sum_{m} | \tilde{c} | a_{m} ζ_{j m}^{2} - \sum_{m \neq i} (| \tilde{c} | a_{m} ζ_{i m}^{2} + | \tilde{c} | a_{m} ζ_{j m}^{2}) ⩾ 0 . \end{matrix} \end{array}$ Concerning the decomposition property we calculate $\begin{array}{rcl} E {\tilde{D}}^{2} A {\tilde{D}}^{2} - F {\tilde{D}}^{4} \\ = (| \tilde{c} | {\tilde{D}}_{i}^{2} + | \tilde{c} | {\tilde{D}}_{j}^{2} - \frac{c}{a_{i}} {\tilde{D}}_{i} {\tilde{D}}_{j}) (\sum_{m} a_{m} {\tilde{D}}_{m}^{2}) \\ - \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) + \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} \\ = \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) - c {\tilde{D}}_{i}^{3} {\tilde{D}}_{j} - \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} \\ - \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{i}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) + \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} = - c {\tilde{D}}_{i}^{3} {\tilde{D}}_{j} = - C {\tilde{D}}^{4} . \end{array}$ Case 3.

The indices α, β, γ, δ contain two identical entries, the other entries are different. We restrict ourselves to the case $(α, β, γ, δ) = (i, i, j, k)$ with three different indices $i, j, k \in {1, \dots, n}$ , permutations of the indices define the same operator. We define the tensors $E = E^{(α, β, γ, δ)}$ and $F = F^{(α, β, γ, δ)}$ through $\begin{array}{rcl} E_{j k} : = E_{k j} : = - \frac{c}{2 a_{i}} = : \tilde{c}, E_{k k} : = E_{j j} : = | \tilde{c} |, \\ F_{k m k m} : = F_{j m j m} : = | \tilde{c} | a_{m} for all m \in {1, \dots, n}, \\ F_{k m j m} : = F_{j m k m} : = \tilde{c} a_{m} for all m \in {1, \dots, n} with m \neq i . \end{array}$ All other entries of E and F are set to zero.

Properties of E and F. By definition, E and F are symmetric. Concerning the positive semi-definiteness of E and F, we calculate for arbitrary $ξ \in R^{n}$ and $ζ \in R^{n \times n}$ $\begin{array}{rcl} \sum_{l, m} E_{l m} ξ_{l} ξ_{m} = | \tilde{c} | ξ_{k}^{2} + | \tilde{c} | ξ_{j}^{2} + 2 \tilde{c} ξ_{k} ξ_{j} ⩾ 0, \\ \sum_{l, m, p, q} F_{l m p q} ζ_{l m} ζ_{p q} = \sum_{m} | \tilde{c} | a_{m} ζ_{k m}^{2} + \sum_{m} | \tilde{c} | a_{m} ζ_{j m}^{2} + \sum_{m \neq i} 2 \tilde{c} a_{m} ζ_{k m} ζ_{j m} ⩾ 0 . \end{array}$ Regarding the decomposition property we calculate $\begin{array}{rcl} E {\tilde{D}}^{2} A {\tilde{D}}^{2} - F {\tilde{D}}^{4} \\ = (| \tilde{c} | {\tilde{D}}_{k}^{2} + | \tilde{c} | {\tilde{D}}_{j}^{2} - \frac{c}{a_{i}} {\tilde{D}}_{k} {\tilde{D}}_{j}) (\sum_{m} a_{m} {\tilde{D}}_{m}^{2}) \\ - \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{k}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) + \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{k} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} \\ = \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{k}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) - c {\tilde{D}}_{k} {\tilde{D}}_{j} {\tilde{D}}_{i}^{2} - \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{k} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} \\ - \sum_{m} | \tilde{c} | a_{m} ({\tilde{D}}_{k}^{2} {\tilde{D}}_{m}^{2} + {\tilde{D}}_{j}^{2} {\tilde{D}}_{m}^{2}) + \sum_{m \neq i} a_{m} \frac{c}{a_{i}} {\tilde{D}}_{k} {\tilde{D}}_{j} {\tilde{D}}_{m}^{2} = - c {\tilde{D}}_{k} {\tilde{D}}_{j} {\tilde{D}}_{i}^{2} = - C {\tilde{D}}^{4} . \end{array}$

Step 3: Construction of E and Ffor $C = c e_{α} \otimes e_{β} \otimes e_{γ} \otimes e_{δ}$ , where no two indices coincide. We treat now the case $(α, β, γ, δ) = (i, j, k, l)$ with pairwise different indices $i, j, k, l \in {1, \dots, n}$ ; our aim is to rewrite the operator $C {\tilde{D}}^{4} = c {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{k} {\tilde{D}}_{l}$ . We obtain corresponding matrices in two steps: (a) We define a positive semi-definite, symmetric tensor $\hat{F} \in R^{n \times n \times n \times n}$ in such a way that $\bar{C} {\tilde{D}}^{4} : = C {\tilde{D}}^{4} - \hat{F} {\tilde{D}}^{4}$ contains only partial derivatives with respect to repeated indices. (b) We apply Step 2 of this proof to the remainder $\bar{C} {\tilde{D}}^{4}$ , which provides $\bar{E}$ and $\bar{F}$ with $- \bar{C} {\tilde{D}}^{4} = \bar{E} {\tilde{D}}^{2} A {\tilde{D}}^{2} - \bar{F} {\tilde{D}}^{4}$ . The desired tensors $E = E^{(α, β, γ, δ)}$ and $F = F^{(α, β, γ, δ)}$ are then obtained as $E : = \bar{E}$ and $F : = \bar{F} + \hat{F}$ .

We set ${\hat{F}}_{i j k l} : = {\hat{F}}_{k l i j} : = \frac{1}{2} c, {\hat{F}}_{i j i j} : = {\hat{F}}_{k l k l} : = \frac{1}{2} | c | .$ (2.7) All other entries of $\hat{F}$ are set to zero. The symmetry of $\hat{F}$ is obvious; regarding positivity we calculate $\sum_{m, p, q, r} {\hat{F}}_{m p q r} ζ_{m p} ζ_{q r} = \frac{1}{2} | c | ζ_{i j}^{2} + \frac{1}{2} | c | ζ_{k l}^{2} + c ζ_{i j} ζ_{k l} ⩾ 0 .$

It remains to check that Step 2 of this proof can be applied to the remainder $\bar{C} {\tilde{D}}^{4} : = C {\tilde{D}}^{4} - \hat{F} {\tilde{D}}^{4}$ . We evaluate $\begin{array}{rclr} C {\tilde{D}}^{4} - \hat{F} {\tilde{D}}^{4} & = & c {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{k} {\tilde{D}}_{l} - [c {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{k} {\tilde{D}}_{l} + \frac{1}{2} | c | {\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} + \frac{1}{2} | c | {\tilde{D}}_{k}^{2} {\tilde{D}}_{l}^{2}] \\ = & - \frac{1}{2} | c | ({\tilde{D}}_{i}^{2} {\tilde{D}}_{j}^{2} + {\tilde{D}}_{k}^{2} {\tilde{D}}_{l}^{2}) . & (2.8) \end{array}$ This operator has nontrivial entries only for repeated indices, it therefore possesses a decomposition by Step 2. This concludes the proof of the lemma. □

With Lemma 2.6 at hand we are in the position to prove the general decomposition result of Lemma 2.5. Proof of Lemma 2.5.

The symmetry of A implies that A is diagonalizable: there exists an orthogonal matrix $S \in SO (n)$ such that $A = S^{T} \tilde{A} S with \tilde{A} = diag (a_{1}, a_{2}, \dots, a_{n}) .$ Since A is positive definite, the eigenvalues $a_{i}$ are positive. Our aim is to apply Lemma 2.6 with the diagonal matrix $\tilde{A}$ and a tensor $\tilde{C}$ , which is defined from C with the transformation S. Lemma 2.6 provides tensors $\tilde{E}$ and $\tilde{F}$ that can be transformed back into the desired tensors E and F.

Step 1: Construction of $\tilde{C}$ . Here, we denote the space of matrices by $M : = R^{n \times n}$ . The tensor C defines a linear map $C : M \to M$ through ${(C (e_{i} \otimes e_{j}))}_{k l} = C_{i j k l}$ . We define a transformed tensor $\tilde{C} : M \to M$ through $\tilde{C} : M ∋ B \mapsto S \cdot C (S^{T} B S) \cdot S^{T} \in M .$ (2.9) The transformation can be inverted: Given $\tilde{C}$ , the corresponding tensor C is given by $C (b) = S^{T} \tilde{C} (S b S^{T}) S$ .

As we show next, with the transformed differential operator $\tilde{D} : = S D$ , the two fourth order differential operators coincide: ${C D}^{4} = \tilde{C} {\tilde{D}}^{4}$ . We use the convention that sums are over repeated indices. $\begin{array}{rclr} \tilde{C} {\tilde{D}}^{4} & = & \sum {(\tilde{C} (e_{i} \otimes e_{j}))}_{k l} {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{k} {\tilde{D}}_{l} \\ = & \sum S_{k α} C {(S^{T} \cdot e_{i} \otimes e_{j} \cdot S)}_{α β} S_{l β} {\tilde{D}}_{i} {\tilde{D}}_{j} {\tilde{D}}_{k} {\tilde{D}}_{l} \\ = & \sum S_{k α} S_{i γ} S_{j δ} C {(e_{γ} \otimes e_{δ})}_{α β} S_{l β} S_{i ξ} D_{ξ} S_{j ζ} D_{ζ} S_{k η} D_{η} S_{l θ} D_{θ} \\ = & \sum C {(e_{γ} \otimes e_{δ})}_{α β} D_{α} D_{β} D_{γ} D_{δ} = {C D}^{4} . & (2.10) \end{array}$

Step 2: Application of Lemma 2.6. Since $\tilde{A}$ is diagonal, we can apply Lemma 2.6 with the data $\tilde{A}$ , S and $\tilde{C}$ . We find symmetric, positive semi-definite tensors $\tilde{E} \in R^{n \times n}$ and $\tilde{F} \in R^{n \times n \times n \times n}$ such that, with $\tilde{D} : = S D$ , $- \tilde{C} {\tilde{D}}^{4} = \tilde{E} {\tilde{D}}^{2} \tilde{A} {\tilde{D}}^{2} - \tilde{F} {\tilde{D}}^{4} .$ (2.11)

We can now define the desired tensors E and F through $E : = S^{T} \tilde{E} S$ and $F : M ∋ b \mapsto S^{T} \cdot \tilde{F} (S b S^{T}) \cdot S \in M .$ (2.12) We note $\tilde{F}$ and F are related by the same formula as $\tilde{C}$ and C, cf. (2.9) and its inverse transformation. The calculation of (2.10) yields $\tilde{F} {\tilde{D}}^{4} = {F D}^{4}$ .

Regarding the operator $\tilde{E} {\tilde{D}}^{2}$ we calculate $\begin{array}{rcl} \tilde{E} {\tilde{D}}^{2} & = & \sum {(S E S^{T})}_{i j} {\tilde{D}}_{i} {\tilde{D}}_{j} = \sum S_{i α} E_{α β} S_{j β} S_{i k} D_{k} S_{j l} D_{l} \\ = & \sum E_{k l} D_{k} D_{l} = {E D}^{2} . \end{array}$ The calculation can also be applied to $A = S^{T} \tilde{A} S$ and provides $\tilde{A} {\tilde{D}}^{2} = {A D}^{2}$ . We conclude that relation (2.11) coincides with $- {C D}^{4} = {E D}^{2} {A D}^{2} - {F D}^{4} .$ This is the decomposition result for the tensors C and A. We remark that, since $S \in SO (n)$ is an orthogonal matrix, the symmetry and positive semi-definiteness of $\tilde{E}$ and $\tilde{F}$ carry over to E and F. This concludes the proof of the general decomposition lemma. □

3. An ε-independent effective equation

In the preceding sections we have seen that, to approximate solutions at the desired order, the heterogeneous wave equation can be replaced by a fourth order equations with constant coefficients. The result is unusual in the following sense: We are concerned with the limit $ε \to 0$ and propose a limit equation that still contains the small parameter $ε > 0$ .

Our aim in this section is a further analysis of the weakly dispersive effective equation. We will introduce moving frame coordinates which will allow us to follow the main pulse of $w^{ε}$ along rays through the origin. Performing the limit $ε \to 0$ in a distributional sense will provide an ε-independent linear third order equation (a linearized KdV-equation) that describes the effective shape of the pulse in dependence on the ray direction.

We emphasize that, in several respects, the weakly dispersive fourth order approximation is superior to the linearized KdV-approximation: The fourth order approximation is a well-posed system with solution $w^{ε}$ , and we have error estimates for $u^{ε} - w^{ε}$ . In contrast, the linearized KdV-approximation comes without an initial condition and without error terms. Indeed, even in the one-dimensional case, the initial condition for the limiting profile W can only be determined with a decomposition in a right-going and a left-going wave. In the higher-dimensional case, not even the distributional convergence as demanded in Proposition 3.1 is clear. Due to these restrictions, the linearized KdV-approximation can only be seen as a model for the further dispersion of travelling pulses after a long time.

In the following we will restrict ourselves to the analysis of the two-dimensional case with even symmetry, $a_{Y} (y_{1}, y_{2}) = a_{Y} (- y_{1}, y_{2}) = a_{Y} (y_{1}, - y_{2}) for all y \in R^{2} .$ (3.1) The above symmetry assumption, which is in particular satisfied in the case of a laminated structure, guarantees that the effective coefficients A and C have the form (cf. proof of Lemma 2.6 in [12]) $\begin{array}{rclr} A = diag (a_{1}, a_{2}), & (3.2) \\ C_{i i i i} = : α_{i}, C_{i j i j} = C_{i j j i} = C_{i i j j} = : β for i, j \in {1, 2} with i \neq j . & (3.3) \end{array}$ All other entries of C, that are not mentioned above, vanish.

We expect that the main pulse of $w^{ε}$ , solution to the weakly dispersive equation, propagates with a direction-dependent speed according to the anisotropic matrix $A = diag (a_{1}, a_{2})$ . We introduce appropriate elliptic coordinates $(r, φ)$ through $(x_{1}, x_{2}) = (r \sqrt{a_{1}} cos φ, r \sqrt{a_{2}} sin φ) .$ (3.4) The above coordinate transform is chosen in such a way that the main pulse of $w^{ε}$ at time t is located along the ellipse that is given as the level set ${x \in R^{2} ∣ r = t}$ . In order to perform a fine analysis of the dynamics of the pulse, we rewrite $w^{ε}$ as a function of $(r, φ, t)$ and use a moving frame in the radial variable r. More precisely, given the solution $w^{ε}$ , we define $W^{ε}$ through $W^{ε} (r, φ, t) : = \{\begin{matrix} w^{ε} (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) & for r > - t ε^{- 2}, \\ 0 & for r ⩽ - t ε^{- 2} . \end{matrix}$ (3.5) The time scaling $t / ε^{2}$ accounts for the fact that the dispersive effects of $w^{ε}$ are weak, i.e. slow in time. The main result of this section is the following. Provided that $W^{ε}$ has a distributional limit W, then W is characterized by an ε-independent linearized KdV-equation.

Proposition 3.1. (Effective behavior in the moving frame)

Let the medium $a_{Y} : R^{2} \to R^{2 \times 2}$ be evenly symmetric in the sense of (3.1). Let $w^{ε} (r, φ, t)$ be the solution to the weakly dispersive wave Eq. (2.3), expressed in elliptic coordinates. Let $W^{ε} (r, φ, t)$ be defined by (3.5). Assume that $W^{ε}$ has a limit in the sense of distributions, $W^{ε} ⇀ W in D^{'} (R \times R \times (0, T)) .$ Then the distribution $U : = \partial_{r} W$ satisfies the following linearized cylindrical Korteweg–de-Vries-equation(linearized KdV-equation) in distributional sense $\partial_{t} U + \frac{1}{2 t} U - \frac{1}{2} κ (φ) \partial_{r}^{3} U = 0 .$ (3.6) Here, the dispersion coefficient κ is given by $κ (φ) : = 6 β \frac{{cos}^{2} (φ)}{a_{1}} \frac{{sin}^{2} (φ)}{a_{2}} + α_{1} \frac{{cos}^{4} (φ)}{a_{1}^{2}} + α_{2} \frac{{sin}^{4} (φ)}{a_{2}^{2}} .$ (3.7)

The angle-dependent coefficient $κ (φ)$ can also be expressed with C, $κ (φ) = \sum_{i j k l} C_{i j k l} ξ_{i} ξ_{j} ξ_{k} ξ_{l} for ξ = (\frac{1}{\sqrt{a_{1}}} cos (φ), \frac{1}{\sqrt{a_{2}}} sin (φ)) .$ (3.8) It can be understood as a measure for the amount of dispersion along the ray $R_{φ} : = {x = (r \sqrt{a_{1}} cos φ, r \sqrt{a_{2}} sin φ) ∣ r \in (0, \infty)}$ . Note that, due to the negative semi-definiteness of C, the dispersion coefficient is always nonpositive: $κ (φ) ⩽ 0$ for every angle φ.

In the proof of the above proposition we will need some elementary formulas, collected in the following remark.

Remark 3.2. (Derivatives in elliptic coordinates)

Let $(r, φ)$ denote the elliptic coordinates defined in (3.4). Then the following relationship between derivatives in Cartesian and elliptic coordinates holds: $\begin{array}{rclr} \partial_{x_{1}} = \frac{1}{\sqrt{a_{1}}} (cos φ \partial_{r} - \frac{1}{r} sin φ \partial_{φ}), & (3.9) \\ \partial_{x_{2}} = \frac{1}{\sqrt{a_{2}}} (sin φ \partial_{r} + \frac{1}{r} cos φ \partial_{φ}), & (3.10) \\ {A D}^{2} = a_{1} \partial_{x_{1}}^{2} + a_{2} \partial_{x_{2}}^{2} = \partial_{r}^{2} + \frac{1}{r} \partial_{r} + \frac{1}{r^{2}} \partial_{φ}^{2} . & (3.11) \end{array}$ For general derivatives of order $γ = γ_{1} + γ_{2}$ with $γ_{1}, γ_{2} \in N_{0}$ one has $\begin{array}{rclr} \partial_{x_{1}}^{γ_{1}} \partial_{x_{2}}^{γ_{2}} & = & \frac{1}{{(\sqrt{a_{1}})}^{γ_{1}} {(\sqrt{a_{2}})}^{γ_{2}}} P^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}, \frac{1}{r}) \\ = & \frac{1}{{(\sqrt{a_{1}})}^{γ_{1}} {(\sqrt{a_{2}})}^{γ_{2}}} ({cos}^{γ_{1}} (φ) {sin}^{γ_{2}} (φ) \partial_{r}^{γ} + {\tilde{P}}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}, \frac{1}{r})), & (3.12) \end{array}$ where $P^{γ_{1}, γ_{2}}$ and ${\tilde{P}}^{γ_{1}, γ_{2}}$ are polynomials and ${\tilde{P}}^{γ_{1}, γ_{2}}$ can be written as ${\tilde{P}}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}, \frac{1}{r}) = \sum_{k = 1}^{γ} Q_{k}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}) \frac{1}{r^{k}},$ where $Q_{k}^{γ_{1}, γ_{2}}$ is a polynomial in $cos φ$ , $sin φ$ , $\partial_{r}$ , $\partial_{φ}$ of degree at most $2 γ$ .

Equations (3.9)–(3.11) are elementary identities for elliptic coordinates; (3.12) follows easily by an induction argument. In the decomposition of $P^{γ_{1}, γ_{2}}$ , we wrote the term of highest order in $\partial_{r}$ explicitly, with the result that all other terms (collected in ${\tilde{P}}^{γ_{1}, γ_{2}}$ ) contain non-vanishing powers of $\frac{1}{r}$ .

A consequence of the formulas (3.9)–(3.12) is the following lemma.

Lemma 3.3.
Let $w^{ε}$ and $W^{ε}$ be related by moving frame coordinates as in (3.5) and let W be a distributional limit of $W^{ε}$ as in Proposition 3.1. Let $γ_{1}, γ_{2} \in N_{0}$ and $γ = γ_{1} + γ_{2}$ . Then the following distributional convergence holds $(\partial_{x_{1}}^{γ_{1}} \partial_{x_{2}}^{γ_{2}} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) ⇀ \frac{{cos}^{γ_{1}} (φ) {sin}^{γ_{2}} (φ)}{{(\sqrt{a_{1}})}^{γ_{1}} {(\sqrt{a_{2}})}^{γ_{2}}} \partial_{r}^{γ} W$ for $ε \to 0$ in $D^{'} (R \times R \times (0, T))$ .
Proof.
We recall the definition $W^{ε} (r, φ, t) = w^{ε} (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})$ . The formula for $\partial_{x_{1}}^{γ_{1}} \partial_{x_{2}}^{γ_{2}}$ of Remark 3.2 provides $\begin{array}{rclr} (\partial_{x_{1}}^{γ_{1}} \partial_{x_{2}}^{γ_{2}} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) \\ = \frac{1}{{(\sqrt{a_{1}})}^{γ_{1}} {(\sqrt{a_{2}})}^{γ_{2}}} [{cos}^{γ_{1}} (φ) {sin}^{γ_{2}} (φ) (\partial_{r}^{γ} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) + R^{ε} (r, φ, t)], & (3.13) \end{array}$ where $R^{ε} (r, φ, t)$ is of the form $\begin{array}{rcl} R^{ε} (r, φ, t) & : = & \sum_{k = 1}^{γ} (Q_{k}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}) w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) \frac{1}{{(r + \frac{t}{ε^{2}})}^{k}} \\ = & \sum_{k = 1}^{γ} Q_{k}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}) W^{ε} (r, φ, t) {(\frac{ε^{2}}{ε^{2} r + t})}^{k} . \end{array}$ The last equality holds, since the derivatives $\partial_{r}$ , $\partial_{φ}$ act on $w^{ε}$ as they act on $W^{ε}$ .

Since $W^{ε}$ converges to W in the sense of distributions, also all derivatives converge in the distributional sense. We conclude $Q_{k}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}) W^{ε} (r, φ, t) ⇀ Q_{k}^{γ_{1}, γ_{2}} (cos φ, sin φ, \partial_{r}, \partial_{φ}) W$ in the distributional sense on $R \times R \times (0, T)$ for every k. In particular, exploiting that $\frac{ε^{2}}{ε^{2} r + t}$ is of order $ε^{2}$ for every $t > 0$ , we obtain $R^{ε} ⇀ 0$ in $D^{'} (R \times R \times (0, T))$ .

The same argument implies for the term containing only r-derivatives $(\partial_{r}^{γ} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) = (\partial_{r}^{γ} W^{ε}) (r, φ, t) ⇀ \partial_{r}^{γ} W in D^{'} (R \times R \times (0, T)) .$ This allows to pass to the distributional limit in (3.13), which concludes the proof of the lemma. □
Proof of Proposition 3.1.
We start from the weakly dispersive equation $\partial_{t}^{2} w^{ε} - {A D}^{2} w^{ε} - ε^{2} \partial_{t}^{2} {E D}^{2} w^{ε} + ε^{2} {F D}^{4} w^{ε} = 0 .$ The general idea of the proof is to transform the above equation into the elliptic coordinates of (3.4), to rewrite it in terms of $W^{ε} (r, φ, t) = w^{ε} (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})$ , and to pass to the distributional limit.

Step 1: The term $(\partial_{t}^{2} w^{ε} - {A D}^{2} w^{ε})$ . We start with the evaluation of time-derivatives of $w^{ε}$ , using the chain rule on $w^{ε} (r, φ, t) = W^{ε} (r - t, φ, t ε^{2})$ , $(\partial_{t}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) = \partial_{r}^{2} W^{ε} (r, φ, t) - 2 ε^{2} \partial_{t} \partial_{r} W^{ε} (r, φ, t) + ε^{4} \partial_{t}^{2} W^{ε} (r, φ, t) .$ (3.14) Combining this result with (3.11), we find $\begin{array}{rcl} (\partial_{t}^{2} w^{ε} - {A D}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) \\ = \partial_{r}^{2} W^{ε} (r, φ, t) - 2 ε^{2} \partial_{t} \partial_{r} W^{ε} (r, φ, t) + ε^{4} \partial_{t}^{2} W^{ε} (r, φ, t) \\ - (\partial_{r}^{2} W^{ε} (r, φ, t) + \frac{1}{r + \frac{t}{ε^{2}}} \partial_{r} W^{ε} (r, φ, t) + \frac{1}{{(r + \frac{t}{ε^{2}})}^{2}} \partial_{φ}^{2} W^{ε} (r, φ, t)) \\ = - ε^{2} (2 \partial_{t} \partial_{r} W^{ε} (r, φ, t) + \frac{1}{ε^{2} r + t} \partial_{r} W^{ε} (r, φ, t)) \\ + ε^{4} (\partial_{t}^{2} W^{ε} (r, φ, t) - \frac{1}{{(ε^{2} r + t)}^{2}} \partial_{φ}^{2} W^{ε} (r, φ, t)) . \end{array}$ We divide by $ε^{2}$ and take the distributional limit, exploiting that, by assumption, $W^{ε} ⇀ W$ in $D^{'} (R \times R \times (0, T))$ . We obtain $\frac{1}{ε^{2}} (\partial_{t}^{2} w^{ε} - {A D}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) ⇀ - 2 \partial_{t} \partial_{r} W - \frac{1}{t} \partial_{r} W$ (3.15) for $ε \to 0$ in $D^{'} (R \times R \times (0, T))$ .

Step 2: The term $- ε^{2} \partial_{t}^{2} {E D}^{2} w^{ε} + ε^{2} {F D}^{4} w^{ε}$ . In view of the scaling in (3.15), we have to analyze the distributional limit of $- {E D}^{2} \partial_{t}^{2} w^{ε} + {F D}^{4} w^{ε}$ .

We recall that, by our construction of E and F, the term can be written as ${C D}^{4} w^{ε} + R^{ε}$ with some remainder $R^{ε}$ that vanishes in the distributional limit as $ε \to 0$ . Indeed, exploiting the solution property of $w^{ε}$ and the decomposition (1.11) one obtains $- {E D}^{2} \partial_{t}^{2} w^{ε} + {F D}^{4} w^{ε} = {C D}^{4} w^{ε} + ε^{2} (- {E D}^{2} {E D}^{2} \partial_{t}^{2} w^{ε} + {E D}^{2} {F D}^{4} w^{ε}) .$ (3.16) In particular, inserting an argument in scaled variables, $\begin{array}{rclr} - ({E D}^{2} \partial_{t}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) + ({F D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) \\ = ({C D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) + R^{ε}, & (3.17) \end{array}$ where the remainder $R^{ε}$ is given by $\begin{array}{rcl} R^{ε} & = & ε^{2} [- ({E D}^{2} {E D}^{2} \partial_{t}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) + ({E D}^{2} {F D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})] \\ = & - ε^{2} (\partial_{r}^{2} - 2 ε^{2} \partial_{t} \partial_{r} + ε^{4} \partial_{t}^{2}) (({E D}^{2} {E D}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})) \\ + ε^{2} ({E D}^{2} {F D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) . \end{array}$ In the second equality we used relation (3.14). In this form we can apply Lemma 3.3 to $R^{ε}$ . We obtain that $({E D}^{2} {E D}^{2} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})$ and $({E D}^{2} {F D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})$ have distributional limits. In particular, taking into account the $ε^{2}$ -factor in the above formula, we conclude the convergence $R^{ε} ⇀ 0$ in $D^{'} (R \times R \times (0, T))$ as $ε \to 0$ .

Regarding the term $({C D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}})$ in (3.17) one can directly apply Lemma 3.3 with $γ = 4$ . Using the explicit form of C from (3.3) , we find $({C D}^{4} w^{ε}) (r + \frac{t}{ε^{2}}, φ, \frac{t}{ε^{2}}) ⇀ (6 β \frac{{cos}^{2} (φ)}{a_{1}} \frac{{sin}^{2} (φ)}{a_{2}} + α_{1} \frac{{cos}^{4} (φ)}{a_{1}^{2}} + α_{2} \frac{{sin}^{4} (φ)}{a_{2}^{2}}) \partial_{r}^{4} W$ in $D^{'} (R \times R \times (0, T))$ as $ε \to 0$ .

Step 3: Conclusion. Summing up the various terms, we find that the distribution W satisfies the equation $- 2 \partial_{t} \partial_{r} W - \frac{1}{t} \partial_{r} W + (6 β \frac{{cos}^{2} (φ)}{a_{1}} \frac{{sin}^{2} (φ)}{a_{2}} + α_{1} \frac{{cos}^{4} (φ)}{a_{1}^{2}} + α_{2} \frac{{sin}^{4} (φ)}{a_{2}^{2}}) \partial_{r}^{4} W = 0 .$ For $U : = \partial_{r} W$ we obtain the equation $\partial_{t} U + \frac{1}{2 t} U - \frac{1}{2} (6 β \frac{{cos}^{2} (φ)}{a_{1}} \frac{{sin}^{2} (φ)}{a_{2}} + α_{1} \frac{{cos}^{4} (φ)}{a_{1}^{2}} + α_{2} \frac{{sin}^{4} (φ)}{a_{2}^{2}}) \partial_{r}^{3} U = 0 .$ This concludes the proof of Proposition 3.1. □

4. Calculation of approximate solutions

In this section we discuss practical aspects of determining the coefficients in the effective dispersive Eq. (1.10). We present a numerical method and compare solutions $u^{ε}$ of the original wave Eq. (1.1) with solutions $w^{ε}$ of (1.10) .

4.1. Cell problems

The decomposition Lemmas 2.5 and 2.6 show that the coefficient tensors A, E and F of the weakly dispersive effective equation can be calculated from A and C, i.e. by the second and fourth derivatives of the Bloch eigenvalue $μ_{0} (k)$ at $k = 0$ . This fact makes the effective tensors computable in terms of cell-problems. In the following calculation we differentiate (1.5) with respect to k and integrate in y to obtain formulas for A and C. The result will be a practical algorithm to determine A and C. We will also transform Lemmas 2.5 and 2.6 into an algorithm that can be used to calculate E and F.

For a wave parameter $k \in Z = {(- 1 / 2, 1 / 2)}^{n}$ we consider the Bloch eigenvalue problem (1.5) for $m = 0$ . The smallest eigenvalue is $μ_{0} (k) \in R$ for every $k \in Z$ . A corresponding eigenfunction is $ψ_{0} = ψ_{0} (y, k)$ . Since a different normalization is more convenient in this section, we now use the letter Ψ; (1.5) then reads $- (\nabla_{y} + i k) \cdot (a_{Y} (y) (\nabla_{y} + i k) Ψ (y, k)) = μ_{0} (k) Ψ (y, k) .$ (4.1) While $ψ_{0}$ was normalized in the $L^{2} (Y)$ -norm, Ψ has a normalized average.

Remark 4.1. (Normalization)

For $k = 0$ , the eigenvalue is $μ_{0} (0) = 0$ and the eigenfunction $ψ_{0} (\cdot, 0)$ is a constant function. We normalize by setting $Ψ (\cdot, k) : = ψ_{0} (\cdot, k) {[\frac{1}{| Y |} \int_{Y} ψ_{0} (y, k) d y]}^{- 1}$ so that ${⟨ Ψ (\cdot, k) ⟩}_{Y} : = \frac{1}{| Y |} \int_{Y} Ψ (y, k) d y = 1 .$ (4.2) In particular, we obtain $Ψ (\cdot, 0) \equiv 1$ . The normalization is possible in a neighborhood of $k = 0$ , since averages of the first eigenfunction $ψ_{0} (\cdot, k)$ do not vanish for small $| k |$ . Due to the simplicity of $μ_{0}$ , the eigenvalue map $k \mapsto μ_{0} (k) \in R$ and the eigenfunction map $k \mapsto ψ_{0} (\cdot, k) \in L^{2} (Y)$ are analytic, see [11].

Due to Remark 4.1, it is legitimate to determine derivatives of $μ_{0}$ by differentiating the eigenvalue problem (4.1) in k. We use standard multi-index notation with $N_{0} = {0, 1, 2, \dots}$ : a multi-index $α = (α_{1}, \dots, α_{n}) \in N_{0}^{n}$ has length $| α | : = α_{1} + α_{2} + \dots + α_{n}$ and defines a differential operator $\partial^{α}$ of order $| α |$ (derivatives $\partial_{j} = \partial_{k_{j}}$ are with respect to $k \in R^{n}$ ), $\partial^{α} : = \partial_{1}^{α_{1}} \partial_{2}^{α_{2}} \dots \partial_{n}^{α_{n}}$ . We define, for $α \in N_{0}^{n}$ , $μ_{0}^{α} : = \partial^{α} μ_{0} |_{k = 0}, Ψ^{α} : = \partial^{α} Ψ |_{k = 0} .$ Differentiating the normalization (4.2), we obtain that the averages of the higher order derivatives vanish, ${⟨ Ψ^{α} ⟩}_{Y} = 0 for all α \in N_{0}^{n}, α \neq 0 .$ We additionally define the differential operators $\begin{array}{rcl} A (k) : = - (\nabla_{y} + i k) \cdot (a_{Y} (y) (\nabla_{y} + i k)), \\ A^{α} : = \partial^{α} A |_{k = 0} . \end{array}$

Lemma 4.2. (The operators

A^{α}

)

Let $e_{j}$ denote the jth Euclidean unit vector. The operators $A^{α}$ with $| α | ⩽ 2$ are $\begin{array}{rcl} A^{0} f = - \nabla_{y} \cdot (a_{Y} \nabla_{y} f), \\ A^{e_{j}} f = - i [(a_{Y} \nabla_{y} f) \cdot e_{j} + \nabla_{y} \cdot (a_{Y} e_{j} f)], \\ A^{e_{i} + e_{j}} f = 2 {(a_{Y})}_{i j} f \end{array}$ for $i, j = 1, \dots, n$ , where $f : R^{n} \to R$ is an arbitrary smooth function. All operators $A^{α}$ with $| α | ⩾ 3$ vanish identically.

Proof.
The formula for $A^{0}$ is obtained from the definition of $A (k)$ by setting $k = 0$ . Concerning the first order derivatives of $A (k)$ we calculate $\begin{array}{rcl} \partial_{j} A (k) f & = & - i e_{j} \cdot (a_{Y} (\nabla_{y} + i k) f) - i (\nabla_{y} + i k) \cdot (a_{Y} e_{j} f) \\ = & - i [e_{j} \cdot (a_{Y} (\nabla_{y} + i k) f) + (\nabla_{y} + i k) \cdot (a_{Y} e_{j} f)] . \end{array}$ Inserting $k = 0$ provides the claim about $A^{e_{j}}$ . For the second order derivatives of $A (k)$ we calculate $\partial_{i} \partial_{j} A (k) f = - i^{2} [e_{j} \cdot (a_{Y} e_{i} f) + e_{i} \cdot (a_{Y} e_{j} f)] = ({(a_{Y})}_{j i} + {(a_{Y})}_{i j}) f = 2 {(a_{Y})}_{i j} f,$ where the last equality holds due to the symmetry of $a_{Y}$ . □

We next want to obtain equations that characterize the functions $Ψ^{α}$ . In order to calculate derivatives of products, we will use the Leibniz formula in the following form: For every multi-index $α \in N_{0}^{n}$ and sufficiently smooth functions $f, g : R^{n} \to C$ there holds $\partial^{α} (f g) = \sum_{β \in N_{0}^{n}} (\begin{matrix} α \\ β \end{matrix}) \partial^{β} f \partial^{α - β} g .$ We use the binomial coefficient $(\binom{α}{β}) : = (\binom{α_{1}}{β_{1}}) (\binom{α_{2}}{β_{2}}) \dots (\binom{α_{n}}{β_{n}})$ , the (partial) ordering $β ⩽ α : \Leftrightarrow β_{i} ⩽ α_{i}$ for every $i \in {1, \dots, n}$ , and note that $(\binom{α}{β}) \neq 0$ is non-vanishing only for $β ⩽ α$ . In particular, in the following calculations, all sums over multi-indices are finite sums. Summation of multi-indices is performed in the standard way as $α \pm β : = (α_{1} \pm β_{1}, \dots, α_{n} \pm β_{n})$ .

With the operator $A (k)$ , we can write Eq. (4.1) as $A (k) Ψ (y, k) = μ_{0} (k) Ψ (y, k)$ . Taking partial derivatives with respect to k with the Leibniz formula, we find the following result.
Lemma 4.3. (Cell problems for $Ψ^{α}$ )

Let $α \in N_{0}^{n}$ be a multi-index. Then the function $Ψ^{α}$ satisfies the relation $A^{0} Ψ^{α} + \sum_{j = 1}^{n} α_{j} A^{e_{j}} Ψ^{α - e_{j}} + \sum_{i ⩽ j = 1}^{n} (\begin{matrix} α \\ e_{i} + e_{j} \end{matrix}) A^{e_{i} + e_{j}} Ψ^{α - e_{i} - e_{j}} = \sum_{β \in N_{0}^{n}} (\begin{matrix} α \\ β \end{matrix}) μ_{0}^{β} Ψ^{α - β} .$ Inserting the operators from Lemma 4.2 and using $μ_{0}^{0} = μ_{0} |_{k = 0} = 0$ , this equation reads $\begin{array}{rclr} - \nabla_{y} \cdot [a_{Y} \nabla_{y} Ψ^{α}] & = & i \sum_{j = 1}^{n} α_{j} [(a_{Y} \nabla_{y} Ψ^{α - e_{j}}) \cdot e_{j} + \nabla_{y} \cdot (Ψ^{α - e_{j}} a_{Y} e_{j})] \\ - 2 \sum_{i ⩽ j = 1}^{n} (\begin{matrix} α \\ e_{i} + e_{j} \end{matrix}) {(a_{Y})}_{i j} Ψ^{α - e_{i} - e_{j}} + \sum_{0 \neq β \in N_{0}^{n}} (\begin{matrix} α \\ β \end{matrix}) μ_{0}^{β} Ψ^{α - β} . & (4.3) \end{array}$

In our next step we obtain formulas for $μ_{0}^{α}$ in terms of $Ψ^{β}$ and $μ_{0}^{β}$ with $β < α$ .

Lemma 4.4. (Formulas for $μ_{0}^{α}$ )

Let $α \in N_{0}^{n}$ be a multi-index. In the case $α = 0$ we have $μ_{0}^{α} = μ_{0}^{0} = 0$ . For $| α | ⩾ 1$ there holds $\begin{array}{rclr} μ_{0}^{α} = 0 if | α | is odd, & (4.4) \\ μ_{0}^{α} = 2 \sum_{i ⩽ j = 1}^{n} (\begin{matrix} α \\ e_{i} + e_{j} \end{matrix}) {⟨ {(a_{Y})}_{i j} Ψ^{α - e_{i} - e_{j}} ⟩}_{Y} - i \sum_{j = 1}^{n} α_{j} {⟨ a_{Y} \nabla_{y} Ψ^{α - e_{j}} ⟩}_{Y} \cdot e_{j} if | α | is even . & (4.5) \end{array}$

Proof.
The statement for $μ_{0}^{0}$ has already been observed, cf. Remark 4.1. For odd $| α |$ , the symmetry $μ_{0} (- k) = μ (k)$ for all $k \in Z$ implies $μ_{0}^{α} = 0$ (compare e.g. Remark 2.7 in [12]). For even $| α |$ , the formula is a consequence of relation (4.3). Indeed, we can reorganize (4.3), writing the term with $β = α$ in the last sum explicitly, and find $\begin{array}{rcl} μ_{0}^{α} Ψ^{0} & = & - \nabla_{y} \cdot [a_{Y} \nabla_{y} Ψ^{α}] - i \sum_{j = 1}^{n} α_{j} [(a_{Y} \nabla_{y} Ψ^{α - e_{j}}) \cdot e_{j} + \nabla_{y} \cdot (Ψ^{α - e_{j}} a_{Y} e_{j})] \\ + 2 \sum_{i ⩽ j = 1}^{n} (\begin{matrix} α \\ e_{i} + e_{j} \end{matrix}) {(a_{Y})}_{i j} Ψ^{α - e_{i} - e_{j}} - \sum_{1 ⩽ | β | ⩽ | α | - 1} (\begin{matrix} α \\ β \end{matrix}) μ_{0}^{β} Ψ^{α - β} . \end{array}$

We integrate this relation over the periodicity cell Y, exploiting ${⟨ Ψ^{0} ⟩}_{Y} = 1$ and ${⟨ Ψ^{α} ⟩}_{Y} = 0$ for all $α \neq 0$ . Furthermore, we use that the integral over the two terms in divergence form vanishes by periodicity, and obtain (4.5). □

The above formulas allow to calculate all unknowns $Ψ^{α}$ and $μ_{0}^{α}$ in a recursive scheme. From Lemmas 4.3 and 4.4, we extract the following algorithm for the computation of the tensors A and C.
Algorithm 1. (Computation of A and C)

The tensors A and C can be computed in five steps.

For $j \in {1, \dots, n}$ , solve (4.3) for $α = e_{j}$ , i.e. $- \nabla_{y} \cdot [a_{Y} \nabla_{y} Ψ^{e_{j}}] = i \nabla_{y} \cdot (a_{Y} e_{j}) .$

For $i, j \in {1, \dots, n}$ , evaluate (4.5) for $α = e_{i} + e_{j}$ , i.e. $A_{i j} = \frac{1}{2} μ_{0}^{e_{i} + e_{j}} = {⟨ {(a_{Y})}_{i j} - \frac{i}{2} ((a_{Y} \nabla_{y} Ψ^{e_{j}}) \cdot e_{i} + (a_{Y} \nabla_{y} Ψ^{e_{i}}) \cdot e_{j}) ⟩}_{Y} .$

For $i, j \in {1, \dots, n}$ , solve (4.3) for $α = e_{i} + e_{j}$ , i.e. $\begin{array}{rcl} - \nabla_{y} \cdot [a_{Y} \nabla Ψ^{e_{i} + e_{j}}] & = & i [(a_{Y} \nabla_{y} Ψ^{e_{j}}) \cdot e_{i} + (a_{Y} \nabla_{y} Ψ^{e_{i}}) \cdot e_{j} \\ + \nabla_{y} \cdot (Ψ^{e_{j}} a_{Y} e_{i}) + \nabla_{y} \cdot (Ψ^{e_{i}} a_{Y} e_{j})] - 2 {(a_{Y})}_{i j} + μ_{0}^{e_{i} + e_{j}} . \end{array}$

For $i, j, k \in {1, \dots, n}$ , solve (4.3) for $α = e_{i} + e_{j} + e_{k}$ , $\begin{array}{rcl} - \nabla_{y} \cdot [a_{Y} \nabla_{y} Ψ^{e_{i} + e_{j} + e_{k}}] \\ = i [(a_{Y} \nabla_{y} Ψ^{e_{j} + e_{k}}) \cdot e_{i} + (a_{Y} \nabla_{y} Ψ^{e_{i} + e_{k}}) \cdot e_{j} + (a_{Y} \nabla_{y} Ψ^{e_{i} + e_{j}}) \cdot e_{k} \\ + \nabla_{y} \cdot (Ψ^{e_{j} + e_{k}} a_{Y} e_{i}) + \nabla_{y} \cdot (Ψ^{e_{i} + e_{k}} a_{Y} e_{j}) + \nabla_{y} \cdot (Ψ^{e_{i} + e_{j}} a_{Y} e_{k})] \\ - 2 ({(a_{Y})}_{i j} Ψ^{e_{k}} + {(a_{Y})}_{i k} Ψ^{e_{j}} + {(a_{Y})}_{j k} Ψ^{e_{i}}) \\ + μ_{0}^{e_{i} + e_{j}} Ψ^{e_{k}} + μ_{0}^{e_{i} + e_{k}} Ψ^{e_{j}} + μ_{0}^{e_{j} + e_{k}} Ψ^{e_{i}} . \end{array}$

For $i, j, k, l \in {1, \dots, n}$ , evaluate (4.5) for $α = e_{i} + e_{j} + e_{k} + e_{l}$ to find $\begin{array}{rcl} C_{i j k l} & = & \frac{1}{24} μ_{0}^{e_{i} + e_{j} + e_{k} + e_{l}} \\ = & \frac{1}{24} ⟨ 2 ({(a_{Y})}_{i j} Ψ^{e_{k} + e_{l}} + {(a_{Y})}_{i k} Ψ^{e_{j} + e_{l}} + {(a_{Y})}_{i l} Ψ^{e_{j} + e_{k}} \\ + {(a_{Y})}_{j k} Ψ^{e_{i} + e_{l}} + {(a_{Y})}_{j l} Ψ^{e_{i} + e_{k}} + {(a_{Y})}_{k l} Ψ^{e_{i} + e_{j}}) \\ - i ((a_{Y} \nabla_{y} Ψ^{e_{j} + e_{k} + e_{l}}) \cdot e_{i} + (a_{Y} \nabla_{y} Ψ^{e_{i} + e_{k} + e_{l}}) \cdot e_{j} \\ + (a_{Y} \nabla_{y} Ψ^{e_{i} + e_{j} + e_{l}}) \cdot e_{k} + (a_{Y} \nabla_{y} Ψ^{e_{i} + e_{j} + e_{k}}) \cdot e_{l}) ⟩_{Y} . \end{array}$

The elliptic problems in Steps 1,3 and 4 are posed on Y with periodic boundary conditions and with the condition of zero mean, i.e.

{⟨ Ψ^{e_{j}} ⟩}_{Y} = {⟨ Ψ^{e_{i} + e_{j}} ⟩}_{Y} = {⟨ Ψ^{e_{i} + e_{j} + e_{k}} ⟩}_{Y} = 0 .

In the above algorithm we have used the fact that $μ_{0}^{α} = 0$ for all α with $| α | = 1$ and $| α | = 3$ . The above cell-problems are complex valued, but the resulting tensors A and C are real. The fact that values of $Ψ^{α}$ and $μ_{0}^{α}$ do not change upon permutations of the entries $α_{1}, α_{2}, \dots, α_{n}$ , allows to reduce the number of problems to be solved in Steps 3 and 4: The relevant number is $(\binom{n + 1}{2})$ and $(\binom{n + 2}{3})$ rather than $n^{2}$ and $n^{3}$ . Moreover, spatial symmetries of $a_{Y} (\cdot)$ can reduce the number of problems further: If $a_{Y} (\cdot)$ is even in both variables, $a_{Y} (y_{1}, y_{2}) = a_{Y} (- y_{1}, y_{2}) = a_{Y} (y_{1}, - y_{2})$ for all $y \in Y$ , then all derivatives of $μ_{0}$ at $k = 0$ involving an odd number of derivatives in one direction vanish, and only $A_{i i}$ , $C_{i i i i}$ , and $C_{i i j j} = C_{i j i j} = C_{j i i j}$ for $i, j \in {1, \dots, n}$ are potentially nonzero.

At this point, we have presented a scheme to compute the effective tensors A and C with the help of a sequence of cell-problems. We conclude this section with the outline of an algorithm (based on the proofs of Lemmas 2.5 and 2.6) that provides formulas for the effective coefficient tensors E and F.

Algorithm 2. (Computation of E and F)

The tensors E and F can be computed in four steps.

Determine $S \in SO (n)$ such that $A = S^{T} \tilde{A} S$ with $\tilde{A} = diag (a_{1}, a_{2}, \dots, a_{n})$ . Define $\tilde{C}$ by (2.9) .

Loop over all indices $1 ⩽ α, β, γ, δ ⩽ n$ such that no two indices coincide (an empty set in dimension $n ⩽ 3$ ), use (2.7) to compile $\hat{F}$ . Define a tensor $\bar{C}$ such that $\bar{C} {\tilde{D}}^{4} = (C - \hat{F}) {\tilde{D}}^{4}$ . The tensor $\bar{C}$ can be chosen such that entries with four different indices vanish, see (2.8) .

Loop over all remaining indices $1 ⩽ α, β, γ, δ ⩽ n$ . Corresponding to $\tilde{A}$ and $\bar{C}$ , compile $\bar{E}$ and $\bar{F}$ from the explicit formulas of Cases 1–3 in the proof of Lemma 2.6 .

Set $\tilde{E} : = \bar{E}$ and $\tilde{F} = \bar{F} + \hat{F}$ . Obtain $E : = S^{T} \tilde{E} S$ and F from (2.12) .

4.2. Numerical results in 2D

In the following, we present numerical results for all the three parts of the homogenization problem: The computation of the original ε-problem, the computation of effective coefficients with the help of cell-problems, and the computation of the weakly dispersive effective problem. The methods vary, we use finite element schemes and finite difference schemes, see e.g. [22]. We refer also to [1] for a recent analysis of numerical methods and further references. Our results show an excellent agreement between solutions $u^{ε}$ of (1.1) and solutions $w^{ε}$ of (1.10).

The initial conditions for all the tests below are $u^{ε} (x, 0) = w^{ε} (x, 0) = e^{- 4 (x_{1}^{2} + x_{2}^{2})}, \partial_{t} u^{ε} (x, 0) = \partial_{t} w^{ε} (x, 0) = 0 .$ (4.6) The periodicity cell is $Y = {[- π, π]}^{2}$ in the y-variables and $ε Y$ in the x-variables. All the tests are carried out for the case of even symmetry in $a_{Y}$ , i.e. $a_{Y} (- y_{1}, y_{2}) = a_{Y} (y_{1}, - y_{2}) = a_{Y} (y_{1}, y_{2})$ for all $y \in Y$ , so that for the even initial data in (4.6) problems (1.1) and (1.10) can be reduced to one quadrant with homogeneous Neumann boundary conditions along the coordinate axes $x_{1} = 0$ , $x_{2} = 0$ . The computational domain Ω is rectangular with two sides coinciding with the negative $x_{1}$ and $x_{2}$ axes. At the remaining two sides of the rectangle we use homogeneous Dirichlet boundary conditions. The size of Ω is chosen so that the solution remains localized in Ω until the final computation time. The function $a_{Y}$ is chosen piecewise constant in all tests; the observed error convergence agrees with (2.1) although Theorem 2.2 treats only differentiable fields $a_{Y}$ .

The even symmetry of $a_{Y}$ in both $y_{1}$ and $y_{2}$ causes that A is diagonal and C has only eight nonzero entries, see (3.2) and (3.3): $\begin{array}{rcl} A_{i i} = : a_{i}, C_{i i i i} = : α_{i}, i \in {1, 2}, \\ C_{1122} = C_{2211} = C_{1212} = C_{2121} = C_{1221} = C_{2112} = : β, \end{array}$ with all other entries of A and C zero. A choice of E and F according to Algorithm 2 is $\begin{matrix} E_{11} = \frac{{- α_{1}}_{+}}{a_{1}} + 3 \frac{{- β}_{+}}{a_{2}}, E_{22} = \frac{{- α_{2}}_{+}}{a_{2}} + 3 \frac{{- β}_{+}}{a_{1}}, E_{12} = E_{21} = 0, \\ F_{1111} = {α_{1}}_{+} + 3 \frac{a_{1}}{a_{2}} {- β}_{+}, F_{2222} = {α_{2}}_{+} + 3 \frac{a_{2}}{a_{1}} {- β}_{+}, \\ F_{2121} = \frac{a_{1}}{a_{2}} {- α_{2}}_{+} + 3 {β}_{+}, F_{1212} = \frac{a_{2}}{a_{1}} {- α_{1}}_{+} + 3 {β}_{+} \end{matrix}$ (4.7) with all other entries of F being zero. We use these tensors in the numerical tests further. Note that $α_{i} ⩽ 0$ holds for $i \in {1, 2}$ , hence ${- α_{i}}_{+} = - α_{i}$ and ${α_{i}}_{+} = 0$ .

4.2.1. Numerical method

The values $a_{j}$ , $α_{j}$ and β for $j = 1, 2$ are computed via Algorithm 1, where the elliptic equations in Steps 1, 3 and 4 are discretized by linear finite elements using the PDE-Toolbox of Matlab. The periodic boundary conditions are implemented by modifying the stiffness matrix and the load vector corresponding to homogeneous Neumann boundary conditions. In all tests a uniform discretization conforming to the material geometry is generated with the Matlab function poimesh, i.e. the elements are all right angled, have equal size and the discontinuity lines of $a_{Y}$ intersect no elements. The element size is given by specifying the spacings $h_{1}$ and $h_{2}$ , i.e. the lengths of the triangle legs in poimesh.

The original wave Eq. (1.1) is discretized in space also via linear finite elements using the PDE-Toolbox of Matlab with uniform elements (generated by poimesh) conforming to the geometry. The values of $h_{1}$ and $h_{2}$ are given for each example further.

For the weakly dispersive problem (1.10), which has constant coefficients, we use the fourth order centered finite difference discretization in space for both the second order and fourth order derivatives. In all tests we use the spacing $d x_{1} = d x_{2} = 0.2$ for (1.10) .

The time discretization of both (1.1) and (1.10) is done via the second order Leap–Frog method in two step formulation, e.g. for (1.1) the semi-discrete problem is thus $u^{(n + 1)} = 2 u^{(n)} - u^{(n - 1)} + {(d t)}^{2} \nabla \cdot (a^{ε} \nabla u^{(n)}),$ where $u^{(n)} \approx u^{ε} (t = n \cdot d t)$ . For the initialization of the scheme we use the second order Taylor expansion $u^{(1)} = u^{(0)} + {(d t)}^{2} \nabla \cdot (a^{ε} \nabla u^{(0)}),$ according to the initial condition $\partial_{t} u^{ε} (x, 0) = 0$ . For (1.1) we use the time step $d t = min {0.01, h_{1} / 4, h_{2} / 4}$ , for (1.10) we use $d t = 0.02$ .

4.2.2. ε-Convergence of the error

Here, our aim is to determine experimentally the convergence rate of the approximation error ${∥ u^{ε} (t) - w^{ε} (t) ∥}_{L^{2} (Ω)}$ . To this end we fix a periodic coefficient matrix field $a_{Y}$ . With the identity matrix $I : = {id}_{R^{2}}$ we set $a_{Y} (y) = {\tilde{a}}_{Y} (y) I$ , where ${\tilde{a}}_{Y}$ is defined through a rectangular geometry, $\begin{matrix} {\tilde{a}}_{Y} (y) = \frac{1}{2} + b (y) - \frac{1}{| Y |} \int_{Y} b (y) d y, \\ b (y) = \{\begin{matrix} 1.6 & for y \in [- \frac{11 π}{13}, \frac{11 π}{13}] \times [- \frac{π}{3}, \frac{π}{3}], \\ 0.2 & otherwise, \end{matrix} \end{matrix}$ (4.8) which is illustrated in Fig. 3(a).

Fig. 3.

(a) Illustration of the geometry and of the function $a_{Y} (\cdot)$ from (4.8) . (b) The solution $u^{ε}$ of the wave equation in a highly oscillatory medium. (c) The solution $w^{ε}$ of the weakly dispersive wave equation with constant coefficients. The geometry is as in (4.8), the plots are for $ε = 0.07$ and $t = 100$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

Our best numerical approximation of the effective coefficients is $a_{1} \approx 0.281, a_{2} \approx 0.179, α_{1} \approx - 0.273, α_{2} \approx - 0.044, β \approx 0.024 .$ These values have been computed with $h_{1} = 2 π / 208$ and $h_{2} = 2 π / 192$ . Within rounding to three decimal places the values do not change with a further mesh refinement.

We solve (1.1) and (1.10) for the values $ε = 0.2, 0.17, 0.15, 0.12, 0.1, 0.07$ up to times $t = 12.5, 17.3, 22.5, 35, 50, 100 \approx ε^{- 2} / 2$ , respectively. To keep the computational expense within limits, we use the discretization $h_{1} = 2 π ε / 13$ , $h_{2} = 2 π ε / 12$ in the simulations of (1.1). Using the same number of uniform elements in the cell problems as in each periodic cell in (1.1), i.e. discretizing Y by $2 \times 13 \times 12 = 312$ uniform elements, we obtain the values $a_{1} \approx 0.2784, a_{2} \approx 0.1506, α_{1} \approx - 0.369, α_{2} \approx - 0.034, β \approx 0.032 .$ (4.9) We calculate solutions to the weakly dispersive Eq. (1.10) with the values (4.9) rather than with the converged coefficient values; this is justified by the fact that (4.9) are the effective coefficients for the particular discretization. Formulas (4.7) and the values in (4.9) produce the effective coefficients $\begin{array}{rcl} E_{11} = 1.3256, E_{22} = 0.2257, \\ F_{1111} = F_{2222} = 0, F_{2121} = 0.1588, F_{1212} = 0.2957 . \end{array}$ We can now investigate the convergence of the error in ε: The $L^{2} (Ω)$ -error ${∥ u^{ε} (t) - w^{ε} (t) ∥}_{L^{2} (Ω)}$ at the final time $t = ε^{- 2} / 2$ is plotted in Fig. 2. As Fig. 2(b) shows, the convergence is close to linear in accordance with estimate (2.1).

In Fig. 3 the solutions $u^{ε}$ and $w^{ε}$ are plotted for $ε = 0.07$ at $t = 100$ . In both plots we see clearly the main pulse located along an ellipse, we see the dispersive oscillations behind the main pulse, and we see that the dispersion is weakest along a ray that has an approximate angle $π / 4$ . Due to the weak dispersion along this ray, the main pulse has its maximal amplitude in this direction (compare also Fig. 6). An excellent agreement of the two calculations is observed.

Fig. 4.

(a) Illustration of the geometry that defines the function $a_{Y} (\cdot)$ in (4.10) . (b) The solution $u^{ε} (x, t = 40)$ for $ε = 0.07$ . (c) The solution $w^{ε} (x, t = 40)$ for $ε = 0.07$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

Fig. 5.

(a) Illustration of the function $a_{Y} (\cdot)$ of (4.12). (b) The solution $u^{ε} (x, t = 40)$ for $ε = 0.07$ . (c) The solution $w^{ε} (x, t = 40)$ for $ε = 0.07$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

4.2.3. Further numerical examples

We consider two more geometries. For a cross-shaped geometry as illustrated in Fig. 4(a), we set $a_{Y} (y) : = {\tilde{a}}_{Y} (y) I$ with ${\tilde{a}}_{Y} (y) = \{\begin{matrix} 2 & for y \in [- \frac{7 π}{9}, \frac{7 π}{9}] \times [- \frac{2 π}{9}, \frac{2 π}{9}] \cup [- \frac{2 π}{9}, \frac{2 π}{9}] \times [- \frac{7 π}{9}, \frac{7 π}{9}], \\ 0.2 & otherwise . \end{matrix}$ (4.10)

The additional symmetry $a_{Y} (y_{1}, y_{2}) = a_{Y} (y_{2}, y_{1})$ for all $y \in Y$ implies the relations $a_{1} = a_{2}$ and $α_{1} = α_{2}$ , see Lemma 2.6 in [12]. Algorithm 1 provides, with Y discretized by $2 \times 18^{2} = 648$ uniform elements of size $h_{1} = h_{2} = 2 π / 18$ , the values $a_{1} = a_{2} \approx 0.3816, α_{1} = α_{2} \approx - 0.1970, β \approx 0.0394 .$ (4.11) To check the accuracy, we calculated the values also using a fine resolution with $2 \times 360 \times 360$ uniform elements. The fine resolution provides $a_{1} = a_{2} \approx 0.406$ , $α_{1} = α_{2} \approx - 0.235$ , $β \approx 0.044$ , and the first three decimal places do not change upon further refinement. Similarly to the example in Section 4.2.2, the discretization error in (4.11) is quite large. Nevertheless, we use these coefficients for the calculation of the solution $w^{ε}$ of (1.10). The solution is compared to $u^{ε}$ , which is computed with the corresponding spatial discretization $h_{1} = h_{2} = 2 π ε / 18$ . The results at time $t = 40$ for $ε = 0.07$ are plotted in Fig. 4(b) and (c).

We finally consider a laminated structure $a_{Y} (y) = {\tilde{a}}_{Y} (y) I$ with ${\tilde{a}}_{Y} (y) = \{\begin{matrix} 2 & for y \in [- π, π] \times [- \frac{2 π}{5}, \frac{2 π}{5}], \\ 0.2 & otherwise, \end{matrix}$ (4.12) compare Fig. 5(a). Effective coefficients are obtained by Algorithm 1, Y is discretized with $2 \times 12 \times 16 = 384$ uniform elements ( $h_{1} = 2 π / 12$ , $h_{2} = 2 π / 16$ ): $\begin{matrix} a_{1} \approx 0.8750, a_{2} \approx 0.3019, \\ α_{1} \approx - 1.9185, α_{2} \approx - 0.0933, β \approx 0.1448 . \end{matrix}$ (4.13) The converged values (with four reliable digits, computed with $2 \times 140 \times 180$ uniform elements) are: $a_{1} \approx 0.9200$ , $a_{2} \approx 0.3125$ , $α_{1} \approx - 1.9645$ , $α_{2} \approx - 0.1170$ , $β \approx 0.1599$ . The effective coefficients determined using (4.7) and (4.13) are $\begin{array}{rcl} E_{11} = 2.1925, E_{22} = 0.3091, \\ F_{1111} = F_{2222} = 0, F_{2121} = 0.7050, F_{1212} = 1.0964 . \end{array}$ Equation (1.1) was discretized in space using $h_{1} = 2 π ε / 12$ and $h_{2} = 2 π ε / 16$ , the solutions $u^{ε}$ and $w^{ε}$ are computed using the above coefficients; they are plotted for $t = 40$ and $ε = 0.07$ in Fig. 5(b) and (c).

Fig. 6.

Solutions $u^{ε}$ and $w^{ε}$ for $ε = 0.07$ at $t = 100$ for the rectangle geometry (4.8) along two rays: (a) $ϕ = 0$ . (b) $ϕ = ϕ_{m} \approx π / 4 + 0.002$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

4.2.4. Dependence of the dispersion on the propagation angle

The rate of dispersion depends on the angle of propagation. Here, we have to distinguish the angle φ of the elliptic coordinates and the corresponding angle ϕ in polar coordinates (describing the observable angle of the ray). We measure the angles such that $ϕ = 0$ corresponds to the negative $x_{1}$ -axis. The two angles are related by $tan (ϕ) = \sqrt{a_{2} / a_{1}} tan (φ)$ . The dependence of the dispersion on the angle can be observed in the numerical results: we see few oscillations in a propagation angle of approximately $ϕ = π / 4$ , more oscillations along rays that are aligned with the coordinate axes, i.e. at angles $ϕ = 0$ and $ϕ = π / 2$ . The angle of minimal dispersion can be obtained by minimizing $κ = κ (φ)$ of (3.7) .

To illustrate the angular dependence, we consider the rectangular geometry (4.8) and the two rays corresponding to $φ = 0$ and $φ = φ_{m}$ , the minimizer of $κ (\cdot)$ . We plot $u^{ε}$ and $w^{ε}$ for $ε = 0.07$ and $t = 100$ in Fig. 6. One can see clearly a much smaller dispersion at the angle $ϕ_{m} \approx π / 4 + 0.002$ corresponding to $φ_{m}$ . Additionally, we observe a much larger error in the aligned direction $φ = 0$ . The values of the dispersion coefficient are $κ (0) \approx - 4.762$ and $κ (φ_{m}) \approx - 0.175$ .

Fig. 7.

Solutions $u^{ε}$ and $w^{ε}$ for $ε = 0.07$ at $t = 40$ for the laminated geometry (4.12) along three rays: (a) $ϕ = 0$ . (b) $ϕ = π / 2$ . (c) $ϕ = ϕ_{m} \approx π / 4 - 0.153$ . (Colors are visible in the online version of the article; http://dx.doi.org/10.3233/ASY-141280.)

For the laminate structure (4.12) we compare the solutions along three directions in Fig. 7, namely the horizontal direction $ϕ = 0$ along which the structure is constant, the vertical direction $ϕ = π / 2$ , orthogonal to the laminates, and an intermediate direction $ϕ_{m}$ that corresponds to the minimizer $φ_{m} \approx π / 4 - 0.153$ of κ. The values of the dispersion coefficient are $κ (0) \approx - 2.506$ , $κ (π / 2) \approx - 1.024$ and $κ (φ_{m}) \approx 0.02$ . The analytical values of κ are always non-positive, hence the positive value $κ (φ_{m})$ is caused by discretization errors.

References

[1]

Abdulle,

M.J.

Grote and

Stohrer, FE heterogeneous multiscale method for long-time wave propagation, C. R. Math. Acad. Sci. Paris 351(11,12) (2013), 495–499.

[2]

Allaire, Dispersive limits in the homogenization of the wave equation, Ann. Fac. Sci. Toulouse Math. (6) 12(4) (2003), 415–431.

[3]

Allaire,

Conca and

Vanninathan, The Bloch transform and applications, in: Actes du 29ème Congrès D’Analyse Numérique: CANum’97, Larnas, 1997, ESAIM: Proc., Vol. 3, Soc. Math. Appl. Indust., Paris, 1998, pp. 65–84 (electronic).

[4]

Allaire,

Palombaro and

Rauch, Diffractive behavior of the wave equation in periodic media: Weak convergence analysis, Ann. Mat. Pura Appl. (4) 188(4) (2009), 561–589.

[5]

Allaire,

Palombaro and

Rauch, Diffractive geometric optics for Bloch wave packets, Arch. Ration. Mech. Anal. 202(2) (2011), 373–426.

[6]

Brahim-Otsmane,

G.A.

Francfort and

Murat, Correctors for the homogenization of the wave and heat equations, J. Math. Pures Appl. (9) 71(3) (1992), 197–231.

[7]

Castro and

Zuazua, Concentration and lack of observability of waves in highly heterogeneous media, Arch. Ration. Mech. Anal. 164(1) (2002), 39–72.

[8]

Conca,

Orive and

Vanninathan, Bloch approximation in homogenization and applications, SIAM J. Math. Anal. 33(5) (2002), 1166–1198 (electronic).

[9]

Conca,

Orive and

Vanninathan, Bloch approximation in homogenization on bounded domains, Asymptot. Anal. 41(1) (2005), 71–91.

10.

[10]

Conca,

Orive and

Vanninathan, On Burnett coefficients in periodic media, J. Math. Phys. 47(3) (2006), 032902.

11.

[11]

Conca and

Vanninathan, Homogenization of periodic structures via Bloch decomposition, SIAM J. Appl. Math. 57(6) (1997), 1639–1659.

12.

[12]

Dohnal,

Lamacz and

Schweizer, Bloch-wave homogenization on large time scales and dispersive effective wave equations, Multiscale Model. Simul. 12(2) (2014), 488–513.

13.

[13]

Fish and

Chen, Space–time multiscale model for wave propagation in heterogeneous media, Comput. Methods Appl. Mech. Engrg. 193(45–47) (2004), 4837–4856.

14.

[14]

Fish,

Chen and

Nagai, Uniformly valid multiple spatial–temporal scale modeling for wave propagation in heterogeneous media, Mechanics of Composite Materials and Structures 8 (2001), 81–99.

15.

[15]

Fish,

Chen and

Nagai, Non-local dispersive model for wave propagation in heterogeneous media: Multi-dimensional case, Internat. J. Numer. Methods Engrg. 54(3) (2002), 347–363.

16.

[16]

Fish,

Chen and

Nagai, Non-local dispersive model for wave propagation in heterogeneous media: One-dimensional case, Internat. J. Numer. Methods Engrg. 54(3) (2002), 331–346.

17.

[17]

G.A.

Francfort and

Murat, Oscillations and energy densities in the wave equation, Comm. Partial Differential Equations 17(11,12) (1992), 1785–1865.

18.

[18]

Lamacz, Dispersive effective models for waves in heterogeneous media, Math. Models Methods Appl. Sci. 21(9) (2011), 1871–1899.

19.

[19]

Lombardo and

Askes, Elastic wave dispersion in microstructured membranes, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 466(2118) (2010), 1789–1807.

20.

[20]

Orive,

Zuazua and

A.F.

Pazoto, Asymptotic expansion for damped wave equations with periodic coefficients, Math. Models Methods Appl. Sci. 11(7) (2001), 1285–1310.

21.

[21]

Pichugin,

Askes and

Tyas, Asymptotic equivalence of homogenisation procedures and fine-tuning of continuum theories, Journal of Sound and Vibration 313(35) (2008), 858–874.

22.

[22]

Quarteroni,

Sacco and

Saleri, Numerical Mathematics. Texts in Applied Mathematics, Springer, 2007.

23.

[23]

Santosa and

W.W.

Symes, A dispersive effective medium for wave propagation in periodic composites, SIAM J. Appl. Math. 51(4) (1991), 984–1005.

Dispersive homogenized models and coefficient formulas for waves in general periodic media

Abstract

Keywords

1. Introduction

Decomposition lemma and approximation result

An ε-independent third order dispersive equation

Algorithms to compute effective quantities

2. A weakly dispersive effective equation

2.1. Main approximation result

Theorem 2.3. (Corollary 2.5 of [12])

Lemma 2.6. (Decomposability for diagonal matrices A)

4.1. Cell problems

Lemma 4.4. (Formulas for μ 0 α )

Algorithm 2. (Computation of E and F)

4.2. Numerical results in 2D

4.2.2. ε-Convergence of the error

References

Lemma 4.4. (Formulas for $μ_{0}^{α}$ )