Comp. Appl. Math. https://doi.org/10.1007/s40314-018-0601-y

Theoretical analysis and numerical simulation for a hyperbolic equation with Dirichlet and acoustic boundary conditions Adriano A. Alcântara1 · Haroldo R. Clark2 · Mauro A. Rincon1

Received: 26 October 2017 / Revised: 17 January 2018 / Accepted: 24 February 2018 © SBMAC - Sociedade Brasileira de Matemática Aplicada e Computacional 2018

Abstract This paper is concerned with a theoretical and numerical study for the initialboundary value problem for a linear hyperbolic equation with variable coefficient and acoustic boundary conditions. On the theoretical results, we prove the existence and uniqueness of global solutions, and the uniform stability of the total energy. Numerical simulations using the finite element method associated with the finite difference method are employed, for one-dimensional and two-dimensional cases, to validate the theoretical results. In addition, numerically the uniform decay rate for energy and the order of convergence of the approximate solution are also shown. Keywords Hyperbolic equation · Acoustic boundary conditions · Existence and uniqueness of solutions · Energy decay · Finite element method · Finite difference method Mathematics Subject Classification 35L20 · 65M60 · 65M06

1 Introduction Let Ω ⊂ Rn is an open, bounded and connected set with smooth boundary Γ . Suppose also Γ to be made up of two disjoint parts Γ0 and Γ1 (Γ = Γ0 ∪ Γ1 and Γ0 ∩ Γ1 = ∅), both connected and with positive measure.

Communicated by Luz de Teresa.

B

Adriano A. Alcântara [email protected] Haroldo R. Clark [email protected] Mauro A. Rincon [email protected]

1

Institute of Mathematics, Federal University of Rio de Janeiro, Rio de Janeiro, RJ, Brazil

2

Department of Mathematics, Federal University of Piauí, Parnaíba, PI, Brazil

123

A. A. Alcântara et al.

We want to determine a pair of functions, u : Ω × (0, ∞) → R and δ : Γ1 × (0, ∞) → R solutions for the linear coupled system with Dirichlet and acoustic boundary conditions    u (x, t) − α(t)Δu(x, t) = 0 in Q = Ω × (0, T ),   u(x, t) = 0 on Σ0 = Γ0 × (0, T ),    u (x, t) + f 1 δ  (x, t) + f 2 δ  (x, t) + f 3 δ(x, t) = 0 on Σ1 = Γ1 × (0, T ),   ∂u     ∂ν (x, t) = δ (x, t) + g(x)u (x, t) on Σ1 = Γ1 × (0, T ),   u(x, 0) = u 0 (x), u  (x, 0) = u 1 (x) in Ω,   δ(x, 0) = δ (x), δ  (x, 0) = ∂u 0 (x) − g(x)u (x) on Γ ,  0 1 1 ∂ν

(1)

where α : [0, ∞) → R and {g, f i } : Γ1 → R for i = 1, 2, 3 are given functions, and ν denotes the unit outward normal vector on Γ1 . The prime symbol,  , denotes the partial derivative with respect to t, and all derivatives are in the sense of the distributions of Schwartz. In this work, we present results of existence, uniqueness of solution and asymptotic behaviour for the of global solutions for the initial-boundary value problem (1). In addition, we developed a numerical method and we show some examples of numerical solutions for different types of meshes and the asymptotic decay of the energy, to verify the efficiency and feasibility of the method and verify consistency with the theoretical results, in one-dimensional dimensions and two-dimensional cases. We also show that the numerical convergence rate from the numerical results is consistent with the order of theoretical convergence. As far as we know, the numerical simulation for the model (1) is unknown. Since the Eq. (1)1 has time-dependent coefficient α = α(t), one cannot use semigroup theory to study the solvability of system (1). So, to prove the existence and uniqueness of solutions we employ Faedo–Galerkin–Lions’ method and the energy method, respectively, with no restrictions on the geometry of the set Ω. However, to prove that the energy associated with the problem is uniformly stable, we assume the subsets Γ0 and Γ1 to have a special geometry which exclude sets having smooth connected boundary or simply connected regions. The Dirichlet boundary conditions (1)2 and the acoustic boundary conditions (1)3 and (1)4 appear naturally due to the acoustic wave motion. As Beale and Rosencrans (1974), we suppose that a fluid fills the region Ω and it is at rest except for an acoustic wave motion with potential velocity u = u(x, t), and that each point on the surface Γ1 reacts like a resistive harmonic oscillator to the excess pressure of the acoustic wave. In this case, the normal displacement δ = δ(x, t) of the boundary at time t satisfies Eq. (1)3 . The coupling between the acoustic pressure and the displacement of the boundary, called impenetrability condition, is given by Eq. (1)4 with the presence of the damping term “g(x)u  (x, t)”, in which ∂u  ∂ν (x, t) = δ (x, t) is obtained from the continuity of the velocity at the boundary. Acoustic boundary conditions were discussed in Morse et al. (1969) and its corresponding time-dependent version was first formulated by Beale and Rosencrans (1974). This pioneer work was followed by Beale (1976), where a detailed analysis of the acoustic boundary conditions for the linear hyperbolic equation in bounded domains was given. In Beale (1976), the problem was formulated as an initial value problem in a Hilbert space and semigroup methods were used to solve it. Moreover, the semigroup generator has essential spectrum and there is no uniform rate of decay for solutions of the problem. We also refer to the readers the following related works: Cousin et al. (2004), Frigeri (2010), Frota and Goldstein (2000), Frota and Larkin (2005), Graber (2012), Kobayashi and Tanaka (2008), Límaco et al. (2011), Mugnolo (2006), and the references therein. Recently, in Frota et al. (2011) and Vicente and Frota (2013a, b), the authors studied acoustic boundary

123

Theoretical analysis and numerical simulation…

conditions for non-locally reacting boundaries, and in Silva et al. (2017), the mixed nonlinear acoustic boundary conditions was studied for a non-linear parabolic–hyperbolic problem. In the theoretical aspects the main improvement addressed in this paper is related with stabilization of the energy of problem (1). Looking at some previous articles on this subject, for instance, Clark et al. (1998) and Medeiros and Miranda (1996), it is assumed that α  (t) ≤ 0 holds a. e. in (0, +∞) in order to prove that the energy of the system decays exponentially, here we show that this assumption is disposable as one can see in Sect. 3. In Bailly and Juve (2000), was studied some numerical solutions of acoustic propagation problems using linearized Euler’s equations, in that, such solutions were obtained by using a dispersion-relation-preserving scheme in space, combined with a fourth-order Runge–Kutta algorithm in time. Numerical solution of the acoustic wave equation using Raviart–Thomas elements was presented in Jenkins (2007). The rest of this paper is organized as follows: In Sect. 2, we establish the existence and uniqueness of global solutions of the model (1). In Sect. 3, we establish the analysis of the uniform stability of the energy. In Sect. 4, we present the numerical method from the variational formulation, with the finite element method for spatial discretization and the finite difference method for solving a system of differential equations of second order in time. In Sect. 5, we show some numerical results to validate the numerical method employed and some tables for various meshes showing the order of convergence, as well as, the graphs of numerical solutions.

2 Existence and uniqueness theory Regarding the functional spaces, we will use the standard notation for the usual Lebesgue and Sobolev spaces. In particular, we will need to use the two following spaces:   1. The space V = v ∈ H 1 (Ω); γ0 (v) = 0 a. e. on Γ0 is a Hilbert space with norm induced by the of H 1 (Ω), i.e., for all v ∈ V , n   2  ∂v |v|2V = (x) dx = |∇v|2L 2 (Ω) Ω ∂ xi i=1

is a norm equivalent to the standard norm of H 1 (Ω) and the trace map of order zero γ0 : H 1 (Ω) → H 1/2 (Γ  ) is linear and continuous. 2. The space HΔ (Ω) = v ∈ H 1 (Ω); Δv ∈ L 2 (Ω) with inner product (v, w) HΔ (Ω) = (v, w) H 1 (Ω) + (Δv, Δw) L 2 (Ω) is a Hilbert space and |v|2HΔ (Ω) = |v|2H 1 (Ω) + |Δv|2L 2 (Ω) is a norm equivalent to the standard norm of H 2 (Ω) and the trace map of order one γ1 : HΔ (Ω) → H −1/2 (Γ ) is linear and continuous. Therefore, γ0 (v) = v|Γ

and γ1 (v) =

∂v  ∂ν Γ

for all v ∈ D(Ω).

Thus, the trace of the normal derivative in (1)4 makes precise.

123

A. A. Alcântara et al.

Henceforth, the symbols (·, ·), (·, ·)Γ1 , |·|, |·|Γ1 denote the inner products and corresponding norms of the Hilbert spaces L 2 (Ω) and L 2 (Γ1 ), respectively. The concept of global solution for the initial-boundary value problem (1) is given by the following: Definition 1 A global solution for system (1) is a pair of functions (u, δ) with u : Q → R and δ : Σ1 → R belong to the class u ∈ L ∞ (0, ∞; V ∩ HΔ (Ω)), u  ∈ L ∞ (0, ∞; V ), u  ∈ L ∞ (0, ∞; L 2 (Ω)), γ0 (u  ), γ0 (u  ) ∈ L 2 (0, ∞; L 2 (Γ1 )), δ, δ  , δ  ∈ L ∞ (0, ∞; L 2 (Γ1 )), and satisfy the following equations:  ∞  ∞ [u  ϕ + α∇u · ∇ϕ] dx dt = α[δ  − gu  ]ϕ dx dt, 0 0 Ω Γ1 ∞ [u  + f 1 δ  + f 2 δ  + f 3 δ]ψ dx dt = 0, 0

(2)

(3)

Γ1

for all ϕ ∈ and for all ψ ∈ L 2 (0, T ; L 2 (Γ1 )). Moreover, u and δ satisfy the initial conditions (1)5 and (1)6 . L 2 (0, T ; V )

Assume that α(t) ∈ C 1 ([0, ∞)), α  (t) ∈ L 1 (0, ∞) ∩ L ∞ (0, ∞) and α(t) ≥ α0 > 0; g(x) ∈ C(Γ 1 , R), such that g(x) ≥ g0 > 0; f i (x) > 0, i = 1, 2, 3 and f i ∈ C(Γ 1 , R).

(4)

Theorem 1 Suppose (4) hold. Let u 0 ∈ V ∩ HΔ (Ω), u 1 ∈ V and δ0 ∈ L 2 (Γ1 ) such that ∂u 0 − δ  (x, 0) + gu 1 = 0. Then there exists a unique pair of functions (u, δ) solution of ∂ν problem (1) in the sense of Definition 1. Proof Existence of solutions: Since Eq. (1)1 has time-dependent coefficient a natural method to show the existence of solution for system (1) is the classical Faedo–Galerkin–Lions method. For this, we choose below an appropriate Hilbertian basis that satisfies the feedback condition: (1)4 . Remark 1 There exists an orthonormal basis (w j ) j∈N in L 2 (Ω) which is an orthogonal and complete system in V ∩ HΔ (Ω). In fact, since the norms of HΔ (Ω) and H 2 (Ω) are equivalents the proof of this fact is a direct adaptation of Propositions 1–3 of the reference Medeiros and Miranda (1996). Therefore, Proposition 3 of Medeiros and Miranda (1996) would be read thus. Proposition 1 Let u 0 ∈ V ∩ HΔ (Ω), u 1 ∈ V and δ0 ∈ L 2 (Γ1 ) such that ∂u 0 (x) − δ  (x, 0) − g(x)u 1 (x) = 0. ∂ν Then, for each > 0 there are p, q ∈ V ∩ HΔ (Ω) such that || p − u 0 ||2V ∩HΔ (Ω) < , ||q − u 1 ||2V < and

123

∂p (x) − δ  (x, 0) − g(x)q(x) = 0, ∂ν

Theoretical analysis and numerical simulation…

where p ∈ V ∩ HΔ (Ω) is solution of ⎧ −Δp = −Δu 0 ⎪ ⎨ p=0 ⎪ ⎩ ∂ p = δ  + gq ∂ν Proof First, since δ  (x, 0) =

∂u 0 ∂ν (x) −

in Ω on Γ0 on Γ1 .

g(x)u 1 (x) then

∂p ∂u 0 (x) − (x) + g(x)[u 1 (x) − q(x)] = 0. ∂ν ∂ν

()

Now, according to the ideas contained in Medeiros and Miranda (1996)

∂p ∂u 0 ()

2 p − u 0 2V ∩HΔ (Ω) = |Δp − Δu 0 |2 + = g[u 1 − q] 2H −1/2 (Γ ) −

1 ∂ν ∂ν H −1/2 (Γ1 ) 2 2 ≤ C u 1 − q H 1/2 (Γ ) ≤ C u 1 − q V. 1

Since V ∩ HΔ (Ω) is dense in V then for each > 0 there exists q ∈ V ∩ HΔ (Ω) such that u 1 − q 2V < . Therefore, p − u 0 2V ∩HΔ (Ω) < . 

Thus, from Proposition 1, let (w j ) j∈N be a orthonormal basis in L 2 (Ω) which is an orthogonal and complete system in V ∩ HΔ (Ω) and (z j ) j∈N an orthonormal basis of L 2 (Γ1 ) . For each m ∈ N, set Vm = span {w1 , . . . , wm } and Z m = span {z 1 , . . . , z m }. Associated with the initial data (u 0 , u 1 , δ0 ) we take sequences (u 0m , u 1m , δ0m )m∈N such that u 0m ∈ Vm , for all m ∈ N and u 0m → u 0 in (V ∩ HΔ (Ω)), u 1m ∈ Vm , for all m ∈ N and u 1m → u 1 in V, δ0m ∈ Z m , for all m ∈ N and δ0m → δ0 in L 2 (Γ1 ).

(5)

For each m ∈ N, there exists a real number Tm > 0 and real functions {c j m , d j m } defined on [0, Tm ] for all 1 ≤ j ≤ m which defines u m (x, t) for all x ∈ Ω and t ∈ [0, Tm ] and δm (x, t), for all x ∈ Γ1 and t ∈ [0, Tm ] such that {u m , δm } satisfies the approximated problem    (t) − gu m (t), w)Γ1 = 0, (u m (t), w) + α(t) (∇u m (t), ∇w) − (δm     γ0 (u m (t)) + f 1 δm (t) + f 2 δm (t) + f 3 δm (t), z Γ = 0, 1

u m (x, 0) = u 0m (x), u m (x, 0) = u 1m (x), for all x ∈ Ω, ∂u 0m (x)  (x, 0) = δm (x, 0) = δ0m (x), δm − g(x)u 1 (x), for all x ∈ Γ1 , ∂ν

(6)

for all w ∈ Vm and for all z ∈ Z m . The existence of a approximated solution {u m , δm } is a consequence of the standard ordinary differential equation theory. In what follows, we get appropriate estimates for {u m , δm } in order to pass to the limit as m → ∞.  in (6) and (6) , respectively, we have Estimate I: Taking w = 2u m and z = 2δm 1 2      2(u m (t), u m (t)) + 2α(t) (∇u m (t), ∇u m (t)) − (δm (t), u m (t))Γ1  + (gu  (t), u m (t))Γ1 = 0,   m           2 u m (t), δm (t) Γ + 2 f 1 δm (t), δm (t) Γ + 2 f 2 δm (t), δm (t) Γ 1 1 1    + 2 f 3 δm (t), δm (t) Γ = 0. 1

123

A. A. Alcântara et al.

Substituting the second equation in the first one, we get    2(u m (t), u m (t)) + 2α(t) (∇u m (t), ∇u m (t)) + ( f 1 δm (t), δm (t))Γ1

    +( f 2 δm (t), δm (t))Γ1 + ( f 3 δm (t), δm (t))Γ1 + (gu m (t), u m (t))Γ1 = 0.

(7)

In (7), observing (4) we have 2g0 |u m (t)|2Γ1 ≤ 2(gu m (t), u m (t))Γ1 and α(t) ≥ α0 > 0, and using the rule of derivation of the product of functions, we find  E m (t) + 2α0 | f 2 δm (t)|2Γ1 + 2g0 α0 |u m (t)|2Γ1   1/2  1/2 ≤ |α  (t)|R |∇u m (t)|2 + | f 1 δm (t)|2Γ1 + | f 3 δm (t)|2Γ1 , 1/2

where

  1/2  1/2 E m (t) = |u m (t)|2 + α(t) |∇u m (t)|2 + | f 1 δm (t)|2Γ1 + | f 3 δm (t)|2Γ1 .

(8)

From the definition of E m (t) and assumption (4)1 , we have 2  |α  (t)|R   1/2  E m (t) + 2α0  f 2 δm (t) + 2g0 α0 |u m (t)|2Γ1 ≤ E m (t), ∀ t ∈ [0, Tm ). Γ1 α0

(9)

 (t)|2 ≥ 0 and integrating (9) on [0, T ), we get Dropping 2α0 | f 2 δm m Γ1 1/2

 E m (t) + 2g0 α0 0

t

|u m (s)|2Γ1 ds ≤ E m (0) +



t

0

|α  (s)|R E m (s) ds. α0

From convergence (5), there exists a positive constant C0 such that E m (0) ≤ C0 for all m ∈ N. Therefore,  t  t 1 E m (t) + 2g0 α0 |u m (s)|2Γ1 ds ≤ C0 + |α  (s)|R E m (s) ds. α0 0 0 As α  ∈ L 1 (0, ∞) by the hypothesis (4)1 , and applying Gronwall’s Lemma we find a positive ( α  L 1 (0,∞) )/α0 constant C1 = C0 e such that  t E m (t) + 2g0 α0 |u m (s)|2Γ1 ds ≤ C1 , 0

  for all m ∈ N and t ∈ [0, Tm ). From this, we can extend the approximated solutions u m , δm to the whole interval [0, ∞) and find a positive constant C2 such that  ∞ |u m (s)|2Γ1 ds) ≤ C2 , (10) E m (t) + 2g0 α0 0

for all t ≥ 0 and independent of m. Estimate II: Differentiating the approximated Eqs. (6)1 and (6)2 with respect to t, setting w = 2u m (t) and z = 2δ m (t),         2(u  m (t), u m (t)) + 2α (t) (∇u m (t), ∇u m (t)) − (δm (t) − gu m (t), u m (t))Γ1    (t) − gu m (t), u m (t))Γ1 = 0, + 2α(t) (∇u m (t), ∇u m (t)) − (δm        2 u m (t) + f 1 δm (t) + f 2 δm (t) + f 3 δm (t), δm (t) Γ = 0, 1

123

Theoretical analysis and numerical simulation…

and substituting the second equation above in the first one, it yields   1/2  (t)|2Γ1 Fm (t) + 2α(t) |g 1/2 u m (t)|2Γ1 + | f 2 δm   1/2  1/2  − α  (t) |∇u m (t)|2 + | f 1 δm (t)|2Γ1 + | f 3 δm (t)|2Γ1    (t) − gu m (t), u m (t))Γ1 , = −2α  (t) (∇u m (t), ∇u m (t)) − (δm where

  1/2  1/2  (t)|2Γ1 + | f 3 δm (t)|2Γ1 . Fm (t) = |u m (t)|2 + α(t) |∇u m (t)|2 + | f 1 δm

(11)

(12)



Taking w = 2u m (t) in identity (6)1 and multiplying both sides of the resulting equation by α  (t)/α(t), we obtain   α  (t)   2 −2α  (t) (∇u m (t), ∇u m (t)) − (δ m (t) − gu m (t), u m (t))Γ1 = 2 u (t) . α(t) m Substituting this identity in inequality (11),   1/2  (t)|2Γ1 Fm (t) + 2α(t) |g 1/2 u m (t)|2Γ1 + | f 2 δm   α  (t)   2 1/2  1/2  u m (t) + α  (t) |∇u m (t)|2 + | f 1 δm (t)|2Γ1 + | f 3 δm (t)|2Γ1 . = 2 α(t)

(13)

 (t)|2 ≥ 0 and noting the definition of In (13), observing (4)1,2 , dropping 2α(t)| f 2 δm Γ1 Fm (t), we have 1/2

Fm (t) + 2α0 g0 |u m (t)|2Γ1 ≤

3  |α (t)|Fm (t). α0

(14)

Now, fix T > 0 arbitrarily and integrate inequality (14) over the interval [0, t[ for t ≤ T . Observing that Fm (0) is bounded and thereafter applying the Gronwall inequality, we obtain a constant C3 > 0, independent of m, such that for all 0 ≤ t ≤ T ,  T Fm (t) + 2α0 g0 |u m (t)|2Γ1 dt ≤ C3 . (15) 0

From hypothesis in (4) on f i one has 0 < f 1 = min f 1 (x) and 0 < f 3 = min f 3 (x) . x∈Γ 1

x∈Γ 1

Using this, (10), (15), and taking C4 = min{ f 1 , f 3 } > 0, it yields    (t)|2Γ1 + |δm (t)|2Γ1 ≤ C2 /C4 and |δm (t)|2Γ1 + |δm (t)|2Γ1 ≤ C3 /C4 . |δm

(16)

Passage to the limit in the approximated system: From estimates (10), (15) and (16), we can extract subsequences (u μ )μ∈N and (δμ )μ∈N of the sequences of approximations (u m )m∈N and (δm )m∈N , respectively, such that ∗

∗

u μ  u and u μ  u  , weak star L ∞ (0, T ; V ), ∗

u μ  u  , weak star L ∞ (0, T ; L 2 (Ω)), γ0 (u μ )  γ0 (u  ) and γ0 (u μ )  γ0 (u  ), weak L 2 (0, T ; L 2 (Γ1 )), ∗

∗

(17)

∗

δμ  δ, δμ  δ  and δμ  δ  , weak star L ∞ (0, T ; L 2 (Γ1 )). With these convergence we can take directly the limit on all terms of the approximate system (6)1, 2 and conclude that the pair of functions (u, δ) satisfy (3).

123

A. A. Alcântara et al.

The regularities in (2)2 –(2)5 are gotten from convergence (17), and the (2)1 is also obtained from this convergence and standard elliptic regularity results. Uniqueness of solutions: We derive the uniqueness of solutions through of the usual energy method. Thus, we sketch the central ideas. Suppose that (u, δ) and ( u,  δ ) are solutions of (1). Let v = u −  u and ϕ = δ −  δ . Thus, v and ϕ have the regularities contained in (2) and, moreover,   v  − αΔv = 0 in L ∞ 0, ∞; L 2 (Ω) , v  + f 1 ϕ  + f 2 ϕ  + f 3 ϕ = 0 in L 2 (0, ∞; L 2 (Γ1 )), −1 ∂v − ϕ  + gv  = 0 in L 2 (0, ∞; H 2 (Γ1 )), ∂ν v(x, 0) = 0, v  (x, 0) = 0 and ϕ(x, 0) = 0, and the variational formulation associated is given by (v  (t), w) + α(t)[(∇v(t), ∇w) − (z  (t), w)Γ1 + (gv  (t), w)Γ1 ] = 0, for all w ∈ V,    v (t) + f 1 ϕ  (t) + f 2 ϕ  (t) + f 3 ϕ(t), σ Γ = 0, for all σ ∈ L 2 (Γ1 ). 1

Taking w = 2v  and σ = 2ϕ  , using the hypotheses of (4), and proceeding in a similar way as done in Estimate I, we find Λ (t) + 2α(t)| f 2 ϕ  (t)|2Γ1 + 2α(t)|g 1/2 v  (t)|2Γ1 ≤ 1/2

|α  (t)| Λ(t), α0

where   1/2 1/2 Λ(t) = |v  (t)|2 + α(t) |∇v(t)|2 + | f 1 ϕ  (t)|2Γ1 + | f 3 ϕ(t)|2Γ1 . 

Since 2α(t)| f 2 ϕ  (t)|2Γ1 + 2α(t)|g 1/2 v  (t)|2Γ1 ≥ 0, then Λ (t) ≤ |αα(t)| Λ(t) for all 0 t ≥ 0. Integrating this inequality from 0 to t ≥ 0, using the hypothesis (4)1 and the Gronwall’s Lemma, we obtain 1/2

Λ(t) ≤ e

( α  L 1 (0,∞) )/α0

Λ(0) = 0,

for all t ≥ 0. From this we conclude the uniqueness of solutions.



3 Asymptotic behaviour The aim of this section is to show that the total energy   1  2 1/2 1/2 |u (t)| + α(t) |∇u(t)|2 + | f 1 δ  (t)|2Γ1 + | f 3 δ(t)|2Γ1 , E(t) = 2

(18)

associated with the strong global solutions of system (1) is asymptotically stable as t → +∞, as long as the subsets Γ0 and Γ1 have a special geometry and the feedback condition g in (4)2 has an specific form. Throughout this section we assume that both Γ0 and Γ1 are closed, connected and disjoints   sets (Γ0 ∩ Γ1 ) = ∅ with positive measure, satisfying Γ0 = {x ∈ Γ ; m(x) · ν(x) ≤ 0}

123

and

Γ1 = {x ∈ Γ ; m(x) · ν(x) > 0},

Theoretical analysis and numerical simulation…

where m(x) := x − x 0 for some fixed x 0 ∈ Rn . We also assume the feedback control in (4)2 to be given by g(x) := m(x) · ν(x), for all x ∈ Γ1 ,

(19)

0 < g0 := min g(x).

(20)

where x∈Γ1

f i , for i = 1, 2, 3. Thus, Set 0 < f i := min f i (x) ≤ max f i (x) =:  x∈Γ1

x∈Γ1

0 < f i ≤ f i (x) ≤  fi

⇐⇒

   f (x) fi 1 ≤ i ≤ , for i = 1, 2, 3;  fi fi    1 ≥ f i (x) ≥ f i , for i = 1, 2, 3.    fi fi

(21)

Using (21), we obtain  1/2  2    f δ (t) ≥ f 2  f 1/2 δ  (t)2 . 2 Γ1 Γ1  f1 1

(22)

Taking the scalar product in L 2 (Ω) and L 2 (Γ1 ), respectively, of (1)1 with u  and of (1)2 with δ  , adding the resulting expression, applying Green’s identity and using the assumption (4), we have E  (t) + α(t)| f 2 δ  (t)|2Γ1 + α(t)|g 1/2 u  (t)|2Γ1  |α  (t)|  1/2 1/2 |∇u(t)|2 + | f 1 δ  (t)|2Γ1 + | f 3 δ(t)|2Γ1 . ≤ 2 1/2

(23)

Note that (23) does not guarantee a priori that the energy E is decreasing. However, assuming that α0 |α  (t)|R ≤ ε, for all t ≥ 0, (24) 2 where the positive real constant ε is chosen such that (28) is satisfied. Using conditions (19), (22) and (24) in inequality (23), we obtain f 2 1/2  2 | f δ (t)|Γ1 + α(t)|(m · ν)1/2 u  (t)|2Γ1  f1 1   ε 1/2 1/2 ≤ α(t) |∇u(t)|2 + | f 1 δ  (t)|Γ1 + | f 3 δ(t)|Γ1 . 4

E  (t) + α(t)

(25)

Finally, in addition to hypotheses (4)1 , we suppose 1 =: α0 ≤ α(t) ≤ α1 , a.e. in R+ . σ

(26)

Now we are in conditions to state and prove the main result of this section. Theorem 2 In addition to hypothesis (4) we suppose (19), (24) and (26). Then the energy E defined in (18) satisfy   E(t) ≤ 3E(0) e− ε/3 t , for all t ≥ 0, (27)

123

A. A. Alcântara et al.

where



1 4f2 1 , 0 < ε = min 1, ,  f 1 (1 + 4a3 ) a2 4K 2

 ,

(28)

{ f 1 , f 2 } and {a3 , a4 } are defined in (21) and (38), respectively, and K 2 is the constant defined in (44). Proof Consider ε a positive real number defined in (28) and E ε (t) := E(t) + ερ(t),

(29)

where ρ(t) = 2(u  (t), (m · ∇)u(t)) + (n − 1/2)(u  (t), u(t)) +[2α(t) − 3α1 ]( f 1 δ  (t), δ(t))Γ1 .

(30)

Hence, using (1)1 and (1)2 , we obtain ρ  (t) = 2α(t)(Δu(t), (m · ∇)u(t)) + 2(u  (t), (m · ∇)u  (t)) +(n − 1/2)α(t)(Δu(t), u(t)) + (n − 1/2)|u  (t)|2  −[2α(t) − 3α1 ] (u  (t), δ(t))Γ1 + ( f 2 δ  (t), δ(t))Γ1  1/2 +( f 3 δ(t), δ(t))Γ1 + [2α(t) − 3α1 ]| f 1 δ  (t)|2Γ1 +2α  (t)( f 1 δ  (t), δ(t))Γ1 =: P1 + · · · + P9 .

(31)

Now we want to get an upper bound for the terms on the right-hand side of (31). Step 1. P1 = 2α(t)(Δu(t), (m · ∇)u(t)): In this first step, we use Rellich’s identity, given by   ∂v (m · ∇)v dΓ, 2 (Δv, (m · ∇)v) = (n − 2)|∇v|2 − (m · ν)|∇v|2 d Γ + 2 Γ Γ ∂ν for all v ∈ V ∩ H 2 (Ω). Since Γ = Γ0 ∪ Γ1 [see, for example, Remark 5.1 in Clark et al. (1998) for more details], we obtain  2  ∂u 1 2α(t) (Δu, (m · ∇)u) ≤ α(t)(n − 2)|∇u|2 + m 2 α(t) dΓ, (32) Γ1 m · ν ∂ν where m =m (x 0 ) = max m(x) Rn . From (1)4 and (19), we see that x∈Ω



∂u ∂ν

2

  2 2 ≤ 2δ  R + 2(m · ν)2 u  (t)R .

Plugging in this inequality into (32) and using (20) and (21), we have P1 ≤ α(n − 2)|∇u|2 +

2α m2 f 1 g0

| f 1 δ  |2Γ1 + 2α m 2 |(m · ν)1/2 u  |2Γ1 . 1/2

(33)

Step 2. P2 = 2(u  (t), (m · ∇)u  (t)): From Green’s identity, rule the derivative of the product and assumption (4)1 , we get 2 α P2 ≤ −n|u  |2 + (m · ν)1/2 u  Γ . (34) 1 α0

123

Theoretical analysis and numerical simulation…

Step 3. P3 = (n −1/2)α(t)(Δu(t), u(t)): From Green’s identity and the boundary conditions (1)4 , we obtain   P3 = −(n − 1/2)α|∇u|2 + (n − 1/2)α δ  u d x − (n − 1/2)α (m · ν)u  u dx. Γ1

Γ1

 2 The continuity of the trace mapping give us (m ·ν)1/2 v(t)Γ ≤ (a1 /2)|∇v(t)|2 , for all v ∈ V and a1 > 0. Using this inequality, (20), (21) and Young’s inequality, we can write,   (n − 1/2)α δ  u dx − (n − 1/2)α (m · ν)u  u dx ≤

2(n

Γ1 − 1/2)2 a

f 1 g0

Γ1

1α

| f 1 δ  |2Γ1 + 2a1 (n − 1/2)2 α|(m · ν)1/2 u  |2Γ1 + 1/2

α |∇u|2 . 8

Thus, P3 ≤ −(n − 5/8)α|∇u|2 +

2(n − 1/2)2 a1 α

f 1 g0   2  1/2  2 + 2a1 (n − 1/2) α (m · ν) u .

| f 1 δ  |2Γ1 1/2

(35)

Γ1

Step 4. We now derive an upper bound for the terms P4 + · · · + P9 . • From Cauchy–Schwarz and Young inequalities, (20) and (21), we get   2 1 3α1    −[2α − 3α1 ](u , δ)Γ1 ≤ α + (m · ν)1/2 u   Γ1 σ 2σ α0    1 3α1  1/2 2 + σα +  f3 δ , Γ1 f 3 g0 2 f 3 g0 α0 where σ > 0 is a constant to be conveniently chosen later. • Again, from Cauchy–Schwarz and Young inequalities and (21), we obtain    f2 3 f 2 α1  1/2  2  + −[2α − 3α1 ]( f 2 δ , δ)Γ1 ≤ α  f1 δ  Γ1 σ f1 2σ f 1 α0     σ f2 3σ α1  f 2  1/2 2 +α +  f3 δ . Γ1 f3 2 f 3 α0 • Here, the choice of σ = 1/α0 in (26) is justified. Indeed,  1/2 2  1/2 2 α1  1/2 2 −[2α − 3α1 ] f 3 δ Γ ≤ −2α  f 3 δ Γ + 3α  f 3 δ Γ 1 1 1 α0  1/2 2  1/2 2 = −2α  f 3 δ Γ + 3αα1 σ  f 3 δ Γ . 1

1

• From assumption (26), we have  1/2 2  1/2 2  1/2 2 [2α − 3α1 ] f 1 δ  Γ ≤ [2α − 3α] f 1 δ  Γ = −α  f 1 δ  Γ . 1

1

1

• From Cauchy–Schwarz and Young inequalities, (4)1 , (21) and (24), we get 1  1/2 2 σ  1/2 2 2α  ( f 1 δ  , δ)Γ1 ≤ |α  |  f 1 δ  Γ + |α  |  f 3 δ Γ 1 1 σ f3  1  1/2 2 f 1  1/2 2 ≤ α  f 1 δ  Γ + σ α  f 3 δ Γ . 1 1 σ f3

123

A. A. Alcântara et al.

Therefore,



 2 1 3α1  (m · ν)1/2 u  Γ + 1 σ 2σ α0     1/2 2 f2 3 f 2 α1 1  1/2  2 +α + + f 1 δ Γ − α  f 1 δ  Γ 1 1 σ σ f1 2σ f 1 α0        1/2 2 1 3α1 f2 3α1  f1 f2 +α σ + + + + + 3α1 − 2  f 3 δ Γ . (36) 1 f 3 g0 2 f 3 g0 α0 f3 2 f 3 α0 f3

P4 + · · · + P9 ≤ (n − 1/2)|u  |2 + α

Plugging (33)–(36) into (31), we obtain 1 1 ρ  (t) ≤ − α(t)|∇u(t)|2 − |u  (t)|2 2 2   1/2 2 2 + α(t)a2 (m · ν)1/2 u  (t)Γ + α(t)a3  f 1 δ  (t)Γ 1 1  1/2 2 2   1/2 − α(t) f δ  (t) − α(t) 2 − σ a4  f δ(t) , 1

Γ1

Γ1

3

(37)

where 1 1 3α1 + 2a1 (n − 1/2)2 + + , α0 σ 2σ α0  2 m2 2a1 (n − 1/2)2 f2 3 f 2 α1 1 a3 := + + + + , σ g0 f 1 g0 f 1 σ f1 2σ f 1 α0   1 3α1 f2 3α1  f1 f2 + + + + 3α1 + . a4 := f 3 g0 2 f 3 g0 α0 f3 2 f 3 α0 f3 a2 := 2 m2 +

(38)

Multiplying (37) by ε and adding with (25), we get  1/2 2 ε ε E ε (t) ≤ − |u  (t)|2 − α(t) |∇u(t)|2 − α(t)ε  f 1 δ  (t)Γ 1 2 4   1/2 2 f2 ε −α(t) f 1 δ  (t)Γ − (1 + 4a3 ) 1  4 f1  2 −α(t)(m · ν)1/2 u  (t)Γ [1 − a2 ε] 1    1/2 2 7 − σ a4  f 3 δ(t)Γ . −α(t)ε 1 4

(39)

Using in (39) the definition of ε (see (28)), we have  1/2 2 ε ε E ε (t) ≤ − |u  (t)|2 − α(t) |∇u(t)|2 − α(t)ε  f 1 δ  (t)Γ 1 2  4  1/2 2 7 −α(t)ε − σ a4  f 3 δ(t)Γ . 1 4

(40)

We now make the appropriate choices for σ . If we set p=

1 f 3 g0

+

 f2 f3

+ 3α1 +

 f1 f3

and

q=

3α1 2 f 3 g0

+

3α1 2f3

,

then (38)3 can be rewritten as 0 . Choosing σ such that σ a4 ≤ 1, since σ = 1/α0   a4 = p +q/α (see (26)), we have (1/α0 ) p + q/α0 ≤ 1. As a consequence, α02 − pα0 − q ≥ 0. These

123

Theoretical analysis and numerical simulation…

inequalities imply that α0 either satisfy  p − p 2 + 4q α0 ≤ 2

or

α0 ≥

p+



p 2 + 4q . 2

From hypothesis (4)1 , α0 have to be positive; with this we discard the first possibility above. Therefore, since α0 = 1/σ , then 0<σ ≤

p+



2 p 2 + 4q

, which implies that σ a4 ≤ 1.

Therefore, since σ a4 ≤ 1, one has 7/4 − σ a4 ≥ 3/4, and thus    1/2  1/2 2 2 7 3 − εα(t) − σ a4  f 3 δ(t)Γ ≤ − εα(t) f 3 δ(t)Γ . 1 1 4 4

(41)

Using the definition of E(t) in (18) and (41) in (40), we get ε E  (t) ≤ − E(t). 2

(42)

From identities (29) and (30) it is easy to prove that there are positive constants K 1 and K 2 such that − 2K 2 E(t) ≤ ρ(t) ≤ 2K 2 E(t), where    f1 K 1 = max 1, ; f3    (n − 1/2)2 m C 2 n 2 3α1 + , K1 1 + K 2 = max 1, 1 + , . 2 α0 2α0 2α0 Multiplying (43) by ε and adding E(t) and, moreover, from (28), we have (ε < can write, 1 3 E(t) ≤ E ε (t) ≤ E(t). 2 2

(43)

(44)

1 4K 2 ). So we

(45) 

From (42) and (45), we obtain (27).

4 Numerical simulation The purpose of this section is to develop a numerical method, based on the finite element method and the finite difference method to make the numerical simulations in order to obtain a pair of approximate solutions (u m , δm ) of system (1).

4.1 Finite element method We now present a semi-discrete formulation for problem (1) using the Galerkin’s method to discretize the spatial variable. Let {Th } be a family of polygonalization Th = {K } of Ω, satisfying the minimum angle condition (see, for instance, Ciarlet 2002) and indexed by the parameter h representing the maximum diameter of elements K ∈ Th . For a given integer l ≥ 1, we introduce the finite element spaces Vhl = {qh ∈ C 0 (Ω); qh | K ∈Pl (K ), ∀K ∈ Th } and Z hl := Vhl |Γ1 = {qh |Γ1 ; qh ∈ Vhl },

123

A. A. Alcântara et al.

where Pl (K ) is the set of polynomials on K of degree less than or equal to l, i.e., Vhl and Z hl are the spaces of piecewise continuous polynomial functions of degree l. For the numerical simulation presented in this paper, we consider in (6), Vm = Vhl and Z m = Z hl with l = 1, i.e., linear basis functions ϕi for the one-dimensional and twodimensional cases. We also consider h being the uniform spatial discretization for the linear and quadrangular domains for the one-dimensional and two-dimensional cases, respectively. The discretization problem of system (1) on Vm and Z m consists to find, using the Faedo– Galerkin–Lions method, the functions u m : (0, T ) → Vm and δm : (0, T ) → Z m given by m m   u m (x, t) = ci (t)ϕi (x) and δm (x, t) = di (t)ϕi (x), i=1

i=1

as solutions of approximate problem (6). Substituting the terms into the equation, we obtain: m 

ci (t)(ϕi , ϕ j ) + α(t)

i=1

m 

ci (t)(∇ϕi , ∇ϕ j ) + α(t)

i=1

− α(t)

m 

m 

ci (gϕi , ϕ j )Γ1

i=1

di (t)(ϕi , ϕ j )Γ1 = ( f, ϕ j );

i=1 m 

di (t)( f 1 ϕi , ϕ j )Γ1 +

i=1

+

m  i=1

m 

di (t)( f 2 ϕi , ϕ j )Γ1 +

m 

di (t)( f 3 ϕi , ϕ j )Γ1

i=1

ci (t)(ϕi , ϕ j )Γ1 = (h, ϕ j )Γ1 , for j = 1, . . . , m.

(46)

i=1

Consider the following notation for the matrices and vectors:       A = ϕi (x), ϕ j (x) , K = ∇ϕi (x), ∇ϕ j (x) , E 1 = ϕi (x), ϕ j (x) Γ , 1       H1 = f 1 ϕi (x), ϕ j (x) Γ , H2 = f 2 ϕi (x), ϕ j (x) Γ , H3 = f 3 ϕi (x), ϕ j (x) Γ , 1 1 1       H4 = gϕi (x), ϕ j (x) Γ , F = f (x, t), ϕ j (x) and G = f (x, t), ϕ j (x) Γ . 1

1

Then, from (46), we obtain a second-order system of ordinary differential equations:    Ac (t) + α(t)K c(t) + α(t)H4 c (t) − α(t)E 1 d  (t) = F(t),   H1 d  (t) + H2 d  (t) + H3 d(t) + E 1 c (t) = G(t),   c(0) = u 0m (x), c (0) = u 1m (x), d(0) = δ0m (x), d  (0) = δ  (x, 0), m

(47)

where c(t) = [c1 (t), . . . , cm (t)]T and d(t) = [d1 (t), . . . , dm (t)]T .

4.2 Finite difference method As, in the general case, the analytical solution of (47) is unknown, so this system will be numerically solved by the finite differences method. For this purpose, consider the discrete time tk = kΔt for k = 0, . . . , N , where Δt = T /N and N is the number of discrete points in the interval [0, T ] satisfying 0 = t0 < t1 < · · · < t N = T . Denoting φ(tk ) = φ k and taking the system (47) at discrete time t = tk+1 and t = tk−1 and doing the arithmetic mean, we obtain

123

Theoretical analysis and numerical simulation…

 K  k+1 k+1  H4  A  (k+1) c α c α(tk+1 )c (tk+1 ) + c(k−1) + + α k−1 ck−1 + 2 2 2  E 1  k+1 (k+1)  1  k+1  + α k−1 d (k−1) = + F k−1 ; + α(tk−1 )c (tk−1 ) − α d F 2 2  H2  (k+1)  H3  k+1  H1  (k+1) d d d + d (k−1) + + d (k−1) + + d k−1 2 2 2  1  E 1  (k+1) c (48) + + c(k−1) = G k+1 + G k−1 . 2 2 Using two-order central difference approximation for terms with second derivative, delayed and advanced difference of order two for terms with first derivative, substituting in (48), and doing some algebraic calculations, we have     Δt  k+1 Δt 2 A + (Δt)2 α k+1 K + 3α (3α k+1 − α k−1 H4 ck+1 − 2 2     − α k−1 )E 1 dk+1 = 4 A + 2Δt (α k+1 − α k−1 )H4 ck − 2 A + (Δt)2 α k−1 K     Δt  k+1 − 3α k−1 H4 ck−1 − 2Δt (α k+1 − α k−1 )E 1 d k α 2   Δt k+1 (α + − 3α k−1 )E 1 d k−1 + (Δt)2 (F k+1 + F k−1 ); 2        Δt E 1 ck+1 + 2H1 + Δt H2 + (Δt)2 H3 dk+1 = 4H1 d k + − 2H1    + Δt H2 − (Δt)2 H3 d k−1 + Δt E 1 ck−1 + (Δt)2 (G k+1 + G k−1 ). +

(49)

Thus, we obtain the following coupled block matrix system, where each matrix has order m × m, which yields a system with order 2m × 2m:    k+1    M1 M2 c L1 = , (50) M3 M4 L2 dk+1 where    Δt  k+1 3α M1 = 2 A + (Δt)2 α k+1 K + − α k−1 H4 , 2   Δt k+1 k−1 (3α − α )E 1 , M2 = − 2     M4 = 2H1 + Δt H2 + (Δt)2 H3 , M3 = Δt E 1 , and the independent terms are given by      L 1 := 4 A + 2Δt (α k+1 − α k−1 )H4 ck − 2 A + (Δt)2 α k−1 K ck−1 − 2Δt (α k+1  Δt k+1   (α − α k−1 )E 1 d k + − 3α k−1 )E 1 d k−1 + (Δt)2 (F k+1 + F k−1 ); 2       L 2 := 4H1 d k + − 2H1 + Δt H2 − (Δt)2 H3 d k−1 + Δt E 1 ck−1 + (Δt)2 (G k+1 + G k−1 ).

123

A. A. Alcântara et al.

To initialize the iterative method (50), that is, for k = 0, the terms c (0) and d  (0) are approximated by central difference, maintaining the order of quadratic convergence.

5 Numerical results Since the pair of exact solution for the model (1) is unknown a priori, we introduce two sufficiently regular forces f (x, t) and f (x, t) on the right side of the equations (1)1 and (1)3 , respectively, such that the exact solution (u, δ), subject to the initial and boundary conditions, are known, in order to validate the used method, measuring the error and the convergence order. Furthermore, we return to the original model (1), that is, we consider f and f null, but in this case the exact solution is unknown. After that, we will take as an exact solution the approximate solution with the best refinement achieved in the discretization, both in time and space. From now on, we will call homogeneous model when f and f are null and non-homogeneous when f and f are non-zero. In Wheeler (1973), error estimates were made for the numerical solutions of the heat and wave equations, both linear and decoupled. In both equations, an error of order O(h 2 +(Δt)2 ) was obtained, using the standard norm L ∞ (0, T ; L 2 (Ω)), functions of piecewise linear bases and finite difference approximations of order two in time. Let us show numerically that the numerical method developed in this work, also has order of quadratic convergence in space and time. In this way, we use the norms L ∞ (0, T ; L 2 (Ω)) and L ∞ (0, T ; L 2 (Γ1 )), piecewise linear bases functions and approximations by finite differences of order two in time to perform the numerical simulations. In order to estimate the convergence order, denoted by p, we consider uniform discretizations identical for time and space, denoted by h i , where for each h i , E i is the error associated with the numerical solution, in the norms indicated on Tables 1, 2, 3 and 4. Thus, given h i and h i+1 , with h i+1 = h i /2, the numerical convergence order is given by p = ln( E i / E i+1 )/ln(2), (see Rincon and Quintino 2016). The results of some numerical examples are presented below, for the one-dimensional and two-dimensional cases.

5.1 One-dimensional case For the non-homogeneous and homogeneous models in the one-dimensional case we consider the final time T = 1, Ω = (0, 1), ∂Ω = Γ0 ∪ Γ1 , where Γ0 = {0} and Γ1 = {1}, and the following conditions: u 0 (x) = −0.5x(x − 1), u 1 (x) = 0, δ0 = 1 and δ  (1, 0) = −0.5. f 1 (1) = f 2 (1) = f 3 (1) = g(1) = 1 and α(t) = exp(−t). Tables 1 and 2 show the error relation, when space and time discretizations are equals, and the order of numerical convergence. Non-homogeneous model In this case, choosing conveniently forces f and f , the exact solution pair is known and given by: u(x, t) = −0.5x(x − 1) cos(t) and δ(1, t) = 1 − 0.5 sin(t). Homogeneous model For f = f = 0, we consider how “exact solution” the numerical solution obtained using a sufficiently fine mesh. In this experiment, we used h = 2−10 which corresponds to N x = 1024 elements.

123

Theoretical analysis and numerical simulation… Table 1 Error of (u m , δm ) with order of convergence p

Table 2 Error of (u m , δm ) with order of convergence p

h = Δt

E

2−5 2−6

um L ∞ (0,1;L 2 (0,1))

δ

m E abs = |δ − δm |

pu m

pδ m

0.5358 × 10−5

0.1096 × 10−4

–

–

0.1307 × 10−5

0.0271 × 10−4

2.0348

2.0144

2−7

0.0323 × 10−5

0.0067 × 10−4

2.0167

2.0091

2−8

0.0080 × 10−5

0.0016 × 10−4

2.0082

2.0046

2−9

0.0020 × 10−5

0.0004 × 10−4

2.0040

2.0023

h = Δt

E

m E abs = |δ − δm |

pu m

pδ m

2−5

0.4340 × 10−3

0.8847 × 10−4

–

–

2−6

0.1232 × 10−3

0.2196 × 10−4

1.8160

2.0100

2−7

0.0357 × 10−3

0.0541 × 10−4

1.7838

2.0211

2−8

0.0102 × 10−3

0.0128 × 10−4

1.7986

2.0723

2−9

0.0025 × 10−3

0.0025 × 10−4

1.9880

2.3228

um L ∞ (0,1;L 2 (0,1))

δ

It is important to note that Tables 1 and 2 show that the order of convergence p is quadratic, that is, the error between the numerical solution and the exact solution decays quadratically when we refine the mesh. In Fig. 1 is shown the graph of the pair of the numerical solution (u m , δm ) and the graphs of the errors, for the non-homogeneous model in the one-dimensional case.

5.2 Two-dimensional case For the one-dimensional case the boundary Γ1 is reduced to a single point. In the twodimensional one we can take a region (a side) and have a greater influence of the acoustic. For the non-homogeneous and homogeneous models, we consider the final time T = 1, Ω = (0, 1) × (0, 1) and ∂Ω = Γ0 ∪ Γ1 , where Γ0 = (2) ∪ (3) ∪ (4) and Γ1 = (1), and the following conditions:

Tables 3 and 4 show the error relation when space and time discretizations are equals, and the order of numerical convergence. Non-homogeneous model In this case, choosing conveniently forces f and f , the exact solution pair of model (1) is known and given by: u(x, y, t) = γ x(x − 1)(y − 1)2 exp(t) and δ(x, t) = 2γ x(x − 1)(exp(t) − 1). Homogeneous model For f = f = 0, we consider how “exact solution” the numerical solution obtained using a sufficiently fine mesh. In this experiment we used h = 2−8 which corresponds to N x × N y = 28 × 28 = 65536 elements in space and Nt = 256 in time.

123

A. A. Alcântara et al.

Fig. 1 Pair of numerical solution (u m , δm ), h = Δt = 2−9 , with absolute error and in the norm L ∞ (0, 1; L 2 (Ω)) Table 3 Error of (u m , δm ) with h = Δt and order of convergence p

Table 4 Error of (u m , δm ) with h = Δt and order of convergence p

um L ∞ (0,1;L 2 (Ω))

h = Δt

E

2−3

0.00336993

2−4 2−5

δm L ∞ (0,1;L 2 (Γ1 ))

pu m

pδ m

0.02583410

–

–

0.00076558

0.00644954

2.1380

2.0020

0.00018447

0.00161226

2.0531

2.0001

2−6

0.00004541

0.00040313

2.0220

1.9997

2−7

0.00001127

0.00010079

2.0097

1.9998

pu m

pδ m

um L ∞ (0,1;L 2 (Ω))

E

δm L ∞ (0,1;L 2 (Γ1 ))

h = Δt

E

2−3

0.00867500

0.00357054

–

–

2−4

0.00250516

0.00102138

1.7919

1.8056

2−5

0.00076351

0.00027234

1.7141

1.9070

2−6

0.00023366

0.00006770

1.7082

2.0080

2−7

0.00006237

0.00001383

1.9053

2.2907

E

It is important to observe that Tables 3 and 4 show that the order of convergence is quadratic, that is, p  2. In Fig. 2 is shown the graph of the pair of the numerical solution (u m , δm ) and the graphs of the errors, for the non-homogeneous model in the two-dimensional case.

123

Theoretical analysis and numerical simulation…

Fig. 2 Pair of numerical solution (u m , δm ) with h = Δt = 2−7 , and with discrete norms error of L ∞ (0, 1; L 2 (Ω)) and L ∞ (0, 1; L 2 (Γ1 )), respectively

Fig. 3 Asymptotic energy decay E(t) for the original system (1). a Case 1D, with h = Δt = 2−9 and T = 10. b Case 2D, with h = Δt = 2−7 and T = 5

Decay of energy As shown in Theorem 2, the energy associated with the system (1), defined in (18), has a decay when t → ∞. In Fig. 3, we show the decay of the energy associated with the numerical solution of the model (1), using the same input data for the one-dimensional (1D) and two-dimensional (2D) cases, and using the Gaussian quadrature method for the calculation of numerical integration to obtain the energy at each discrete time tk .

123

A. A. Alcântara et al.

6 Conclusions In this paper, we concentred on the study of the linear coupled system with Dirichlet and acoustic boundary conditions. We established the existence and uniqueness of global solutions via the Faedo–Galerkin–Lions and energy methods, respectively. In addition, we show that the energy decays asymptotically when t → ∞. The numerical simulations of the model made from the numerical methods developed and presented in the one-dimensional and two-dimensional cases, show satisfactory numerical results, and the numerical convergence order in time and space is quadratic. The numerical simulation for energy (18) justifies the theoretical result of the asymptotic decay of energy shown in Theorem 2. We also showed some numerical examples to illustrate the results presented in the theoretical part. Acknowledgements A. A. Alcantara was partially supported by CAPES-Brazil. M. A. Rincon was partially supported by CNPq-Brazil.

References Bailly C, Juve D (2000) Numerical solution of acoustic propagation problems using linearized euler equations. AIAA J 38(1):22–29 Beale JT (1976) Spectral properties of an acoustic boundary condition. Indiana Univ Math J 25(9):895–917 Beale JT, Rosencrans SI (1974) Acoustic boundary conditions. Bull Am Math Soc 80(6):1276–1278 Ciarlet PG (2002) The finite element method for elliptic problems. SIAM, Philadelphia Clark H, Jutuca LSG, Miranda MM (1998) On a mixed problem for a linear coupled system with variable coefficients. Eletron J Differ Equ 1(04):1–20 Cousin AT, Frota CL, Larkin NA (2004) On a system of Klein–Gordon type equations with acoustic boundary conditions. J Math Anal Appl 293(1):293–309 Frigeri S (2010) Attractors for semilinear damped wave equations with an acoustic boundary condition. J Evol Equ 10(1):29–58 Frota C, Medeiros L, Vicente A et al (2011) Wave equation in domains with non-locally reacting boundary. Differ Integral Equ 24(11/12):1001–1020 Frota CL, Goldstein JA (2000) Some nonlinear wave equations with acoustic boundary conditions. J Differ Equ 164(1):92–109 Frota CL, Larkin NA (2005) Uniform stabilization for a hyperbolic equation with acoustic boundary conditions in simple connected domains. In: Cazenave T et al (eds) Contributions to nonlinear analysis. Progress in nonlinear differential equations and their applications, vol 66. Birkhäuser Basel, pp 297–312 Graber PJ (2012) Uniform boundary stabilization of a wave equation with nonlinear acoustic boundary conditions and nonlinear boundary damping. J Evol Equ 12(1):141–164 Jenkins EW (2007) Numerical solution of the acoustic wave equation using Raviart–Thomas elements. J Comput Appl Math 206(1):420–431 Kobayashi Y, Tanaka N (2008) An application of semigroups of locally Lipschitz operators to carrier equations with acoustic boundary conditions. J Math Anal Appl 338(2):852–872 Límaco J, Clark HR, Frota CL, Medeiros LA (2011) On an evolution equation with acoustic boundary conditions. Math Methods Appl Sci 34(16):2047–2059 Medeiros L, Miranda M (1996) On a boundary value problem for wave equations: existence, uniquenessasymptotic behavior. Revista de Matemáticas Aplicadas 17:47–73 Morse PM, Ingard KU, Beyer R (1969) Theoretical acoustics. J Appl Mech 36:382 Mugnolo D (2006) Abstract wave equations with acoustic boundary conditions. Mathematische Nachrichten 279(3):299–318 Rincon M, Quintino N (2016) Numerical analysis and simulation for a nonlinear wave equation. J Comput Appl Math 296:247–264 Silva PB, Clark H, Frota C (2017) On a nonlinear coupled system of thermoelastic type with acoustic boundary conditions. Comput Appl Math 36(1):397–414 Vicente A, Frota C (2013a) Nonlinear wave equation with weak dissipative term in domains with non-locally reacting boundary. Wave Motion 50(2):162–169

123

Theoretical analysis and numerical simulation… Vicente A, Frota C (2013b) On a mixed problem with a nonlinear acoustic boundary condition for a non-locally reacting boundaries. J Math Anal Appl 407(2):328–338 Wheeler MF (1973) l∞ estimates of optimal orders for galerkin methods for one-dimensional second order parabolic and hyperbolic equations. SIAM J Numer Anal 10(5):908–913

123