On the Accuracy and Efficiency of Transient Spectral Element Models for Seismic Wave Problems

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. On the Accuracy and Efficiency of Transient Spectral Element Models for Seismic Wave Problems Mönkölä, Sanna Mönkölä, S. (2016). On the Accuracy and Efficiency of Transient Spectral Element Models for Seismic Wave Problems. Advances in Mathematical Physics, 2016, Article 9431583. https://doi.org/10.1155/2016/9431583 2016 Research Article On the Accuracy and Efficiency of Transient Spectral Element Models for Seismic Wave Problems Sanna Mönkölä Department of Mathematical Information Technology, University of Jyv¨ askyl¨ a, Agora, P.O. Box 35, 40014 University of Jyv¨ askyl¨ a, Jyv¨ askyl¨ a, Finland Correspondence should be addressed to Sanna M¨ onk¨ ol¨ a; sanna.monko[email protected] Received 30 November 2015; Revised 11 April 2016; Accepted 13 April 2016 Academic Editor: Christian Engstrom Copyright © 2016 Sanna M¨ onk¨ ol¨ a. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This study concentrates on transient multiphysical wave problems for simulating seismic waves. The presented models cover the coupling between elastic wave equations in solid structures and acoustic wave equations in fluids. We focus especially on the accuracy and efficiency of the numerical solution based on higher-order discretizations. The spatial discretization is performed by the spectral element method. For time discretization we compare three different schemes. The efficiency of the higher-order time discretization schemes depends on several factors which we discuss by presenting numerical experiments with the fourthorder Runge-Kutta and the fourth-order Adams-Bashforth time-stepping. We generate a synthetic seismogram and demonstrate its function by a numerical simulation. 1. Introduction Several formulations exist for modeling the seismic vibrations as interaction between acoustic and elastic waves. Typically, the displacement is solved in the elastic structure. The fluid can be modeled using finite element formulations basedonfluidpressure,displacement,velocitypotential,or displacement potential [1]. Two approaches, in which the displacement is solved in the elastic structure, predominate in modeling the interaction between acoustic and elastic waves. Expressing the acoustic wave equation by the pressure in the fluid domain leads to a nonsymmetric formulation (see, e.g., [2, 3]), while using the velocity potential results in a symmetric system of equations (see, e.g., [4–7]). Recently, a velocitystrain formulation has been considered by Wilcox et al. [8]. Essentially, the difference between the different models is only in the choice of the variables presenting the wave propagation in different domains. Nevertheless, from the numerical point of view, this choice determines the features of the coupled problem. Hence, it also gives the guideline for using appropriate discretization and solution methods. Discretization methods play a crucial role in the efficiency. The key factor in developing efficient solution methods is the use of high-order approximations without computationally demanding matrix inversions. We attempt to meet these requirements by using the spectral element [9] method (SEM) for space discretization. The SEM was pioneered in the mid 1980s by Patera [10] and Maday and Patera [11], and it combines the geometric flexibility of finite elements with the high accuracy of spectral methods. The method is widely used for simulating seismic waves, in both frequency and time domains (see, e.g., [12, 13]). We have applied it to time-harmonic acoustoelastic equations in [14, 15], while this paper concentrates on the accuracy and efficiency of the time domain solutions. In time domain simulations, the efficiency of the method depends also strongly on the time discretization. The secondorder central finite difference (CD), or leap-frog, scheme isalow-ordermethodwhichiseasytoimplementand thus in general use for discretizing in time domain. This method is used with the SEM by Komatitsch et al. [16], and a modification for utilizing larger time steps in the fluid than in the solid media is presented by Madec et al. [17]. Since the method is only of second-order accuracy, higher-order time discretizations are needed for efficient computer simulations based on higher-order space discretization. The fourth-order central finite differences are not usable at this stage, since the scheme is unconditionally unstable. Instead, a secondorder accurate time discretization scheme can be converted Hindawi Publishing Corporation Advances in Mathematical Physics Volume 2016, Article ID 9431583, 15 pages http://dx.doi.org/10.1155/2016/9431583 2Advances in Mathematical Physics to a fourth-order accurate one, as is done by Tarnow and Simo [18]. This symplectic scheme is compared with the CD method by Nissen-Meyer et al. [19]. The importance of higher-order time discretizations, in the field of seismology, has recently been considered by, for example, De Basabe and Sen [20], Peter et al. [21], and Liu et al. [22] and references therein. Lax and Wendroff [23] introduced a higher-order generalization of the finite difference scheme. The method has been applied for acoustic wave equations to provide the fourthorder accuracy, with the finite difference space discretization by Dablain [24] and with the spectral element space discretization by Cohen and Joly [25]. However, this approach does not fit to a case, in which an absorbing boundary condition is used for truncating the computational domain, unless velocity-stress or velocity-displacement formulation is used [26]. By Kubatko et al. [27], Runge-Kutta time discretization methods were used in conjunction with discontinuous Galerkin (DG) finite element spatial discretization, and Antonietti et al. [28] have applied the scheme for elastic wave propagation problems. The second- and third-order Runge-Kutta methods are presented, with a fourth-order Hamiltonian-based space discretization for acoustic and elastic waves in separate domains by Ma et al. [29] and with ℎ𝑝-adaptive discontinuous Galerkin method for simulating tsunami propagation by Blaise and St-Cyr [30], respectively. Since there were no research results considering the accuracy and efficiency of the numerical simulation of fluid and solid waves, in which the fourth-order Runge-Kutta (RK) time discretization is combined with higher-order space discretization, we have applied it for acoustic [31] and elastic [32] waves. The method is used with the eighth-order finite difference space discretization for solving an elastodynamic problem also by Martin et al. [33]. With respect to the time step Δ𝑡,theCDmethodis second-order accurate, while the RK method is fourthorder accurate. Although the computational effort of the RK method is approximately four times that of the CD scheme at each time step, the results of using higher-order timestepping scheme were promising [31, 32]. Both methods lead to an explicit time-stepping scheme. The drawback is that the schemes need to satisfy the stability condition, which limits the length of the time step. For acoustic waves, Yang et al. [34] show that a fourth-order Runge-Kutta method is more accurate for time discretization than the fourth-order Lax- Wendroff method. In this paper, we apply the CD and the RK methods to multiphysical wave problems. We also present the fourthorder Adams-Bashforth (AB) scheme to discuss the efficiency factors of the higher-order time discretization schemes. In what follows, we first consider the mathematical models and boundary conditions in Section 2. Further, we discretize the modelswithrespecttospaceandtimeinSections2.2and2.3. Numerical examples are considered for the validation of the accuracy of the approximations and for demonstrating the synthetic seismogram in Section 3. The concluding remarks are presented in Section 4. Γes Γes Γi Γes Γef Γef Γef ΩsΩf x1=0 Figure 1: The domain Ωis divided into the solid part Ω𝑠and the fluid part Ω𝑓. 2. Numerical Models In what follows, we concentrate essentially on two multiphysical linear models for simulating the transient propagation of seismic waves. In both models, the domain Ω⊂R2is divided into the solid part Ω𝑠and the fluid part Ω𝑓(see Figure 1). First, we present the model for elastic deformations in Ω𝑠 and use it as a starting point for deriving the different models presenting the propagation of acoustic vibrations in the fluid domain Ω𝑓. Then, we apply the coupling conditions to the models and discretize the problems. We see that, despite the fact the difference between the models is only in the choice of the variable presenting the wave propagation, it affects the structureofthecoupledproblematthediscretelevel.Forthe sake of the efficiency of the solving process and the accuracy of the numerical solution, the numerical methods, relying on the model-specific features, are presented. 2.1. Coupled Problems. The deformation of the elastic structure is modeled by the displacement uin the solid domain Ω𝑠, such that 𝜌𝑠(x)𝜕2u 𝜕𝑡2−∇⋅𝜎(u)=f,(1) where 𝜌𝑠(x)is the density of the structure and f=(f1,f2)𝑇is thesourcefunction.Weconcentrateonisotropicmaterials, forwhichthenumberofessentialmaterialconstantsreduces to the two Lam´ e parameters, 𝜇=𝐸/(2(1+]))and 𝜆= 𝐸]/((1+])(1−2])),where𝐸is the Young modulus describing thestiffnessofthesolidand]is the Poisson ratio presenting the compressibility of the solid as the ratio of lateral to longitudinal strain in a uniaxial tensile stress, such that 0< ]<1/2. Hence, the stress tensor can be presented as 𝜎(u)= 𝜆(∇⋅u)I+2𝜇𝜀(u),whereIis the identity matrix. In general, we assume the medium to be heterogeneous. Thus, the partial derivatives in ∇⋅𝜎apply to 𝜆and 𝜇as well as to the displacement. The speed of pressure waves (P-waves) 𝑐𝑝=√(𝜆(x)+2𝜇(x))/𝜌𝑠(x)and shear waves (S-waves) 𝑐𝑠= √𝜇(x)/𝜌𝑠(x)are presented as functions of the Lam´ e parameters and density. The P-waves move in a compressional motion, while the motion of the S-waves is perpendicular to the direction of wave propagation [35]. In principle, the wave equation in fluid media can be derived as a special case of (1), by taking into account the fact that there are no S-waves in fluids. That is because fluids Advances in Mathematical Physics 3 are inviscid and, thus, have no internal friction and can not support shear stresses. Thus, taking the divergence of (1) leads to modeling acoustic waves in fluid domain Ω𝑓by the pressure field 𝑝(see, e.g., [16, 36]), 1 𝜌𝑓(x)𝑐(x)2𝜕2𝑝 𝜕𝑡2−∇⋅(1 𝜌𝑓(x)∇𝑝)=𝑓, (2) where 𝜌𝑓(x)is the density of fluid, 𝑐(x)=√𝜆(x)/𝜌𝑓(x)is the speed of the wave, 𝑝=−𝜆(x)∇⋅udepends on the vectorvalued displacement uin the fluid domain, and 𝑓=−∇⋅f isthesourcefunction.Byusingthevelocitypotential𝜙,the corresponding model is 1 𝑐(x)2𝜕2𝜙 𝜕𝑡2−∇2𝜙=𝑓𝜙,(3) where 𝑓𝜙isthesourcefunction(see,e.g.,[4,5,17]). The equations need to be completed by the initial and boundary conditions to get a well-posed and physically meaningful problem. In this paper, we use the first-order absorbing boundary condition [37] for truncating the exterior domain on the boundaries Γ𝑒𝑠 and Γ𝑒𝑓. On the interface Γ𝑖betweenfluidandsoliddomains,thenormalcomponents of displacements and forces are balanced. That is presented, depending on the choice of variable in the fluid domain, by 𝜎(u)n𝑠−𝑝n𝑓=0, 𝜌𝑓(x)𝜕2u 𝜕𝑡2⋅n𝑠−𝜕𝑝 𝜕n𝑓=0, (4) or 𝜎(u)n𝑠+𝜌𝑓(x)𝜕𝜙 𝜕𝑡n𝑓=0, 𝜕u 𝜕𝑡⋅n𝑠+𝜕𝜙 𝜕n𝑓=0, (5) where n𝑠and n𝑓are the outward pointing unit normal vector on the boundary Γ𝑖,withrespecttothesolidandfluid domain, respectively. The interface conditions coupling the two domains are needed to model the following fact: when a wave coming from the fluid domain confronts the elastic domain,itisnottotallyreflected,butpartofitpassestothe elastic domain and turns to elastic vibrations. The analogous action is seen when the elastic wave propagates to the fluid, although usually the reflections from fluid to solid are minor when compared with the reflections from solid to fluid. To be more precise, the magnitude of the reflection depends on the difference between the densities and the wave speeds of thematerials.Thus,thereflectionsaremoresignificantfrom comparatively stiff structures than for more flexible obstacles. Fortheweakformulationofthesystemconsistingof(1)- (2) and (4) we introduce the function spaces 𝑉={V∈ 𝐻1(Ω𝑓)}and W={w∈(𝐻1(Ω𝑠))2}. By multiplying (1) with any test function win the space Wand (2) with any test function Vin the space 𝑉,usingGreen’sformula,and substituting the boundary conditions, we get the following weak formulation: find (𝑝,u)satisfying (𝑝(𝑡),u(𝑡))∈(𝑉×W) for any 𝑡∈[0,𝑇]and ∫Ω𝑓1 𝜌𝑓(x)𝑐(x)2𝜕2𝑝 𝜕𝑡2V𝑑𝑥+∫Ω𝑓1 𝜌𝑓(x)∇𝑝⋅∇V𝑑𝑥 −∫Γ𝑖𝜕2u 𝜕𝑡2⋅n𝑠V𝑑𝑠+∫Γ𝑒𝑓 1 𝑐(x)𝜌𝑓(x)𝜕𝑝 𝜕𝑡V𝑑𝑠 =∫Ω𝑓𝑓V𝑑𝑥+∫Γ𝑒𝑓 1 𝜌𝑓(x)𝑦extV𝑑𝑠, ∫Ω𝑠𝜌𝑠(x)𝜕2u 𝜕𝑡2⋅w𝑑𝑥+∫Ω𝑠𝜎(u):𝜀(w)𝑑𝑥−∫Γ𝑖𝑝n𝑓 ⋅w𝑑𝑠+∫Γ𝑒𝑠 𝜌𝑠(x) ⋅(𝑐𝑝(x)𝑛2 𝑠1 +𝑐𝑠(x)𝑛2 𝑠2 𝑛𝑠1𝑛𝑠2 (𝑐𝑝(x)−𝑐𝑠(x)) 𝑛𝑠1𝑛𝑠2 (𝑐𝑝(x)−𝑐𝑠(x))𝑐 𝑝(x)𝑛2 𝑠2 +𝑐𝑠(x)𝑛2 𝑠1 )𝜕u 𝜕𝑡 ⋅w𝑑𝑠=∫Ω𝑠 f⋅w𝑑𝑥+∫Γ𝑒𝑠 gext ⋅w𝑑𝑠, (6) where 𝑦ext and gext =(gext1,gext2)𝑇are the source term acting on the artificial absorbing boundaries. Respectively, we can define the weak formulation for the system of (1), (3), and (5). 2.2. Spatial Discretization. The spectral element method is obtained from the weak formulation of the coupled problem by restricting the problem presented in the infinitedimensional spaces into finite-dimensional subspaces. Thus, in order to produce an approximate solution for the problem, the given domain Ωis discretized into a collection of 𝑁𝑒 quadrilateral elements Ω𝑖,𝑖=1,...,𝑁𝑒,suchthatΩ= ⋃𝑁𝑒 𝑖=1 Ω𝑖. Each element is associated with four nodes. For the discrete formulation, we define the reference element Ωref = [0,1]2and invertible affine mappings G𝑖:Ωref →Ω𝑖such that G𝑖(Ωref )=Ω𝑖.TheboundaryofΩ𝑖is a union of its four edges, which are the images of the four edges of the reference element, obtained by the mapping G𝑖.Respectively, the vertex nodes of Ω𝑖are the images of the four vertices of the reference element, obtained by the mapping G𝑖.These properties are essential for the mappings G𝑖to be sufficiently smooth. Each of 𝑁𝑒elements is individually mapped to the reference element, and we make use of the affine mapping to make transformations from the physical domain to the reference domain and vice versa. The mapping between the reference element and the 𝑖th element is defined such that G𝑖(𝜉,𝜁)=x=(𝑥1,𝑥2)∈Ωi,where𝜉and 𝜁are the Gauss- Lobatto points in the reference element (see, e.g., [9]). After the spatial discretization by the spectral elements, thesemidiscreteformofthecoupledproblemisstatedas M𝜕2u 𝜕𝑡2+S𝜕u 𝜕𝑡+Ku=F,(7) where u∈R 𝑁is the global block vector containing the values of the variables in both fluid and solid domains at time 𝑡 4Advances in Mathematical Physics at the Gauss-Lobatto points of the quadrilateral mesh. The variable which we use in the solid domain is the displacement u(x,𝑡). Thus, the number of degrees of freedom (DOF) in the solid domain, 𝑁𝑠,isinthetwo-dimensionaldomaintwice the number of discretization points in the solid domain 𝑁𝑠. The total number of degrees of freedom, that is, 𝑁, depends on the variable of the fluid domain. If the fluid domain is modeled by using pressure 𝑝(x,𝑡)or velocity potential 𝜙(x,𝑡), we have scalar values at each spatial discretization point, and the number of degrees of freedom in the fluid domain, expressed as 𝑁𝑓,isequaltothenumberofdiscretization points in the fluid domain 𝑁𝑓. In practice, the discrete counterparts of the variables are approximated as linear combinations of the corresponding nodal values and the spectral element basis functions 𝜑𝑖, 𝑖=1,...,𝑁𝑓,and𝜓𝑖,𝑖=1,...,𝑁𝑠, which are higher-order Lagrange interpolation polynomials [9]. In the case of the formulation with pressure and displacement, entries of the 𝑁×𝑁matrices M,S,andKand the right-hand side vector Faregivenbytheformulas M=((M𝑠)11 00 0(M𝑠)22 0 (A𝑓𝑠)1(A𝑓𝑠)2M𝑓), S=((S𝑠)11 (S𝑠)12 0 (S𝑠)21 (S𝑠)22 0 00S𝑓), K=((K𝑠)11 (K𝑠)12 (A𝑠𝑓)1 (K𝑠)21 (K𝑠)22 (A𝑠𝑓)2 00K𝑓,), F=(f𝑠 f𝑓), (8) where the 𝑁𝑓×𝑁𝑓matrix blocks and the 𝑁𝑓-dimensional right-hand side vector corresponding to the fluid domain are (M𝑓)𝑖𝑗 =∫Ω𝑓1 𝜌𝑓(x)𝑐(x)2𝜑𝑖𝜑𝑗𝑑𝑥, (S𝑓)𝑖𝑗 =∫Γ𝑒𝑓 1 𝜌𝑓(x)𝑐(x)𝜑𝑖𝜑𝑗𝑑𝑠, (K𝑓)𝑖𝑗 =∫Ω𝑓1 𝜌𝑓(x)∇𝜑𝑖⋅∇𝜑𝑗𝑑𝑥, (f𝑓)𝑖=∫Ω𝑓𝑓𝜑𝑖𝑑𝑥+∫Γ𝑒𝑓 𝑦ext𝜑𝑖𝑑𝑠, (9) where 𝑖,𝑗=1,...,𝑁𝑓.Respectively,the2𝑁𝑠×2𝑁𝑠block matrices and the 2𝑁𝑠-dimensional vector representing the elastic waves are M𝑠=((M𝑠)11 0 0(M𝑠)22), S𝑠=((S𝑠)11 (S𝑠)12 (S𝑠)21 (S𝑠)22), K𝑠=((K𝑠)11 (K𝑠)12 (K𝑠)21 (K𝑠)22), f𝑠=((f𝑠)1 (f𝑠)2), (10) which have the components ((M𝑠)11)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)𝜓𝑗𝜓𝑖𝑑𝑥, ((M𝑠)22)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)𝜓𝑗𝜓𝑖𝑑𝑥, ((S𝑠)11)𝑖𝑗 =∫Γ𝑒𝑠 𝜌𝑠(x)(𝑐𝑝𝑛2 𝑠1 +𝑐𝑠𝑛2 𝑠2)𝜓𝑗𝜓𝑖𝑑𝑠, ((S𝑠)12)𝑖𝑗 =∫Γ𝑒𝑠 𝜌𝑠(x)(𝑐𝑝−𝑐𝑠)𝑛𝑠1𝑛𝑠2𝜓𝑗𝜓𝑖𝑑𝑠, ((S𝑠)21)𝑖𝑗 =∫Γ𝑒𝑠 𝜌𝑠(x)(𝑐𝑝−𝑐𝑠)𝑛𝑠1𝑛𝑠2𝜓𝑗𝜓𝑖𝑑𝑠, ((S𝑠)22)𝑖𝑗 =∫Γ𝑒𝑠 𝜌𝑠(x)(𝑐𝑝𝑛2 𝑠2 +𝑐𝑠𝑛2 𝑠1)𝜓𝑗𝜓𝑖𝑑𝑠, ((K𝑠)11)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)((𝑐2 𝑝−2𝑐2 𝑠)𝜕𝜓𝑗 𝜕𝑥1𝜕𝜓𝑖 𝜕𝑥1 +2𝑐2 𝑠(𝜕𝜓𝑗 𝜕𝑥1𝜕𝜓𝑖 𝜕𝑥1+12𝜕𝜓𝑗 𝜕𝑥2𝜕𝜓𝑖 𝜕𝑥2))𝑑𝑥, ((K𝑠)12)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)((𝑐2 𝑝−2𝑐2 𝑠)𝜕𝜓𝑗 𝜕𝑥2𝜕𝜓𝑖 𝜕𝑥1 +𝑐2 𝑠𝜕𝜓𝑗 𝜕𝑥1𝜕𝜓𝑖 𝜕𝑥2)𝑑𝑥, ((K𝑠)21)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)((𝑐2 𝑝−2𝑐2 𝑠)𝜕𝜓𝑗 𝜕𝑥1𝜕𝜓𝑖 𝜕𝑥2 +𝑐2 𝑠𝜕𝜓𝑗 𝜕𝑥2𝜕𝜓𝑖 𝜕𝑥1)𝑑𝑥, ((K𝑠)22)𝑖𝑗 =∫Ω𝑠𝜌𝑠(x)((𝑐2 𝑝−2𝑐2 𝑠)𝜕𝜓𝑗 𝜕𝑥2𝜕𝜓𝑖 𝜕𝑥2 Advances in Mathematical Physics 5 +2𝑐2 𝑠(12𝜕𝜓𝑗 𝜕𝑥1𝜕𝜓𝑖 𝜕𝑥1+𝜕𝜓𝑗 𝜕𝑥2𝜕𝜓𝑖 𝜕𝑥2))𝑑𝑥, ((f𝑠)1)𝑖=∫Ω𝑠 f1𝜓𝑖𝑑𝑥+∫Γ𝑒𝑠 gext1𝜓𝑖𝑑𝑠, ((f𝑠)2)𝑖=∫Ω𝑠 f2𝜓𝑖𝑑𝑥∫Γ𝑒𝑠 gext2𝜓𝑖𝑑𝑠, (11) where 𝑖,𝑗=1,...,𝑁𝑠. The matrices arising from the coupling between acoustic and elastic wave equations are A𝑓𝑠 and A𝑠𝑓, for which it holds that ((A𝑓𝑠)1)𝑖𝑗 =−∫Γ𝑖𝑛𝑠1𝜓𝑗𝜑𝑖𝑑𝑠, ((A𝑓𝑠)2)𝑖𝑗 =−∫Γ𝑖𝑛𝑠2𝜓𝑗𝜑𝑖𝑑𝑠, ((A𝑠𝑓)1)𝑖𝑗 =−∫Γ𝑖𝑛𝑓1𝜑𝑗𝜓𝑖𝑑𝑠, ((A𝑠𝑓)2)𝑖𝑗 =−∫Γ𝑖𝑛𝑓2𝜑𝑗𝜓𝑖𝑑𝑠. (12) For A𝑓𝑠,𝑖=1,...,𝑁𝑓and 𝑗=1,...,𝑁𝑠,whereas,forA𝑠𝑓, 𝑖=1,...,𝑁𝑠and 𝑗=1,...,𝑁𝑓. The computation of the elementwise matrices and vectors involves the integration over the elementwise subregions. Evaluating these integrals analytically is usually complicated or even impossible. That is why a numerical integration procedure is used. In practice, we replace the integrals by finite sums, in which we use Gauss-Lobatto weights and nodal points.Thevaluesofthesesumsarecomputedelementby element with the Gauss-Lobatto integration rule. Collocation points are now the nodes of the spectral element. All but one of the shape functions will be zero at a particular node. Thus, for 𝑖 =𝑗,(M𝑓)𝑖𝑗 =0and (M𝑠)𝑖𝑗 =0implying that the matrices M𝑓and M𝑠are diagonal. Furthermore, the matrix Mis a lower triangular block matrix with diagonal blocks. Thus, the inverse of the matrix Mis also a lower triangular block matrix with diagonal blocks, M−1 =( (M𝑠)−1 11 00 0(M𝑠)−1 22 0 −M−1 𝑓(A𝑓𝑠)1(M𝑠)−1 11 −M−1 𝑓(A𝑓𝑠)2(M𝑠)−1 22 M−1 𝑓),(13) and explicit time-stepping with central finite differences requires only matrix-vector multiplications. In practice, the stiffness matrix Kis assembled once at the beginning of the simulation. It is stored by using the compressed column storage including only the nonzero matrix elements. The other options would have been using a mixed spectral element formulation [38]. In the case of the formulation with velocity potential and displacement, the entries of the 𝑁×𝑁matrices M,S,andK and the right-hand side vector Faregivenbytheformulas M=((M𝑠)11 00 0(M𝑠)22 0 00M𝑓), S=((S𝑠)11 (S𝑠)12 (A𝑠𝑓)1 (S𝑠)21 (S𝑠)22 (A𝑠𝑓)2 (A𝑓𝑠)1(A𝑓𝑠)2S𝑓), K=((K𝑠)11 (K𝑠)12 0 (K𝑠)21 (K𝑠)22 0 00K𝑓,), F=(f𝑠 f𝜙𝑓), (14) where the 𝑁𝑓×𝑁𝑓matrix blocks and the 𝑁𝑓-dimensional right-hand side vector corresponding to the fluid domain are (M𝑓)𝑖𝑗 =∫Ω𝑓𝜌𝑓(x) 𝑐(x)2𝜑𝑖𝜑𝑗𝑑𝑥, (S𝑓)𝑖𝑗 =∫Γ𝑒𝑓 𝜌𝑓(x) 𝑐(x)𝜑𝑖𝜑𝑗𝑑𝑠, (K𝑓)𝑖𝑗 =∫Ω𝑓𝜌𝑓(x)∇𝜑𝑖⋅∇𝜑𝑗𝑑𝑥, (f𝜙𝑓)𝑖=∫Ω𝑓𝜌𝑓(x)𝑓𝜙𝜑𝑖𝑑𝑥+∫Γ𝑒𝑓 𝜌𝑓(x)𝑦𝜙ext𝜑𝑖𝑑𝑠, (15) for 𝑖,𝑗=1,...,𝑁𝑓.The2𝑁𝑠×2𝑁𝑠block matrices and the 2𝑁𝑠-dimensional vector representing the elastic waves areexactlythesameasinthenonsymmetricformulation. The matrices arising from the coupling between acoustic and elastic wave equations are A𝑓𝑠 and A𝑠𝑓,forwhichitholdsthat ((A𝑓𝑠)1)𝑖𝑗 =∫Γ𝑖𝜌𝑓(x)𝑛𝑠1𝜓𝑗𝜑𝑖𝑑𝑠, ((A𝑓𝑠)2)𝑖𝑗 =∫Γ𝑖𝜌𝑓(x)𝑛𝑠2𝜓𝑗𝜑𝑖𝑑𝑠, ((A𝑠𝑓)1)𝑖𝑗 =∫Γ𝑖𝜌𝑓(x)𝑛𝑓1𝜑𝑗𝜓𝑖𝑑𝑠, ((A𝑠𝑓)2)𝑖𝑗 =∫Γ𝑖𝜌𝑓(x)𝑛𝑓2𝜑𝑗𝜓𝑖𝑑𝑠. (16) For A𝑓𝑠,𝑖=1,...,𝑁𝑓and 𝑗=1,...,𝑁𝑠,whereas,forA𝑠𝑓, 𝑖=1,...,𝑁𝑠and 𝑗=1,...,𝑁𝑓. Using the vector-valued displacement u(x,𝑡)in the fluid domain doubles the number of degrees of freedom in the fluid domain. In other words, 𝑁𝑓=2𝑁𝑓. Therefore also the 6Advances in Mathematical Physics memory consumption for computing the pure displacementdisplacement interaction is higher than in the case of the other couplings considered. Furthermore, if the formulation with displacement in both domains is considered, the divergence of the variable in the fluid domain is involved in the weak formulation, and we would need a scheme that approximates the functions in 𝐻(div,Ω𝑓)better than the spectral element method does. For instance, Raviart- Thomas finite elements, which are used for acoustic wave equations [39], could be used for this purpose. Furthermore, the Raviart-Thomas elements are not a good choice for discretizing the solid domain. Nevertheless, the coupling conditions should be satisfied at the interface of the two domains. That is, we would need to change the discretization approach in the solid domain as well, to make the degrees of freedom coincide or fulfill the coupling conditions, for example, by using Lagrange multipliers [40]. 2.3. Time Discretization. After dividing the time interval [0,𝑇]into 𝑁time steps, each of size Δ𝑡=𝑇/𝑁, applying the appropriate time discretization into semidiscretized form (7), and taking into account the initial conditions, we obtain the matrix form of the fully discrete state equation. Despite the use of higher-order space discretization, the time-stepping for transient wave equations is usually performed by low-order methods like the central finite difference (CD) scheme. By proceeding this way, in the case of the nonsymmetric formulation,ateachtimestep𝑖we compute first the displacement u𝑖andthenthepressure𝑝𝑖from the equations M𝑠u𝑖+1 −2u𝑖+u𝑖−1 Δ𝑡2+S𝑠u𝑖+1 −u𝑖−1 2Δ𝑡 +K𝑠u𝑖 +A𝑠𝑓𝑝𝑖=f𝑖 𝑠, M𝑓𝑝𝑖+1 −2𝑝𝑖+𝑝𝑖−1 Δ𝑡2+S𝑓𝑝𝑖+1 −𝑝𝑖−1 2Δ𝑡 +K𝑓𝑝𝑖 +A𝑓𝑠 u𝑖+1 −2u𝑖+u𝑖−1 Δ𝑡2=f𝑖 𝑓, (17) with the initial conditions u0=e𝑠0, 𝑝0=e𝑓0, u1−u−1 2Δ𝑡 =e𝑠1, 𝑝1−𝑝−1 2Δ𝑡 =e𝑓1, (18) where u𝑖,𝑝𝑖,f𝑖 𝑠,andf𝑖 𝑓are the vectors u,𝑝,f𝑠,andf𝑓at 𝑡=𝑖Δ𝑡. Because the matrix sums (M𝑠+(Δ𝑡/2)S𝑠)and (M𝑓+(Δ𝑡/2)S𝑓)are diagonal, their inverses are obtained simply by inverting each diagonal element. Respectively, the state equation for the formulation with velocity potential and displacement is M𝑠u𝑖+1 −2u𝑖+u𝑖−1 Δ𝑡2+S𝑠u𝑖+1 −u𝑖−1 2Δ𝑡 +K𝑠u𝑖 +A𝑠𝑓 𝜙𝑖+1 −𝜙𝑖−1 2Δ𝑡 =f𝑖 𝑠, M𝑓𝜙𝑖+1 −2𝜙𝑖+𝜙𝑖−1 Δ𝑡2+S𝑓𝜙𝑖+1 −𝜙𝑖−1 2Δ𝑡 +K𝑓𝜙𝑖 +A𝑓𝑠 u𝑖+1 −u𝑖−1 2Δ𝑡 =f𝑖 𝜙𝑓, u0=e𝑠0,𝜙0=e𝜙𝑓0,u1−u−1 2Δ𝑡 =e𝑠1,𝜙1−𝜙−1 2Δ𝑡 =e𝜙𝑓1, (19) where u𝑖,𝜙𝑖,f𝑖 𝑠,andfi 𝜙𝑓 are the vectors u,𝜙,f𝑠,andf𝜙𝑓 at 𝑡=𝑖Δ𝑡.Inthiscase,u𝑖+1 and 𝜙𝑖+1 need to be solved simultaneously, and we need to invert (M+(Δ𝑡/2)S)at each time step. If the boundaries of the computational domain consist only of horizontal and vertical lines, the coefficient matrix (M+(Δ𝑡/2)S)needed for inversion at each time step can be implemented as a block matrix consisting of diagonal blocks. Then, u𝑖and 𝜙𝑖can be solved simply by using matrixvector multiplications from the formulas u𝑖=(I−D−1 𝑠Δ𝑡2A𝑠𝑓D−1 𝑓Δ𝑡2A𝑓𝑠)−1 ⋅D−1 𝑠(𝑦1−Δ𝑡2A𝑠𝑓D−1 𝑓𝑦2), 𝜙𝑖=D−1 𝑓(𝑦2−Δ𝑡2A𝑓𝑠u𝑖), 𝑖=1,...,𝑁−1, (20) where D𝑠=M𝑠+(Δ𝑡/2)S𝑠and D𝑓=M𝑓+(Δ𝑡/2)S𝑓are diagonal matrices and 𝑦1=(2M𝑠−Δ𝑡2K𝑆)u𝑖+(Δ𝑡2S𝑆−M𝑠)u𝑖−1 +Δ𝑡2f𝑖 𝑠+Δ𝑡2A𝑠𝑓𝜙𝑖−1, 𝑦2=(2M𝑓−Δ𝑡2K𝑓)𝜙𝑖+(Δ𝑡2S𝑓−M𝑓)𝜙𝑖−1 +Δ𝑡2f𝑖 𝜙𝑓 +Δ𝑡2A𝑓𝑠u𝑖−1. (21) Remark 1. Itwouldalsobepossibletouncenteroneofthe two first-order derivatives in time to uncouple the problem. That approach would involve first-order difference approximations and introduce some dissipation. That is why we neglect deeper considerations of the uncentered schemes. Sate equation (7) can be presented as a system of differential equations 𝜕y/𝜕𝑡=𝑓(𝑡,y(𝑡)),wherey=(u,k)𝑇is avectoroftime-steppingvariablesuand k=𝜕u/𝜕𝑡and the function 𝑓(𝑡,𝑦(𝑡))=(k,−M−1(Sk+Ku−F))𝑇.To Advances in Mathematical Physics 7 this modified form, we can apply the fourth-order Runge- Kutta method, which is a Taylor series method. In general, the Taylor series methods keep the errors small, but there is the disadvantage of requiring the evaluation of higher derivatives of the function 𝑓(𝑡,y(𝑡)). The advantage of the Runge-Kutta method is that explicit evaluations of the derivatives of the function 𝑓(𝑡,y(𝑡))are not required, but linear combinations of the values of 𝑓(𝑡,y(𝑡))are used to approximate y(𝑡).Inthe fourth-order Runge-Kutta method, the approximate yat the 𝑖th time step is defined as y𝑖=y𝑖−1 +Δ𝑡6(𝑘1+2𝑘2+2𝑘3+𝑘4), (22) where y𝑖=(u𝑖,𝜕u𝑖/𝜕𝑡)𝑇contains the global block vector u𝑖, including the values of the variables in both the fluid and the solid domain at the 𝑖th time step, and its derivative k𝑖=𝜕u𝑖/𝜕𝑡 at time 𝑡=𝑖Δ𝑡,𝑖=1,...,𝑁.Theinitialconditionisgivenby y0=e=(e0,e1)𝑇,and𝑘𝑗=(𝑘𝑗1,𝑘𝑗2)𝑇,𝑗=1,2,3,4,arethe differential estimates as follows: (𝑘11 𝑘12)=(𝑓1(𝑖Δ𝑡,u𝑖,k𝑖) 𝑓2(𝑖Δ𝑡,u𝑖,k𝑖)), (𝑘21 𝑘22)=(𝑓1(𝑖Δ𝑡+Δ𝑡2,u𝑖+𝑘11 2,k𝑖+𝑘12 2) 𝑓2(𝑖Δ𝑡+Δ𝑡2,u𝑖+𝑘11 2,k𝑖+𝑘12 2)), (𝑘31 𝑘32)=(𝑓1(𝑖Δ𝑡+Δ𝑡2,u𝑖+𝑘21 2k𝑖+𝑘22 2) 𝑓2(𝑖Δ𝑡+Δ𝑡2,u𝑖+𝑘21 2k𝑖+𝑘22 2)), (𝑘41 𝑘42)=(𝑓1(𝑖Δ𝑡+Δ𝑡,u𝑖+𝑘31,k𝑖+𝑘32) 𝑓2(𝑖Δ𝑡+Δ𝑡,u𝑖+𝑘31,k𝑖+𝑘32)). (23) In other words, in order to get differential estimates (23), the function 𝑓is evaluated at each time step four times, and then the successive approximation of yis calculated by formula (22). If the matrix Mis diagonal, as is in the formulation with the velocity potential, the only matrix inversion needed in time-stepping is computed simply by inverting each diagonal element in the matrix M.Thisrequiresonly𝑛floating point operations, which is the number of diagonal elements in the matrix Mand known as the number of degrees of freedom in the space discretization. Since the matrix Scontains only diagonal blocks and coupling terms, the operation countofthematrix-vectorproductSkis of order 𝑛.Inthe matrix-vector multiplication involving the sparse stiffness matrix K, only nonzero matrix entries are multiplied, which requires the order of 𝑟2𝑛operations, where 𝑟is the order of the polynomials in the spectral element basis. Besides these, 2𝑛additions and 3𝑛multiplications are needed for a single evaluation of the function 𝑓. According to (22), the computation of y𝑖needs 14𝑛floating point operations. Thus, thecomputationalcostforeachtimestepofthestateequation is of order 𝑂(𝑟2𝑛)also with the RK time-stepping. Although the computational cost is of the same order for both the CD and the RK time-steppings, the number of floating point operators needed for the RK is nearly four times that of the CD. Next, we make an effort to decrease the computing time and still maintain the high accuracy provided by higher-order time discretizations. For this purpose, we present the fourthorder Adams-Bashforth (AB) method based on approximating the functions 𝑓(𝑡,y(𝑡))by interpolating polynomials. It gives the solution yat the 𝑖th time step as y𝑖=y𝑖−1 +ℎ24(55𝑓((𝑖−1)Δ𝑡,y𝑖−1) −59𝑓((𝑖−2)Δ𝑡,y𝑖−2)+37𝑓((𝑖−3)Δ𝑡,y𝑖−3) −9𝑓((𝑖−4)Δ𝑡,y𝑖−4)), (24) where y𝑖=( u𝑖,𝜕u𝑖/𝜕𝑡)𝑇contains the vector u𝑖and its derivative k𝑖=𝜕u𝑖/𝜕𝑡at time 𝑡=𝑖Δ𝑡,𝑖=4,...,𝑁.Hence, it is a multistep method that requires information at four previous time steps implying that another method is needed for starting the time marching. At this stage, we utilize the fourth-order Runge-Kutta method for computing the startup values y0,y1,y2,andy3to initialize the multistep method. 3. Numerical Examples In this section, the computational accuracy and efficiency of solving seismic wave problems are considered. The focus of Section 3.1 is on different formulations and implementation processes. In Section 3.2, the efficiency is further analysed by using higher-order time discretization. Finally, simulated seismograms and two-dimensional illustrations presenting the seismic wave motion are given in Section 3.3. The numerical experiments presented in this section are carried out on an AMD Opteron 885 processor at 2.6 GHz. 3.1. Comparison between Symmetric and Nonsymmetric Formulations. We illustrate the computational cost of both symmetric and nonsymmetric formulations by solving a time-dependent problem both expressing the acoustic wave equation by the pressure and by using the velocity potential in thefluiddomain.Theright-handsidesandinitialconditions aredefinedtosatisfytheanalyticalsolution𝑝=𝜔𝜌𝑓(x)sin(𝜔⋅ x)cos(𝜔𝑡),𝜙=−sin(𝜔⋅x)sin(𝜔𝑡),andu=(cos(𝜔⋅x/ 𝑐𝑝(x))cos(𝜔𝑡),cos(𝜔⋅x/𝑐𝑠(x))cos(𝜔𝑡))𝑇. The problem is solved in a domain, which consists of the solid part Ω𝑠=[−1,0]×[0,1]and the fluid part Ω𝑓=[0,1]× [0,1](see Figure 1). We use square-element meshes with mesh step size ℎ=0.1, and the element order is increased in bothpartsofthedomainfrom1to5.Themeshesarematching on the coupling interface Γ𝑖set at 𝑥1=0for 𝑥2∈[0,1]. On the other boundaries we have the absorbing boundary conditions. The material parameters in the fluid domain are 𝜌𝑓(x)=1.0and 𝑐(x)=1.0.Inthesoliddomain,weusethe values 𝑐𝑝(x)=6.20,𝑐𝑠(x)=3.12,and𝜌𝑠(x)=2.7.Theangular frequency 𝜔=4𝜋isthesameforbothmedia,andweset 8Advances in Mathematical Physics 1 2 3 4 5 CPU time Order of the polynomial basis Pressure Velocity potential, Lapack LU Velocity potential, SuperLU 102 101 100 (a) CPU time with respect to the element order Pressure Velocity potential, Lapack LU Velocity potential, SuperLU 1 2 3 4 5 Memory consumption Order of the polynomial basis 106 105 (b) Memory consumption with respect to the element order Figure 2: CPU time (in seconds) and memory (in kilobytes) consumed for solving a time-dependent problem with different formulations. The number of time steps is fixed to be 400, and square-element meshes with mesh step size ℎ=0.1are used in both media. the propagation direction (1,0)by the vector 𝜔=(𝜔1,𝜔2)= (1,0)𝜔.Thetimeinterval[0.0,0.5]is divided into 400 steps, each of size Δ𝑡=0.00125, to guarantee the stability condition also for the higher element orders. When the pressure formulation is considered, the time marching involves only matrix-vector multiplications, and actual matrix inversions are not needed. In the implementation of the velocity potential formulation, a “one-shot” method is required for solving the linear system including u𝑖+1 and 𝜙𝑖+1 at each time step 𝑖=1,...,𝑁.Inthatcase, the matrix which is needed to be inverted is stored either as abandmatrixorbyusingthecompressedcolumnstorage including only the nonzero matrix elements. The Lapack LU decomposition routines dgbtrf and dgbtrs usethebandmatrixstoragemode.Withthese routines the memory and CPU time consumption increases rapidly when the element order is increased (see Figure 2). The SuperLU library routines perform an LU decomposition with partial pivoting, and the triangular system solves through forward and backward substitution. At this stage, we alsoutilizethesparsityofthematrixbyusingthecompressed columnstoragemode.Sinceallthenonzeroelementsare not near the diagonal of the matrix, the compressed column storagemoderequireslessstorageandcomputingoperations than the band storage mode. That is why the SuperLU library gives a less demanding procedure for solving the linear system than the Lapack library. This is because the matrices arising from the space discretization and including the coupling terms have, in general, a sufficiently large bandwidth. Thus, a remarkably larger amount of memory is needed for storing the sparse coefficient matrix in the band matrix form used in conjunction with Lapack routines than in the compressed column storage utilized with the SuperLU. Consequently, less time is needed when fewer elements are employed in the solution procedure with the SuperLU. The performance of the Lapack routines could be improved, for instance, by using an appropriate node numbering of the mesh. We conclude that the computational efforts are of the same order of magnitude whether the problem in the fluid domain is solved with respect to pressure or whether we use velocity potential formulation in conjunction with the linear solver provided by the SuperLU library. The comparison between numerical and analytical solution shows that in both media the accuracy improves when the element order grows until a certain error level is reached (see Figure 3). This error level, shown as a horizontal line in Figure 3, reflects the error level of time discretization. Since we use a fixed number of time steps at all element orders, the error of time discretization becomes dominant for higherorder elements. Hence, finer time steps or higher-order time discretizations are needed in conjunction with higher-order elements. In principle, both symmetric and nonsymmetric formulation should lead to the same order of accuracy at each element order. The results obtained in the solid domain and depicted in Figure 3(a) are perfectly in balance with this hypothesis. However, in Figure 3(b) we see that in the fluid domain higher accuracy is obtained with the velocity potential than with the pressure formulation at each element order. In the case of pressure formulation, the interface condition is derived by differentiating the displacement component twice with respect to time. From the physical point of view, some information is lost in that procedure, and linear growth of the displacement with respect to time is not eliminated at the interface. However, the solutions of the pressure formulation converge towards the analytical solution. Although the pressure formulation gives less accurate results than the Advances in Mathematical Physics 15 [28] P. F. Antonietti, I. Mazzieri, A. Quarteroni, and F. Rapetti, “High order space-time discretization for elastic wave propagation problems,” in Spectral and High Order Methods for Partial Differential Equations—ICOSAHOM 2012: Selected Papers from the ICOSAHOM Conference, June 25–29, 2012, Gammarth, Tunisia, M. Aza¨ ıez, H. El Fekih, and J. S. Hesthaven, Eds., vol. 95 of Lecture Notes in Computational Science and Engineering,pp.87– 97, Springer, Berlin, Germany, 2014. [29] X. Ma, D. Yang, and F. Liu, “A nearly analytic symplectically partitioned Runge-Kutta method for 2-D seismic wave equations,” Geophysical Journal International,vol.187,no.1,pp.480–496, 2011. [30] S. Blaise and A. St-Cyr, “A dynamic hp-Adaptive discontinuous Galerkin method for shallow-water flows on the sphere with application to a global tsunami simulation,” Monthly Weather Review,vol.140,no.3,pp.978–996,2012. [31] E. Heikkola, S. M¨ onk¨ ol¨ a, A. Pennanen, and T. Rossi, “Controllability method for the Helmholtz equation with higher-order discretizations,” Journal of Computational Physics,vol.225,no. 2, pp. 1553–1576, 2007. [32] S. M¨ onk¨ ol¨ a, E. Heikkola, A. Pennanen, and T. Rossi, “Timeharmonic elasticity with controllability and higher order discretization methods,” Journal of Computational Physics,vol.227, no. 11, pp. 5513–5534, 2008. [33] R. Martin, D. Komatitsch, S. D. Gedney, and E. Bruthiaux, “A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using Auxiliary Differential Equations (ADE-PML),” Computer Modeling in Engineering & Sciences,vol.56,no.1,pp.17–42,2010. [34] D. Yang, X. Ma, S. Chen, and M. Wang, “A fourth-order Runge– Kutta method with low numerical dispersion for simulating 3D wave propagation,” in WavesinFluidsandSolids,R.P.Vila,Ed., InTech, 2011. [35] C. P. A. Wapenaar and A. J. Berkhout, ElasticWave Field Extrapolation: Redatuming of Single and Multi-Component Seismic Data, Elsevier, Amsterdam, The Netherlands, 1989. [36] X. Wang and K.-J. Bathe, “Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems,” International Journal for Numerical Methods in Engineering,vol.40,no.11,pp.2001–2017,1997. [37] B. Engquist and A. Majda, “Radiation boundary conditions foracousticandelasticwavecalculations,”Communications on Pure and Applied Mathematics,vol.32,no.3,pp.313–357,1979. [38] G. Cohen and S. Fauqueux, “Mixed spectral finite elements for the linear elasticity system in unbounded domains,” SIAM Journal on Scientific Computing, vol. 26, no. 3, pp. 864–884, 2005. [39] S. K¨ ahk¨ onen, R. Glowinski, T. Rossi, and R. A. M¨ akinen, “Solution of time-periodic wave equation using mixed finite elements and controllability techniques,” Journal of Computational Acoustics,vol.19,no.4,pp.335–352,2011. [40] A. Berm´ udez, P. Gamallo, P. Hervella-Nieto, R. Rodr´ ıguez, and D. Santamarina, “Fluid-structure acoustic interaction,” in Computational Acoustics of Noise Propagation in Fluids. Finite and Boundary Element Methods,pp.253–286,Springer,2008. [41] N. Ricker, “The form and laws of propagation of seismic wavelets,” Geophysics,vol.18,no.1,pp.10–40,1953. [42] T. Vdovina, S. E. Minkoff, and S. M. Griffith, “A two-scale solution algorithm for the elastic wave equation,” SIAM Journal on Scientific Computing,vol.31,no.5,pp.3356–3386,2009. Submit your manuscripts at http://www.hindawi.com Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Mathematics Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Mathematical Problems in Engineering Hindawi Publishing Corporation http://www.hindawi.com Differential Equations International Journal of Volume 2014 Applied Mathematics Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Probability and Statistics Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Mathematical Physics Advances in Complex Analysis Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Optimization Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Combinatorics Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 International Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Operations Research Advances in Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Function Spaces Abstract and Applied Analysis Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 International Journal of Mathematics and Mathematical Sciences Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 The Scientific World Journal Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Algebra Discrete Dynamics in Nature and Society Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Decision Sciences Advances in Discrete Mathematics Journal of Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 Stochastic Analysis International Journal of