Numerical analysis of a wave equation for lossy media obeying a frequency power law

We study a wave equation with a nonlocal time fractional damping term that models the effects of acoustic attenuation characterized by a frequency dependence power law. First we prove existence of unique solutions to this equation with particular attention paid to the handling of the fractional derivative. Then we derive an explicit time stepping scheme based on the finite element method in space and a combination of convolution quadrature and second order central differences in time. We conduct a full error analysis of the mixed time discretization and in turn the fully space time discretized scheme. A number of numerical results are presented to support the error analysis for both smooth and nonsmooth solutions.

The interest in this problem stems from the modelling of acoustic attenuation that occurs as a wave propagates through lossy media [28,Chapter 4]. In applications, this is applied to modelling high intensity focused ultrasound therapy (HIFU) where the lossy media is biological tissue. HIFU is a noninvasive and nonionising medical treatment that focuses multiple high intensity acoustic pressure waves on a region, ablating it away ( [30]). It is known that acoustic attenuation obeys a frequency dependence characterized by the following power law, where S is the amplitude, ∆ − → x is the wave propagation distance, ω denotes frequency and α(ω) the attenuation coefficient which defined by α(ω) = α 0 |ω| γ for γ the frequency power exponent and α 0 a constant relating to the media ( [7]). Values of γ have been determined by many experiments and field measurements, their results conclude that for most media γ ∈ (0, 2) ( [10]). In this paper we restrict our attention to γ ∈ (−1, 1). For the range γ ∈ (1, 2) (and potentially also (0, 1)) we prefer changing the weak damping term ∂ γ+1 t with the strong damping of the form (−∆)∂ γ−1 t ( [6,12,13,32]). The techniques developed in this paper, can be used to analyze the strongly damped case for the full range γ ∈ (0, 2); see ( [1]). Another model that we can consider is the nonlinear fractional Westervelt equation ( [12]), this contains the strong damping term mentioned above and describes the nonlinear behaviour that the high intensity focusing causes during HIFU. Similar models appear also in fractional order viscoelasticity, similar to both the strong damping case [14,25] and the weak damping case investigated here [23]. [27] derived a wave equation with a convolution type operator to incorporate the effects of acoustic attenuation and showed that it adheres to the frequency dependence power law (3). [7] praise this model for its simplicity, due to it only containing two parameters, but criticise it for being difficult to implement initial conditions. To overcome this they refined the model to include the Caputo fractional derivative as the damping term. Furthermore within their paper they show that even with this modification the solutions still obey the power law.
Issues that arise within this model largely stem from the convolution based fractional derivative that requires us to store the full history when numerically solving with a time stepping scheme, making it expensive to compute, particularly in 3D. In the literature this is countered by replacing the fractional time derivative with a fractional Laplacian ( [8,31]). However, using efficient quadrature methods ( [2]) will also be sufficient in reducing memory requirements and allows us to remain using the simpler operator.
The paper is structured in the following way. Section 2 begins by outlining necessary definitions and lemmas required to prove existence and uniqueness of solutions to (1). In Section 3 we develop a numerical scheme to approximate solutions of (1) using finite element methods in space and a combination of second order central difference and convolution quadrature in time. In Section 4 we study the convergence analysis of the mixed time dicretization used which allows us to determine the convergence order of the full scheme. Lastly, in Section 5 we show results from the implementation of our scheme to support the theory from the previous section for both smooth and non smooth solutions of (1).

A damped wave equation model
In this section we describe the model with the fractional in time weak damping and prove existence and uniqueness of the solution. In order to do this we will require some properties of fractional calculus, that we list first. For more detail on fractional calculus see textbooks by [9] and [22].
-The space of H 1 (Ω) functions with a zero trace on ∂Ω is denoted by H 1 0 (Ω). We first give the definition of the Riemann Lioville fractional integral. Definition 2.1. For β > 0 the Riemann Liouville fractional integral of a function f is defined by where Γ denotes the Gamma function.
If f (t) = t µ for µ > −1 we have that The Caputo fractional derivative is obtained by applying the Riemann Liouville fractional integral to a classical derivative of the function. Definition 2.2. For γ > −1 and n = ⌈γ⌉ the Caputo fractional derivative of a function f is defined by for t > 0.
Again, an explicit formula can be given in case where f (t) = t µ for µ > 0 we have that (c) The semi group property holds for fractional integrals: for any f ∈ C[0, T ] and α, β ≥ 0 From the definition of the fractional derivative we see that the model (1) is a Volterra integro differential equation. To proceed with the analysis we will need a result on the corresponding ordinary differential equation.
Using the definition of the Caputo derivative in (6) gives where n = ⌈γ⌉. Next, define v = u ′′ , i.e., We first consider the case γ ∈ (0, 1), i.e., n = 1. By substituting v into (7) we obtain a Volterra integral equation for v: If γ ∈ (−1, 0), the equation is given by Again we wish to let v = u ′′ but to do so directly we must write the integral term in the above equation as Now we substitute the above equation and (8) into the ODE (9) and we have, This is now a Volterra integral equation with a continuous kernel. Hence, see [5, Theorem 2.1.5] and Theorem A.2 in the appendix, a unique solution v ∈ C[0, T ] exists and consequently also the solution u ∈ C 2 [0, T ] of the original equation. The asymptotic behaviour of u ′′ again follows from Theorem A.2.
Next we show that the time-fractional term satisfies a positivity result that will in turn imply its damping properties. Lemma 2.3. Let f ∈ C([0, T ]; L 2 (Ω)), g ∈ C 1 ([0, T ]; L 2 (Ω)) and let γ ∈ (0, 1). Then the following hold: Proof. The proof of an analogue of (a) for Caputo fractional derivatives is given in [4]. The result can also be deduced from Lemma 1.7.2 in [29] and Lemma 3.1 in [23], which for our case gives Minimizing the kernel gives the required result This inequality but with a different constant is given in [21,Theorem A.1].
A numerical test indicates that the constant derived here is larger than the constant in [21, Theorem A.1] for all ν ∈ (0, 1). The second part of the lemma follows by using the semigroup property, see and then applying part (a).
Proof. We follow the standard proof for the wave equation as described, e.g., in [11,17], with modifications required due to the fractional damping term.
Let {ω k } m k=1 be an orthogonal basis of H 1 0 (Ω) and an orthonormal basis of L 2 (Ω) with corresponding eigenvalues {λ k } m k=1 . We look for u m of the form and for k = 1, 2, . . . , m. As this problem is equivalent to (11) with ∂ t u m , integrating in time, and using Lemma 2.3b we obtain where the energy is given by E(t; v) = 1 2 ∂ t v 2 + 1 2 ∇v 2 . Applying Cauchy Schwarz inequality, the definition of the energy and the Gronwall inequality in the usual way we obtain that the energy is bounded independently of m for a constant C = C(T ).
In the standard way, see [11], we obtain a bound on It remains to bound the term containing the fractional derivative.
Recalling that γ + 1 ∈ (0, 1), then using the definition of the Caputo derivative, see Definition 2.2, Young's convolution inequality and (13) we deduce that Returning to (14) we obtain the bound Hence u m has a subsequence that converges weakly to a u in the following spaces Further, in the usual way u is a weak solution of the original problem (1) and satisfies the initial condition [11]. Finally we will show that u is a unique weak solution. To do so we let f = 0, u 0 = 0 and v 0 = 0 and show that the weak solution is u ≡ 0. For a fixed s > 0, let Testing (1) with v we obtain Then integrating with respect to time, where we used Lemma 2.1b. Now by the definition of v, We now prove the corresponding result for γ ∈ (0, 1).
Proof. The only difficulty in extending the proof of Theorem 2.4 to the case of γ ∈ (0, 1) is the bound on ∂ γ+1 t in (15). To circumvent this problem we note that due to the additional smoothness assumption on f , the solution u m ∈ C 2 [0, T ] of (11) satisfies the time differentiated equation Now testing with ∂ 2 t u m and using Lemma 2.3(a) we obtain the energy bound we have obtained a bound on ∂ 2 t u and consequently using again Young's inequality Using this bound in the proof of Theorem 2.4 we obtain the result. [23] is of a similar form to the Szabo equation with γ ∈ (0, 1); see equation (2.5) in [23].

Remark 1. The fractional Zener wave equation investigated in
With similar approach to ours, the authors prove uniqueness and existence of the Zener model. The solution in [23] is understood in a weaker sense, namely they require lower regularity of the data but do not give bounds on the second derivative in time of the solution.
Remark 2. As is usual in PDEs with time fractional derivatives, we expect a singularity at t = 0 even for smooth data. Namely, due to Theorem 2.2 we expect for smooth f and γ ∈ (0, 1) The right hand side is then If f is to be smooth, we need to match the term t 1−γ by setting The right hand side is then If f is to be smooth, we need to match the term t −γ by setting z 0 = − aγ Γ(3−γ) v 0 . Alternatively, for u to be smooth we would need f to have a singularity of the type t 1−γ for γ ∈ (0, 1) and t −γ for γ ∈ (−1, 0) with further weaker singularities at t = 0.

Fully discrete system
To obtain the fully discrete system we will use a finite element method in space and a combination of leapfrog and BDF2 based convolution quadrature discretization in time. The motivation for using leapfrog is to obtain an explicit scheme, whereas the main motivation for using convolution quadrature are its excellent stability properties ( [19,18]) and the ability to evaluate it very efficiently ( [2,24]). An alternative discretization of the fractional time derivative is the L1 scheme ( [26,16]). This scheme also has the required stability property, i.e., it preserves the positivity property of the fractional derivative; see [26, Lemma 3.1].

Spatial semidiscretization
To discretize in space, we make use of a piecewise linear Galerkin finite element method. Namely let V h be a family of finite dimensional subspaces of H 1 0 (Ω) parametrized by the meshwidth h > 0. We assume that these spaces satisfy the following approximation property for some constant C > 0 independent of h. Furthermore, we assume that an inverse inequality holds for some constant C inv . The semidiscrete problem then reads: We denote projections of the initial data u 0 and v 0 onto V h by u h 0 and v h 0 respectively; the projections will be specified later on.

Time discretization
In time we discretize using the explicit, leapfrog scheme , where t n = nκ for some fixed time step κ > 0, u h n ∈ V h is an approximation of u(t n ) and w n is an approximation of ∂ γ+1 t u(t n ) which we describe next. First of all we define the central difference operator We extend the definition also to continuous functions g ∈ C[0, T ] with a given time derivative at t = 0: i.e.,∂ t g(t) is exact at t = 0, the central difference quotient for t ≥ κ and is the linear interpolant of∂ t g for t ∈ [0, κ]. It remains to discretize the fractional derivative ∂ γ t . The formula should be computable efficiently and should retain the positivity property of the fractional derivative as described in Lemma 2.3. All this is satisfied by convolution quadrature introduced by Lubich in [18], which we describe next.
Convolution quadrature (CQ) is based on an A(θ)-stable linear multistep method. For a continuous function g, the CQ formula for ∂ γ t is given by where and ω j are convolution weights defined below. The term χ γ ω n g(0) above is added to the standard definition of convolution quadrature in order to correct for the fact that we are using CQ to compute Caputo rather than Riemann Liouville fractional derivatives. Next, we define the convolution weights. As central differences are approximations of order 2, we restrict our discussion to second order, BDF2 based CQ. The corresponding convolution weights are then given by We can also define the approximation at intermediate values of t by where we define g(t) ≡ 0 for t < 0. Further, it is clear that definition (18) extends to sequences by for t ∈ (0, T ]. As mentioned above, most of the results in the literature analyse CQ as an approximation to Riemann Liouville derivatives ( [18]). However, as the Caputo and Riemann Liouville derivatives of order γ are equivalent for γ ∈ (−1, 0) and for γ ∈ (0, 1) are equivalent if g(0) = 0, we can deduce from [20, Theorem 2.2] the following result.
Next we add correction terms that integrate lower order terms exactly and do not destroy the convergence for the higher order terms. For γ ∈ (−1, 0) it is sufficient to correct for constant functions, whereas for γ ∈ (0, 1) we will need to correct for linears as well.
Here, the correction terms w j are chosen so that Furthermore, denote w nj := w j (t n ) and hence ∂ γ κ * g(t n ) = n j=0 ω n−j g(t j ) + w n0 g(0) + w n1 g(t 1 ).
Note that by definition ∂ γ κ already contains the correction for constants if γ ∈ (0, 1). In that case the above just adds the correction for linear functions.

Lemma 3.2. Under the conditions of Lemma 3.1 it holds that for
with C > 0 independent of κ ∈ (0,κ).
Using the above shorthand for the CQ approximation, we can write the fully discretized systems as or when including the correction n = 1, . . . , N − 1. The coupling of the two time discretizations is similar to the FEM-BEM coupling in [3]. As both of the above schemes are explicit, we will see during the course of the analysis that the following CFL condition is required It remains to describe the choice of initial data u h 0 , u h 1 . To do this we require the Ritz projection, denoted by R h : H 1 0 (Ω) → V h and defined by (∇R h u, ∇v) = (∇u, ∇v) for all v ∈ V h and the L 2 projection is denoted by P h : We have the approximation property, see, e.g., [15], We then define the initial data by Using the PDE and the fact that ∂ 2 t u is continuous we see that P h ∂ 2 t u(0) ∈ V h is the solution of Note also that by definition The reason for using the mixed approximation ∂ γ κ∂ t instead of a fully CQ approximation ∂ γ+1 κ is to conserve the sign of the damping term.
Proof. The proof follows directly from the frequency domain estimate for v ∈ L 2 (Ω) and the Herglotz theorem ([3, Theorem 2.3]).

Convergence analysis
We start with analyzing the error of the mixed approximation ∂ γ κ∂t when applied to t β . We will first need a simple technical lemma.
Proof. Consider first the case η ≤ 0. Then t η is a decreasing function and If η ≥ 0, then t η is an increasing function and we have instead Hence in both cases the sum is bounded by C max(T η+1 , κ η+1 ) with the constant depending on η.
Proof. Throughout the proof, C > 0 denotes a generic constant allowed to change from one step to another.
For β = 0, B 0 = 0 and for β = 1 or β ≥ 2 we have from (21) with ∂ t g = βt β−1 that Combining all the cases gives the stated result. Remark 3. With a different extension of∂ t g(t) for t ∈ (0, κ), it may be possible to improve the estimate for β ∈ (2, 3). However, as we will see in Remark 4 and Lemma A.3, this estimate is sufficient to obtain an optimal global error.
We now prove the corresponding lemma for the corrected quadrature.
We next investigate the error for general, smooth functions.
Lemma 4.4. There exists a constant C > 0 independent of κ ∈ (0,κ) for some small enoughκ > 0, but depending on γ and T such that Proof. Letp = ⌈γ⌉ + 3 and consider the Taylor expansion of u with the Peano kernel remainder term where (t − τ ) + = max(t − τ, 0) and t ∈ [0, T ]. The problem then reduces to analysing the error for polynomials , with g(t) = t k . Hence applying the result of Lemma 4.2 finishes the proof.
The corresponding result with the corrected quadrature is stated next.
Theorem 4.6. Let u be the solution of (10) and u n ∈ V h , n = 0, . . . , N , the solution of the fully discrete system (24). If u is sufficiently smooth we have that The constants implicit in the error estimate are independent of κ ∈ (0,κ), for sufficiently smallκ.
Proof. Let e h n = R h u(t n ) − u h n be the error which satisfies, 1 κ 2 e h n+1 − 2e h n + e h n−1 , w + ∇e h n , ∇w + a γ ∂ γ κ∂t e h (t n ), w = (δ n + a γ ε n , w) , and Mimicking the continuous case, we test with where E e n = As∂ t e h (0) = 0, see (28), from Lemma 3.3 it follows that Hence, by taking the sum over n = 1, . . . , J, J ≤ N − 1, of (31) we obtain the estimate Using the Cauchy Schwarz inequality we have, It remains to bound δ n , ε n and the initial energy. For a sufficiently smooth u a Taylor expansion shows that To bound ε n we split the error into two parts .

Using the Poincaré inequality and the definition of the Ritz projection we have
A ≤ C ∂ γ κ∂ t u(t n ) − ∂ γ t ∂ t u(t n ) 1 . Hence by Lemma 4.4, we have for γ ∈ (−1, 0) and for γ ∈ (0, 1) By (27) we have that B ≤ Ch 2 ∂ γ t ∂ t u(t n ) 2 . Thus from Lemma 4.1 we have that if γ ∈ (−1, 0) and for γ ∈ (0, 1) It remains to estimate the initial error Recalling the definition of the initial data (28) and that e h n = R h u(t n ) − u h n we have for some ξ ∈ (0, κ) Properties of the L 2 projection tell us that With this inequality and (27) we conclude that As R h is the Ritz projection we have that (∇e h 1 , ∇e h 0 ) = 0 and therefore Combining the result together with the triangle inequality applied to u(t n )− u h n = u(t n ) − R h u(t n ) + e h n and the approximation properties of the Ritz projection (27) gives the result.
We next turn to the corrected scheme.
The constants implicit in the error estimate are independent of κ ∈ (0,κ), for sufficiently smallκ.
Proof. The proof is a modification of its noncorrected counterpart. The variation lies in the use of Lemma 4.5 and more critically Lemma 3.3.
In particular using the same notation we have that For γ ∈ (−1, 0), w n1 = 0 so the proof can proceed in the same way. For γ ∈ (0, 1), let us first investigate the case N = 2, i.e., for sufficiently small κ > 0. Hence, proceeding as before and using the same notaion, but applying Lemma 4.5 we obtain that and hence e h (t 2 ) = O(κ 3 ) + O(κh 2 ).
For n > 1 we proceed as where in the last step we used the result for n = 1.
Proceeding with the proof with this alteration we arrive at the result.
Remark 4. The above two theorems assume a sufficiently smooth solution. However, as we have seen in Remark 2, unless the right hand side is of a special form we expect that u has the singular behaviour at t → 0 of the form t 2−γ+⌈γ⌉ .
To investigate the expected convergence rate, we need to investigate ε n and δ n , where we used the notation from the the proofs of Theorem 4.6 and 4.7. We first consider δ n and again using Newton's generalized binomial theorem, see Lemma A.3, find that δ n ≤ Ct ⌈γ⌉−γ−2 κ 2 . Hence, by Lemma 4.1 For γ ∈ (0, 1), the same arguments imply that κ N n=1 ε n = O(κ 2−γ ) for both the corrected and uncorrected scheme, which is in both cases equal to the error due to the approximation of the second derivative.

Smooth solution
First we consider the problem of approximating solutions to (1) in 1D on the interval Ω = [0, 1] with h = 6κ using the two schemes (24) and (25). We construct the right hand side f so that the exact solution is given by u(x, t) = (sin(24t) + cos(12t)) sin(πx). (38) We measure the error in the following norm (39) to compare with the theoretical results in Theorem 4.6 and 4.7.
In Fig. 1 we see that the numerical experiments agree with the predicted convergence rates for various values of γ, except that for γ = 0.25, Fig. 1c, we seem to obtain a higher than expected convergence rate O(h 2 ) in contrast to the predicted rate O(h 1.75 ). However, by increasing the value of α 0 in (2) the predicted convergence rate becomes visible; see   Next we perform an experiment in 2D on the domain Ω = [−1, 1] × [−1, 1] with h = 10κ. The right hand side is chosen so that the exact solution is given by u(x, t) = (sin(24t) + cos(12t)) sin(πx) sin(πy).
In 2D experiments we use the L 2 error defined as In Fig. 3 we show the convergence of the L 2 error Error = max n u n − u(t n ) . and achieve the expected convergence orders for γ = 0.7 with and without correction terms.
The error norm is again as in (39). The results shown in Fig. 4 generally agree with our claims from Remark 4, and in some cases we actually achieve a higher convergence rate than expected. More specifically, for γ = 0.25 we observe second order convergence with and without correction terms, and when γ = 0.75 we have a rate of O(h 1.35 ) with corrected CQ. By increasing α 0 we would see the expected convergence rates in these two cases, similarly to the adjustment we see in
For t ≥ 2κ