Wavelet Based Simulation of Elastic Wave Propagation

One of the multiresolution-based schemes reinforced with mathematical background is the wavelet theory. Development of this theory is simultaneously done by scientists, mathematicians and engineers [1]. Wavelets can detect different local features of data; the properties that locally separated in different resolutions. Wavelets can efficiently distinguish overall smooth variation of a solution from locally high transient ones separated in different resolutions. This multiresolution feature has been interested many researchers, especially ones in the numerical simulation of PDEs [1]. Wavelet based methods are efficient in problems containing very fine and sharp transitions in limited zones of a computation domain having an overall smooth structure. In brief, the most performance of such multiresolution-based methods is obtained in systems containg several length scales.


Introduction
Multiresolution-based studying has rapidly been developed in many branches of science and engineering; this approach allows one to investigate a problem in different resolutions, simultaneously. Some of such problems are: signal & image processing; computer aided geometric design; diverse areas of applied mathematical modeling; and numerical analysis.
One of the multiresolution-based schemes reinforced with mathematical background is the wavelet theory. Development of this theory is simultaneously done by scientists, mathematicians and engineers [1]. Wavelets can detect different local features of data; the properties that locally separated in different resolutions. Wavelets can efficiently distinguish overall smooth variation of a solution from locally high transient ones separated in different resolutions. This multiresolution feature has been interested many researchers, especially ones in the numerical simulation of PDEs [1]. Wavelet based methods are efficient in problems containing very fine and sharp transitions in limited zones of a computation domain having an overall smooth structure. In brief, the most performance of such multiresolution-based methods is obtained in systems containg several length scales.
Regarding wavelet-based simulation of PDEs, two different general approaches have been developed; they are: 1) projection methods, 2) non-projection ones.
In the projection schemes, in general, the wavelet functions are used as solution basis functions. There, all of the computations are performed in the wavelet spaces; the results are finally re-projected to the physical space [2,3]. In non-projection schemes, the wavelets are only used as a tool to detect high-gradient zones; once these regions are captured, the other common resolving schemes (e.g., finite difference or finite volume method) are employed to simulate considered problems. In this approach, all computations are completely done in the physical domain, and thereby the corresponding algorithms are straightforward and conceptually simple [4]. There are some other schemes that incorporate these two general approaches. They use wavelets as basis functions in a wavelet-based adapted grid points, e.g. [5].
The advantages of the wavelet-based projection methods are: 1. Wavelets provide an optimal basis set; it can be improved in a systematic way. To improve an approximation, wavelet functions can locally be added; such improvements do not lead to numerical instability [3]. 2. Most of the kernels (operators) have sparse representation in the wavelet spaces and therefore speed of solutions is high. The band width of the sparse operators can also be reduced by considering a pre-defined accuracy. This leads to inherent adaptation which no longer needs to grid adaptation [6][7][8].The matrix coefficients can easily be computed considering wavelet spaces relationship [6][7][8]. 3. The coupling of different resolution levels is easy [3]. The coupling coefficients can easily be evaluated considering multiresolution feature of the wavelet spaces [6][7][8]. 4. Different resolutions can be used in different zones of the computation domain. 5. The numerical effort has a linear relationship with system size [3]. In the wavelet system, the fast algorithms were developed [9]. Another considerable property of the wavelet transform is its number of effective coefficients: it is much smaller than data size, itself (in spit of the Fourier transform). These two features leads to fast and accurate resolving algorithms.
Regarding non-projection approaches, the common method is to study a problem in accordance with the solution variation; i.e., using different accuracy in different computational domains. In this method more grid points are concentrated around highgradient zones to detect high variations, the adaptive simulation. In this case, only the important physics of a problem are precisely studied, a cost-effective modeling. Once the grid is adapted, the solution is obtained by some other common schemes, (e.g., the finite difference [4,[30][31][32][33][34][35][36][37][38], or finite volume [39][40][41][42][43] method) in the physical space. The wavelet coefficients of considerable values concentrate in the vicinity of high-gradient zones. The coefficients have a one to one correspondence with their spatial grid points, and thereby, by considering points of considerable coefficient values, the grid can be adapted. For this purpose, the points of small enough coefficient values are omitted from the computing grid. In these grid-based adaptive schemes, the degrees of freedom are considered as point values in the physical space; this feature leads to a straightforward and easy method. In some cases the two approaches, projection and non-projection ones, are incorporated; e.g., adaptive collocation methods [5,[23][24][25][26], and adaptive Galerkin ones [28,29].
There is also some other approaches using wavelets only to detect local feature locations, without grid adaptation. In one approach, spurious oscillation locations are captured by the wavelets; thereafter, the oscillations are locally filtered out by a post-processing step [44][45].
The filtering can be done by the conjugate filtering method only in the detected points. In the other approach, to control spurious oscillations the spectral viscosity is locally added in high-gradient/dicontinuous regions; such zones are detcted by the wavelets. This approach is suitable for simulation of hyperbolic systems containing discontinuous solutions. There, artificial diffusion is locally added only in high-frequency components [46]. In these two approaches, the wavelet transforms are used as a tool to detect highly non-uniform localized spatial behaviors and corresponding zones.
The two aforementioned general wavelet based outlooks, projection and non-projection ones, have successfully been implemented for simulation of stress wave propagation problems. The wavelet-based projection methods were successfully used for simulation of wave propagation problems in infinite and semi-infinite medias [12,[47][48][49][50][51][52]. Another important usage is wave propagation in structural engineering elements; e.g. wave propagation in the nano-composites [53]. The non-projection methods were also employed for wave-propagation problems, one can refer to [35][36][37][38].
In brief, it should be mentioned that other powerful and common methods exist for simulation of wave propagation problems for engineering problems; some of which are: the finite difference and finite element schemes, e.g. [54,55]. These methods are precisely studied and relevant numerical strength and drawbacks are investigated. Regarding these schemes, some of important numerical features are: 1) source of numerical errors: truncation and roundoff errors, [56]; 2) effect of grid/element irregularities on truncation error and corresponding dissipation and dispersion phenomena [57]; 3) internal reflections from grids/element faces [58][59][60][61][62][63][64][65]; 4) the inherent dissipation property [66,67]. These features lead in general to numerical (artificial) dissipation and dispersion phenomena. In general to control these two numerical drawbacks in wave propagation problems, it is desirable to refine spatio-temporal discritizations [68]. Considering the spatial domain, this can effectively be done by the wavelet theory. In the wavelet-based projection methods, the inherent adaptation is used, while in the non-projection ones, the multiresolution-based grid adaptation is utilized.
This chapter is organized as follows. In section 2, the wavelet-based projction method will be survived. This section includes: 1) a very brief explanation of main concept of multiresolution analysis; 2) in brief review of wavelet-based projection method for solution of PDEs and computation of the spatial derivatives; 3) the issues related to a 2D wave propagation example. In section 3, the wavelet-based non-projection ones will be presented.
It includes: 1) wavelet-based grid adaptation scheme with interpolating wavelets; 2) solution algorithm; 3) smoothing splines; 4) an example: wave propagation in a two layered media. This chapter ends with a brief conclusion about the presented wavelet-based approaches.

Wavelet based projection method in wave propagation problem
In the wavelet based projection methods, the wavelets are used as basis functions in numerical simulation of wave equations. This section has following sub-sections: multiresolution analysis and wavelets; representation of operators in the wavelet spaces; the semi group time integration methods; a SH wave propagation problem.

Multiresolution analysis and wavelet basis
In this subsection, wavelet-based multiresolution analysis and wavelet construction methods will be survived.

Multiresolution analysis
A function or a signal, in general, can be viewed as a set of a smooth background with low frequency component (approximation one) and local fluctuations (local details) of variant high frequency terms. The word "multiresolution" refers to the simultaneous presence of different resolutions in data. In the multiresolution analysis (MRA), the space of functions that belong to square integrable space, 2 ( ) L  , are decomposed as a sequence of detail subspaces, denoted by   k w , and an approximation subspace, indicated with j v . The approximation of ( ) f t at resolution level j , ( ) f t , is in j v and the details ( ) k d t are in k w (detail sub-spaces of level k ). The corresponding scale of resolution level j is usually chosen to be of order 2 j  [69,70]. In orthogonal wavelet systems, the multiresolution analysis of 2 ( ) L  is nested sequences of the subspaces   j v such that: i.
Exists a function ( ) t  , called the scaling function such that set   The sub-space j v denotes the space spanned by family The second equation implies that the scaling function ( ( ) t  ) has unit energy and therefore by multiplying it with data, the energy of signals do not alter. The dilated and shifted version of the scale function, , ( ) j k t  is usually normalized with the coefficient / 2 2 j to preserve the energy conservation concept; namely, The subspace j w itself is spanned by a dilated and shifted wavelet function is usually referred as the mother wavelet. The wavelet function, , ( ) j k t  is localized both in time (or space) and frequency (scale); it oscillates in such a way that its average to be zero, i.e.: , ( ) 0 j k t dt any function in 0 v or 0 w can be expanded in terms of the basis function of 1 v , i.e.: These important equations are known as: dilation equations, refinement equations or two-scale relationships [69][70][71][72]. The k h and k h  are called filter coefficients, and can be obtained, in general, by the following relationships: where N is length of the scaling coefficient filters, As mentioned before, the multiresolution decomposition of 2 ( ) L  leads to a set of subspaces with different resolution levels; i.e., 2 In this regard, by can be expanded as: This equation shows that the function ( ) f x is converted into an overall smooth approximation (the first parenthesis), and a series of local fluctuating (high-frequency details) of different resolutions (the second parenthesis). In the above equation ( , ) c j k and ( , ) d j k are called the scaling (approximation) and wavelet (detail) coefficients, respectively. These transform coefficients are usually stored in an array as follows: ; this storing style is commonly referred as the standard form.
In the orthogonal wavelet systems, the coefficients ( , ) c j k and ( , ) d j k can be determined by: Fast algorithms were developed to these coefficient evaluations and relevant inverse transform [9,69].
In the following, the multiresolution-based decomposition procedure is qualitatively investigated by an example.

Derivation of filter coefficients
Considering the above mentioned necessary properties of scaling function and other possible assumptions for scaling/wavelet functions, the filter coefficients can be evaluated.  One choice is the necessity that the set function   can exactly reconstruct polynomials of order upto but not greater than p [71,72]. The polynomial can be represented as: By taking the inner product of the wavelet function ( ( ) x  ) with the above equation, other conditions can be obtained; since: As i  coefficients are arbitrary, then it is necessary that each of the above integration to be In case that / 2 p N  , where N is even (for filter coefficients of length N ) the resulted wavelet family is known as the Daubechies wavelets. In this case, for N scaling filter coefficients, N independent equations exist, and unique results can be obtained.
In Figure 2  Other choices can be considered for scaling/wavelet functions construction; some of such assumptions are: imposing vanishing moment conditions for both scaling and wavelet functions (e.g., Coiflet wavelets), obtaining maximum smoothness of functions, interpolating restriction, and/or symmetric condition [69]. To fulfill some of these requirements, the orthogonality requirement can be relaxed and the bi-orthogonal system is used [69]. For numerical purposes, some other requirements can also be considered; for example Dahlke et al. [16] designed a wavelet family which is orthogonal to their derivatives. This feature leads to a completely diagonal projection matrix and thereby a fast solution algorithm [14].

Expressing operators in wavelet spaces
In this subsection, multiresolution analysis of operators will be presented [6][7][8].
Assume T denotes an operator of the following form: The aim is to represent the operator T in the wavelet spaces; this can be done by projection the operator in the wavelet spaces.
The projection of the operator in the approximation space of resolution level j ( j v ) can be represented as: In the same way, the projection of the operator in the detail subspace j w , of resolution j , is: ). This representation is the telescopic form of the operator T .
The schematic shape of the operator T in telescopic (multiresolution) form is presented in Figure 3; this form of representation is known as the Non-Standard form (NS form). In this figure, it is assumed that: min 1 J  , max 4 J  . There, the coefficients i d and i s are the scale and detail coefficients, respectively; these coefficients are obtained from the common discrete wavelet transform of data x . The ˆ& i i d s are the NS form of the wavelet coefficients, and should be converted to standard form by a proper algorithm (will be discussed).
The projection of the operator T in the wavelet space results to set   , , , where j denotes resolution levels. This form is called NS from, since both of the scale ( j s ) and detail ( j d ) coefficients are simultaneously appeared in the formulation, see Figure 3.  s . This is because there is the two-scale relationship between the wavelet (detail) and scale functions; for more details see [6,7]. This fact, leads to a simple and fast algorithm for calculation of j A , j B , j  elements. The NS form of the operator / d dx obtained by the Daubechies wavelet of order 12 ( 12 Db ) is presented in Figure 4; there The aforementioned telescopic representation is for 1D data. For higher dimensions, the extension is straightforward: the method can independently be implemented for each dimension.

The semi-group time discretization schemes
The scheme used here for temporal integration is the semi-group methods [74,75]. These schemes have a considerable stability property: corresponding explicit methods have a stability region similar to typical implicit ones.  where: L and N represent the linear and non-linear terms, respectively; . The initial condition is: in  , and the linear boundary Regarding the standard semi-group method, the solution of the above mentioned equations    ; where I is the identity matrix.

Simulation of 2D SH propagating fronts
The governing equation of the SH scalar wave (anti-plane shear wave) is: is the out-of plane displacement;  and  are shear modules and density, respectively. By defining a linear operator L y , the above mentioned equation can be rewritten as: For using the semi-group temporal integration scheme, a new variable / y y v u t    is introduced and consequently the above equation will be represented as a system of vectors: This system can be rewritten in vector notation as follows [48]: The simplest explicit semi-group time integration scheme is obtained for case in this case the discretized form of the wave equation is: . . .
Wave Propagation Theories and Applications 366

The absorbing boundary conditions: infinite boundaries
The absorbing boundaries are usually used for presenting infinite boundaries. The defect of numerical simulations is occurrence of artificial boundaries which reflect incoming energies to the computation domain. In this study, the absorbing boundary introduced in [76] is used to simulate infinite boundaries, where the absorbing boundary condition is considered explicitly. Therefore, the wave equation is modified by a damping term ( , ). ( , , ) where, ( , ) Q x z is an attenuation factor. This factor is zero in computation domain and increases gradually approaching to the artificial boundaries. Consequently, the waves incoming towards these boundaries are gradually diminished. In general, no absorbing boundary can dissipate all incoming energies, i.e. some small reflections will always remain.
The above mentioned modification, performed for SH wave equation is as follows: And the modified vector form of the equation is:

Free boundaries
There are different approaches for imposing the free boundary conditions in finitedifference methods; some of which are: 1) using equivalent surface forces (explicit implementation) [48]. In this method the equivalent forces will be up-dated in each time step; 2) employing artificial grid points by extending the computing domain (a common method); 3) considering nearly zero properties for continuum domain in simulation of the free ones [77]; in this case the boundary is replaced with an internal one. In the following examples (done by the wavelet based projection method) the third approach will be used. The first method are mostly be used for simple geometries.

Example
In the following, a scalar elastic wave propagation problem will be considered. The results confirm the stability and robustness of the wavelet-based simulations. The plane wave condition is simulated by an initial imposed out-of plane harmonic deformation where corresponding wave number is 64 / 12 k  . For time integration, the simplest form of the semi-group temporal integration method is used. The snapshots of results (displacement ( , , ) y u x z t ) at different time steps are presented in Figure 5; there the light gray circle represents the tunnel. The displacement ( , , 0.0048) y u x z t  is illustrated in Figure 6; the total CPU computation time is 569 sec. for two different uniform grids, this problem is re-simulated by the finite difference method (with accuracy of order 2 in the spatial domain); the grid sizes are 143 143

Wavelet based simulation of second order hyperbolic systems (wave equations)
In this section, wavelet-based grid adaptation method is survived for modeling the second order hyperbolic problems (wave equations). The strategy used here is to remove spurious oscillations directly from adapted grids by a post-processing method. The employed stable smoothing method is the cubic smoothing spline, a kind of the Tikhonov regularization method. This section is devoted to the following subsections: interpolating wavelets and corresponding grid adaptation; relevant algorithm for adaptive simulation of wave equations; smoothing splines definition; a 2D P-SV wave propagation example.

Interpolating wavelets and grid adaptation
In multiresolution analysis, each wavelet coefficient (detail or scale) is uniquely linked to a particular point of underlying grid. This distinctive property is incorporated with compression power of the wavelets and therefore a uniform grid can be adapted by grid reduction technique.
In this method a simple criteria is applied in 1D grid, based on the magnitude of corresponding wavelet coefficients. The existing odd grid points at level j should be removed if corresponding detail coefficients are smaller than predefined threshold (  ); wavelet coefficients and grid points have one-to-one correspondence [69].
In this work, Dubuc-Deslauriers (D-D) interpolating wavelet [69] is used to grid adaptation. The D-D scaling function satisfies the interpolating property and has a compact support [69]. In the case of the D-D wavelets, the grid points correspond to the approximation and detail spaces at resolution j are denoted by j V and j W , respectively. These sets are locations of the wavelet transform coefficients: the ( , ) c j k and ( , ) d j k locations belong to j V and j W , respectively. These locations are: According to this formula, the natural cubic spline interpolation is obtained by 1 p  and the least-squares straight line fit by 0 p  . In 1 p  the interpolating property is vanished while the smoothing property is increased. In the above functional, the errors are measured by summation and the roughness by integral. Therefore, the smoothness and accuracy are obtained simultaneously. In the mentioned equations, the trade-off between smoothness and goodness of fit to the data is controlled by smoothing parameter.
The p should be selected properly, otherwise it leads to over smoothed or under smoothed results. The former are seen in the scheme presented in Reinsch [80], according to Hutchinson-Hoog, [79] and the latter in the scheme offered in Craven-Wahba [83], according to Lee [84].
The smoothness and accuracy in fitting should be incorporated in such a way that the proper adapted grid and accurate solution are obtained simultaneously in adaptive simulations. Hence trial-and-error method is effective in finding appropriate range of p .
This study shows that in { } 1 Smoothing spline, being less sensitive to noise in the data, has optimal properties for estimating the function and derivatives. The error bounds in estimating the function, belonging to Sobolev space, and its derivatives are presented by Ragozin [86]. He showed that the estimation of function and its corresponding derivatives are converged as the interpolating properties and the sampled points are increased [86].
The smoothing splines work satisfactory for irregular data; this is because the method is a kind of Tikhnov regularization scheme [82,87,88].

Numerical example
The following example is to study the effectiveness of the proposed method concerning some phenomena in elastodynamic problems. Regarding using multiresolution-based adaptive algorithm, the simulation of wave-fields can properly be performed in the media especially one has localized sharp transition of physical properties. The example of such media is solid-solid configurations. In fact, to be analyzed by traditional uniform grid-based methods, these media show major challenges. The main assumptions in the presented example are: 1-applying D-D interpolating wavelet of order 3; 2-decomposing the grid (sampled at

Conclusion
Multiresolution based adaptive schemes have successfully been used for simulation of the elastic wave propagation problems. Two general approaches are survived: projection and non-projection ones. In the first case the solution grid in not adapted, while in the second one it is done. the results confirm that the projection method is more stable than the common finite differnce schemes; since in the common methods spurious oscillations develop in numerical solutions. In the wavelet-based grid adaptation method, it is shown that grid points concentrate properly in both high-gradient and transition zones. There, for remedy non-physical oscillations the smoothing splines (a regularization method) are used.