Open access peer-reviewed chapter

# Wavelet Based Simulation of Elastic Wave Propagation

Submitted: December 11th 2011Reviewed: August 4th 2012Published: February 13th 2013

DOI: 10.5772/52097

## 1. 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 . 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 . 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 . 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. .

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 .

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 needsto grid adaptation [6-8].The matrix coefficients can easily be computed considering wavelet spaces relationship [6-8].

3. The coupling of different resolution levels is easy . The coupling coefficients can easily be evaluated considering multiresolution feature of the wavelet spaces [6-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 . In the wavelet system, the fast algorithms were developed . 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.

The wavelet-based projection methods, however, have two major drawbacks: 1) projection of non-linear operators; 2) imposing both boundary conditions and corresponding geometries [4, 8].

The most common wavelet based projection methods are: the telescopic representation of operators in the wavelet spaces [6-8], wavelet-Galerkin [2, 3,10-19], wavelet-Taylor Galerkin [20-22], and collocation methods [5,23-26] (in this approach, the wavelet-based grid adaptation scheme is incorporated with the wavelet-based collocation scheme). Some efforts have been done to impose properly boundary conditions in these methods. Some of which are: 1) wavelets on an interval [11, 27], 2) fictitious boundary conditions [12, 13, 28, 29], 3) reducing edge effects by proper extrapolation of data at the edges , 4) incorporation of boundary conditions with the capacitance matrix method [2, 15].

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 high-gradient 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- 38], or finite volume [39-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 onecorrespondence 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-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 . 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-52]. Another important usage is wave propagation in structural engineering elements; e.g. wave propagation in the nano-composites . The non-projection methods were also employed for wave-propagation problems, one can refer to [35-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, ; 2) effect of grid/element irregularities on truncation error and corresponding dissipation and dispersion phenomena ; 3) internal reflections from grids/element faces [58-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 . 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.

## 2. 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.

### 2.1. Multiresolution analysis and wavelet basis

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

#### 2.1.1. 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, L2(), are decomposed as a sequence of detail subspaces, denoted by {wk}, and an approximation subspace, indicated withvj. The approximation of f(t)at resolution level j, f(t), is in vjand the details dk(t)are in wk(detail sub-spaces of level k). The corresponding scale of resolution level jis usually chosen to be of order 2j[69, 70].In orthogonal wavelet systems, the multiresolution analysis of L2()is nested sequences of the subspaces {vj}such that:

1. ...v1v0v1...L2()

2. v={0},v+=L2

3. f(t)vjf(2t)vj+1

4. f(t)v0f(tk)v0

5. Exists a functionϕ(t), called the scaling function such that set {ϕ(tk)}kZis a basis of v0.

The sub-space vjdenotes the space spanned by family {ϕj,k(t)}, i.e., vj=spank{ϕj,k(t)}¯where ϕj,k(t)=2j/2ϕ(2jtk). The ϕj,k(t)is a scaled and shifted version of the ϕ(t); thereby the function ϕ(t)is known as the father wavelet. The scale functions (ϕj,k(t)) are localized in both spatial (or time) and frequency (scale) spaces. The function ϕ(t)is usually designed so that:ϕ(t)dt=1&|ϕ(t)|2dt=1. 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 2j/2to preserve the energy conservation concept; namely, |ϕj,k(t)|2dt=1. Since vjvj+1, there exist a detail space wjthat are complementary of vjin vj+1, i.e., vjwj=vj+1. The subspace wjitself is spanned by a dilated and shifted wavelet function family, i.e. {ψj,k(t)}, where ψj,k(t)=2j/2ψ(2jtk); the function ψ(t)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.:ψj,k(t)dt=0.This is because the wavelet function measures local fluctuations; the variations which are assumed to have zero medium. Similar to scaling function and for the same reason, energy of the wavelet functions are unit, i.e., |ψj,k(t)|2dt=1. The approximate and detail subspaces satisfy orthogonally conditions as follows: vjwj& (wjwjforjj). These relations lead to:ψj,k,ϕj,l=0&ψj,k,ψj,l=δklδjj&φj,k,φj,l=δklwhere <f(t).g(t)>=f*(t)g(t)dt(the inner product).Due to the fact that v0v1and w0v1, then any function in v0orw0can be expanded in terms of the basis function of v1, i.e.:ϕ(t)=2khkϕ(2tk)&ψ(t)=2khkϕ(2tk).These important equations are known as: dilation equations, refinement equationsor two-scale relationships[69-72]. The hkand hkare called filter coefficients, and can be obtained, in general, by the following relationships:hk=<ϕ1,k,ϕ>and hk=<ϕ1,k,ψ>.The orthogonality condition of ϕ(xk)and ψ(xk)leads to relationship:hk=(1)khN1kwhere Nis length of the scaling coefficient filters, {hk},k=1,,N.

As mentioned before, the multiresolution decomposition of L2()leads to a set of subspaces with different resolution levels; i.e.,L2=vjwjwj+1.... In this regard, by using one decomposition level, a function f(x)vJmax(a space with sampling step 1/2Jmax) can be expanded as:f(x)=(l=+c(Jmax1,k)ϕJmax1,,l(x))+(n=d(Jmax1,n)ψJmax1,n(x))

By following thestep by step decomposition of approximation space, if the coarsest resolution level is JminJmax1, then the function f(x)vJmaxcan be represented as:

f(x)=(l=+c(Jmin,k)ϕJmin,l(x))+(j=JminJmax1n=d(j,n)ψj,n(x))E5

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:{{d(Jmax1,n)},{d(Jmax2,n)},,{d(Jmin,n)},{c(Jmin,n)}}; 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:c(j,k)=<f(x),ϕj,k(x)>&d(j,k)=<f(x),ψj,k(x)>. 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. Figure 1, illustrates the horizontal acceleration recorded at the El-Centro substation (f1) and corresponding wavelet-based decompositions. There, the symbol a0refers to the approximation space (a0=l=+c(0,k)ϕ0,l(x)) and d0to d9denote the detail spaces (dj=n=+d(j,n)ψj,n(x)); where, the finest and coarsest resolution levels are Jmax=10and Jmin=0, respectively. The superposition of all projected data, f2(i.e., f2=a0+j=09dj) and the difference f2f1are presented as well. It is clear that a0approximates the overall smooth behavior; the projections d0-d9include localfluctuations in different resolutions. There, the frequency content of djis in accordance with the resolution level j. The wavelet used for the decompositions is the Daubechies wavelet of order 12 (will be explained subsequently).

#### 2.1.2. 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.

In orthogonal systems, necessary conditions for the scaling functions are :

1. Normalization condition: ϕ(t)dt=1which leads to k=1Nhk=1; Ndenotes filter length.

2. Orthogonality condition: ϕ(t).ϕ(t+l)dt=δ0,lor equivalently k=1Nhkhk+2l=δ0,l; the parameter δ0,ldenotes the Kronecker delta. For filter set {hk}:k=1,2,,Nwhere Nis an even number, this condition provides N/2independent conditions.

Essential conditions 1 & 2 provide N/2+1independent equations. Other requirements can be assumed to obtain the remaining equations.

One choice is the necessity that the set function {φ(xk)}can exactly reconstruct polynomials of order upto but not greater than p[71, 72]. The polynomial can be represented as:f(x)=α0+α1x++αp1xp1; on the other hands: f(x)=k=ckφ(xk).

By taking the inner product of the wavelet function (ψ(x)) with the above equation, other conditions can be obtained; since:<f(x),ψ(x)>=k=ck<φ(xk),ψ(x)>0, or:<f(x),ψ(x)>=α0ψ(x)dx+α1x.ψ(x)dx++αp1xp1.ψ(x)dx=0.

As αicoefficients are arbitrary, then it is necessary that each of the above integration to be equal to zero:xl.ψ(x)dx=0,l=0,1,,p1; theseequations lead to pequations where p1of them are independent [69, 71, 72]. These equations mean that the first pmoments of the wavelet function must be equal to zero; this condition in the frequency domain leads to relationship [dlψ^(w)/dwl]w=0=0,l=0,1,,p1. It can be shown that these conditions lead to conditions:kN(1)khkkl=0,l=0,1,,p1.

In case that p=N/2, where Nis even (for filter coefficients of lengthN) the resulted wavelet family is known as the Daubechies wavelets. In this case, for Nscaling filter coefficients, Nindependent equations exist, and unique results can be obtained.

In Figure 2 the Daubechies scaling and wavelet functions of order 12 in spatial (φ(x)&ψ(x)) and frequency (|φ^(w)|&|ψ^(w)|) domains are illustrated. It is evident that functions φ(x)&ψ(x), and |φ^(w)|&|ψ^(w)|have localized feature. In this figure φ(x)denotes the first derivative of the scaling function.

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 . To fulfill some of these requirements, the orthogonality requirement can be relaxed and the bi-orthogonal system is used . For numerical purposes, some other requirements can also be considered; for example Dahlke et al.  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 . Figure 2.The Daubechies scaling and wavelet functions of order 12 in spatial and frequency domains, as well as first derivative of the scaling function.

### 2.2. Expressing operators in wavelet spaces

In this subsection, multiresolution analysis of operators will be presented [6-8].

Assume Tdenotes an operator of the following form:T:L2()L2(). 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 levelj(vj) can be represented as:Pj:L2()vjwhere,(Pjf)(x)=f,ϕj,kϕj,k(x). In the same way, the projection of the operator in the detail subspacewj, of resolutionj, is:Qj:L2()wj;Qj=Pj+1Pjwhere,(Qjf)(x)=f,ψj,kψj,k(x). The Qjdefinition is directly resulted from the multiresolution property, i.e., vj1vjand vj=vj1wj1.

For representing the operator Tin the multiresolution form, firstly, a signal xvJmaxis considered, where Jmaxdenotes the finest resolution level, wheredx=1/2Jmax. The data xcan then be projected into the scaling (approximation) and detail spaces of resolution j=Jmax1by one step wavelet transform, i.e.: x=Pj(x)+Qj(x).

Considering a linear operator (function) Tand multiresolution feature, the function T(x)can be presented as follows:T(x)=TJmax=T(Pj+Qj)=T(Pj)+T(Qj)=TPj+TQj.

However T(Pj)&T(Qj)are no longer orthogonal to each other; so each of them can be re-projected to vjand wjas follows:T(Pj)=PjTPj+QjTPj&T(Qj)=PjTQj+QjTQj.

By substituting these relationships inthe equationT(x)=TJmax, we have:

T(x)=(QjTQj+QjTPj+PjTQj)+PjTPjE6

Each term of the above equation belongs to either vjor wjas follows:

Aj=QjfQjwj;Bj=QjfPjwj;Γj=PjfQjvj;
Tj=PjfPjvjE7

In the above equations, Bjand Γjrepresent interrelationship effects of subspacesvjandwj. Using these symbols, the operator Tcan be rewritten as:TJmax=Tj+1=(Aj+Bj+Γj)+TjBy continuously repeating the above mentioned procedure for operators Tj, finally, theTJmaxcan be expressed in the multiresolution representation as follows:

TJmax=i=JminJmax(Ai+Bi+Γi)+TJminE8

where Jmindenotes the coarsest resolution level (i.e., dx=1/2Jmin). This representation is the telescopic form of the operator T.

The schematic shape of the operator Tin 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: Jmin=1, Jmax=4. There, the coefficientsdiandsiare the scale and detail coefficients, respectively; these coefficients are obtained from the common discrete wavelet transform of data x. The d^i&s^iare 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 Tin the wavelet space results to set {Aj,Bj,Γj}jZ, wherejdenotes resolution levels. This form is called NS from, since both of the scale (sj) and detail (dj) coefficients are simultaneously appeared in the formulation, see Figure 3.

The matrix elements of projected operators Aj, Bj, Γj, and Tjare α'j, β'j, γ'jand s'j, respectively; for the derivative operator of order n, dn/dxn,the element definitions are:

αilj=(2j)nψ(2jxi)ψ(n)(2jxl)(2j)dx=(2j)nαilE9
βilj=(2j)nψ(2jxi)ϕ(n)(2jxl)(2j)dx=(2j)nβilE10
γilj=(2j)nϕ(2jxi)ψ(n)(2jxl)(2j)dx=(2j)nγilE11
silj=(2j)nϕ(2jxi)ϕ(n)(2jxl)(2j)dx=(2j)nsilE12

Where:

βl=+ψ(xl)dndxnϕ(x)dx
αl=+ψ(xl)dndxnψ(x)dx,E13
γl=+ϕ(xl)dndxnψ(x)dx,
sl=+ϕ(xl)dndxnϕ(x)dxE14

The coefficientsαilj, βilj, γiljand siljare not independent; the coefficientsαilj, βilj, γiljcan be expressed in terms ofsilj. 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 Aj, Bj, Γjelements.The NS form of the operator d/dxobtained by the Daubechies wavelet of order 12 (Db12) is presented in Figure 4; there Jmin=7&Jmax=10. It is clear that the projected operator is banded in the wavelet space.

To convert {d^j,s^j}JminjJmax1to the standard form {{dj}JminjJmax1,sJmin}, the vector s^jis expanded for JminjJmax1by the following algorithm :

1. d¯Jmax1=s¯Jmax1=0
2. j=Jmax1,Jmax2,,JminE15

(2.1.) If jJmax1then evaluated¯j1&s¯j1from equations¯j+s^j=d¯j1+s¯j1, where d¯j1=Qj1(s¯j+s^j)and s¯j1=Pj1(s¯j+s^j).

(2.2.) evaluatedj1=d¯j1+d^j1,

1. At level j=Jmin, we have sJmin=s¯Jmin+s^Jmin.

The aforementioned telescopic representation is for 1D data. For higher dimensions, the extension is straightforward: the method can independently be implemented for each dimension.

### 2.3. 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.

The semi-group time integration scheme can be used for solving nonlinear equations of form:u,t=Lu+Nf(u)inΩd

where: LandNrepresent the linear and non-linear terms, respectively; u=u(x,t); xd, d=1,2,3;t[0,T]. The initial condition is:u(x,0)=u0(x)inΩ, and the linear boundary condition is:Bu(x,0)=0onΩd1,t[0,T].

Regarding the standard semi-group method, the solution of the above mentioned equations is a non-linear integral equation of the form: u(x,t)=et.L.u(x,0)+0te(tτ)LN(u(x,τ))dτ.

For numerical simulations, the u(x,t)should be discretized in time; the discretized value at time tn=t0+nΔt(Δtis the time step) will be denoted by unu(x,tn). In the same way the discrete form ofN(u(x,t))att=tnis NnN(u(x,tn)).

If the linear operator is a constant, i.e., L=q, the discretized form of the above equation is :un+1=eq.l.Δtun+1l+Δt(γ.Nn+1+m=0M1βm.Nnm), where M+1is the number of time levels considered in the discretization and lM; the coefficients γandβmare functions ofq.Δt. It is clear that the explicit solution is obtained when γ=0; for other choices the scheme is implicit.For case l=1&γ=0(the explicit method) the coefficients γandβmare presented in Table (1).In this table, for linear operator Lthe coefficient Qkis :

Qk=Qk(L.Δt);Qj(L.Δt)=eL.ΔtEj(L.Δt)(L.Δt)j;Ej(L.Δt)=k=0j1(L.Δt)kk!;j=0,1,E16
For j=0,1,2the above mentioned relations yield:Q0(LΔt)=eLΔt;Q1(LΔt)=(eLΔtI)(LΔt)1;Q2(LΔt)=(eLΔtILΔt)(LΔt)2; where Iis the identity matrix.
 order β2 β1 β0 M 1 0 0 Q1 1 2 0 −Q2 Q1+Q2 2 3 Q2/2+Q3 −2(Q2+Q3) Q1+3Q2/2+Q3 3

### Table 1.

Coefficient values for case γ=0&l=1(the explicit scheme), where Qk=Qk(qΔt).

### 2.4. Simulation of 2D SH propagating fronts

The governing equation of the SH scalar wave (anti-plane shear wave) is:

ρ2uyt2=x(μuyx)+z(μuyz)+fyE17

Where uy=uy(x,z)is the out-of plane displacement; μand ρare shear modules and density, respectively. By defining a linear operator Ly, the above mentioned equation can be rewritten as:2uyt2=Ly.uy+fyρwhereLy=1ρx(μx)+1ρz(μz).

For using the semi-group temporal integration scheme, a new variable vy=uy/tis introduced and consequently the above equation will be represented as a system of vectors:

uyt=vy &  vyt=Ly.uy+fyρE18

This system can be rewritten in vector notation as follows :

U=LU+F whereU=(uyvy);L=(0ILy0);F=1ρ(0fy)E19

The simplest explicit semi-group time integration scheme is obtained for caseγ=0&M=1; in this case the discretized form of the wave equation is:Un+1=eΔtLUn+Δt.β0.Fn.

For utilizing the semi-group method the non-linear term eΔtLis approximated by corresponding Taylor expansion :eΔtL=I+Δt.L+Δt22!L2+Δt33!L3+Δt44!L4+....

The coefficient β0can be evaluated as:β0=Q1(LΔt)=(eLΔtI)(LΔt)1. Similarly, the β0can be approximated by its Taylor expansion, i.e.:β0=I+Δt2.L+Δt26L2+Δt324L3+....

#### 2.4.1. 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  is used to simulate infinite boundaries, where the absorbing boundary condition is considered explicitly. Therefore, the wave equation is modified by a damping term Q(x,z).u˙y(x,z,t)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:

2uyt2+Quyt=1ρ(x(μuyx)+z(μuyz))+fyρE20

And the modified vector form of the equation is:U=(uyvy);L=(0ILyQ);F=1ρ(0fy).

#### 2.4.2. Free boundaries

There are different approaches for imposing the free boundary conditions in finite-difference methods; some of which are: 1) using equivalent surface forces (explicit implementation) . 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 ; 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.

#### 2.4.3. 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.

Example:Here scattering of plane SH waves due to a circular tunnel in an infinite media will be presented. The absorbing boundary is used for simulation of infinite domain; the considered function of Q(x,z)is:Q(X,Z)=ax.(ebx.X2+ebx.(Xnx)2)+az.(ebz.Z2+ebz.(Znz)2)where:X=x/dx; Z=z/dz;dx=dz=1/128; ax=az=10000; bx=bz=0.02; nx=nz=128; x[0,1]; z[0,1]. In simulations it assumed: Jmax=7(the finest resolution level); Jmin=4. The Daubechies wavelet of order 12 is considered in calculations. The assumed mechanical properties are: μ=1.8×104kPa&ρ=2ton/m3.

The plane wave condition is simulated by an initial imposed out-of plane harmonic deformation where corresponding wave number is k=64/12. For time integration, the simplest form of the semi-group temporal integration method is used. The snapshots of results (displacement uy(x,z,t)) at different time steps are presented in Figure 5; there the light gray circle represents the tunnel. The displacement uy(x,z,t=0.0048)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×143and 200×200uniform points. Temporal integrations are done by the 4th Runge-Kutta method. Corresponding displacements at t=0.0048are illustrated in Figure 7.Considering Figures 6 & 7, it is clear that the dispersion phenomenon occurs in the common finite difference scheme. There, in each illustration, total CPU computational time presented in the below of each figure. Figure 5.The snapshots of displacement (uy(x,z,t)) at different times. Figure 6.The displacementuy(x,z,t)at timet=0.0048, obtained by the wavelet-based method. Figure 7.The displacementuy(x,z,t)att=0.0048, obtained by the finite difference method; the right and the left figures correspond to grids of size143×143and200×200, respectively.

## 3. Wavelet based simulation of second order hyperbolic systems (wave equations)

In this section, wavelet-based grid adaptation method is survived formodeling 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.

### 3.1. 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 jshould be removed if corresponding detail coefficients are smaller than predefined threshold (ε); wavelet coefficients and grid points have one-to-one correspondence .

In this work, Dubuc-Deslauriers (D-D) interpolating wavelet  is used to grid adaptation. The D-D wavelet of order 2M1(with support Supp(ϕ)=[2M+1,2M1]), is obtained by auto-correlations of Daubechies scaling function of order M(withMvanishing moments).

The D-D scaling function satisfies the interpolating property and has a compact support . In the case of the D-D wavelets, the grid points correspond to the approximation and detail spaces at resolution jare denoted by Vjand Wj, respectively. These sets are locations of the wavelet transform coefficients: the c(j,k)and d(j,k)locations belong to Vjand Wj, respectively. These locations are:

Vj={xj,k[0,1]:xj,k=k/2j};j,k{0,1,,2j}Wj={xj+1,2k+1(0,1):xj+1,2k+1=(2k+1)/2j+1};j,k{0,1,,2j1}E21

Regarding interpolation property of D-D scaling functions, the approximation coefficients (c(Jmin,k)) at points xJmin,kVJminare equal to sampled values of a considered function f(x)at these points, i.e., c(Jmin,k)=f(xJmin,k). The detail coefficients measured at points xj+1,2k+1Wj(of resolution j) is the difference of the function at points xj+1,2k+1(i.e., f(xj+1,2k+1))and corresponding predicted values (the estimated ones in the approximation space). The predicted values are those obtained from the approximation space of resolution j(the corresponding points belong to Vj); the estimated values are denoted by Pfj+1(xj+1,2k+1). In the D-D wavelets, a simple and physical concept exist for such estimation; the estimation at xj+1,2k+1is attained by the local Lagrange interpolation by the known surrounding grid points {xj+1,2k=xj,k}Vj(namely, the even-numbered grid points in Vj+1). For the D-D wavelet of order2M1, 2Mmost neighbor points, including in Vj, are selected in the vicinity of xj+1,2k+1for interpolation; for points far enough from boundary points, the selected points are: {xj+1,2k2n}n{M+1,M+2,,M}. Using such set, the estimation at the point xj+1,2k+1is denoted byPfj+1(xj+1,2n+1), and the detail coefficients are: dj,n=f(xj+1,2n+1)Pfj+1(xj+1,2n+1).

The above mentioned 1D reduction technique can easily be extended to 2D grid points [30, 34]. The boundary wavelets, introduced by Dohono , are also used around edges of finite grid points.

### 3.2. Wavelet-based adaptive-grid method for solving PDEs

At the time step (t=tn), if the solution of PDE isf(x,t), then the procedure for wavelet-based adaptive solution is:

1. Determining the grids, adapted by adaptive wavelet transform, using f(x,tn1)(stepn1). The values of points withoutf(x,tn1), are obtained by locally interpolation (for example, by the cubic spline method);

2. Computing the spatial derivatives in the adapted grid using local Lagrange interpolation scheme, improved by anti-symmetric end padding method . In this regard, extra non-physical fluctuations, deduced by one sided derivatives, are reduced. Here, five points are locally chosen to calculate derivatives and therefore a high-order numerical scheme is achieved [4, 33];

3. Discretizing PDEs in spatial domain first, and then solving semi-discrete systems. The standard time-stepping methods such as Runge-Kutta schemes can be used to solve ODEs at the time t=tn;

4. Denoising the spurious oscillations directly performed in non-uniform grid by smoothing splines (the post processing stage);

5. Repeating the steps from the beginning.

For 1D data of lengthn, smoothing spline of degree2m1, needsm2.noperations , and a wavelet transform (employing pyramidal algorithm) uses noperations. Therefore both procedures are fast and effective. However for cost effective simulation, the grid is adapted after several time steps (e.g. 10-20 steps) based on the velocity of moving fronts. In this case, the moving fronts can be properly captured by adding some extra points to the fronts of adapted grid at each resolution level (e.g., 1 or 2 points to each end at each level).

### 3.3. Smoothing splines

The noisy data are recommended not to be fitted exactly, causing significant distortion particularly in the estimation of derivatives. The smoothing fit is used to remove noisy components in a signal; therefore,interpolation constraint is relaxed. The discrete values of nobservations yj=y(xj)where j=1,2,...,nand x1<x2<...<xnare assumed in order to determine a functionf(x), that yj=f(xj)+εj. εjare random, uncorrelated errors with zero mean and variance σj2. Here, f(x)is the smoothest possible function in fitting the observations to a specific tolerance. It is well known that the solution to this problem is minimizer,f(x),of the functional:

j=1nWj|yjf(xj)|2+(1p)p|(dmf(x)/dxm)|2dx,0p1E22

where, λ=(1p)/p(0λ+) is a Lagrangian parameter, n is the number of observations,Wjis weight factor at pointxjand mis the derivative order.

It can be shown that spline of degreek=2m1, having 2m2continuous derivatives, is an optimal solution; where, n2m. In this chapter the cubic smoothing spline is chosen to have a minimum curvature property; hence, m=2(2m1=3) and fC2[x1,xn][80-82].

According to this formula, the natural cubic spline interpolation is obtained by p=1and the least-squares straight line fit byp=0. In p<1the 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 pshould be selected properly, otherwise it leads to over smoothed or under smoothed results. The former are seen in the scheme presented in Reinsch , according to Hutchinson-Hoog,  and the latter in the scheme offered in Craven-Wahba , according to Lee .

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{Wj}=1, the approximated proper values of pare 0.75- 0.95. The lower values of pare applicable for non-uniformly weighed data, i.e. Wj1. The values of {Wi}and pcan be constant or variable in{(xi,yi)}sequence . Here, the constant weights and smoothing parameter are studied. The {Wj}is assumed as 1in all considered cases.

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 . He showed that the estimation of function and its corresponding derivatives are converged as the interpolating properties and the sampled points are increased .

The smoothing splines work satisfactory for irregular data; this is because the method is a kind of Tikhnov regularization scheme [82, 87, 88].

### 3.4. 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 1/28spatial step in the finest resolution) in three levels; 3- repeating re-adaptation and smoothing processes every ten time steps.

Example:In this example, the wave-fields are presented in inclined two-layered media with sharp transition of physical properties in solid-solid configuration. The numerical methods which do not increase the number of grid points around the interface, have difficulties with the problems of layered media. In such problems, the speeds of elastic waves are largely different. The incident waves, either P or S, can be reflected and refracted from interface in the form of P and S waves.

Schematic shape of considered computational domain is illustrated in Figure8. It is assumed that the top layer is a soft one, while the other one is a stiff layer.It is considered that at point S, the top layer is subjected to an initial imposed deformation ux(x,z,t=0)which is: uz(x,z,t=0)=exp(500((x0.35)2+(z0.25)2)).In the numerical simulation it is assumed that: p=0.85andε=105. As mentioned before, the absorbing boundary condition is considered explicitly for simulation of infinite boundaries.This modification, performed for P-SV wave equations, is:

((λ+2μ)ux,xx+μux,zz)+((λ+μ)uz,xz)=ρ(u¨x+Q(x,z).u˙x)((λ+2μ)uz,zz+μuz,xx)+((λ+μ)ux,xz)=ρ(u¨z+Q(x,z).u˙z)E23

In the above equation, it is assumed that: Q=ax(ebx.x2+ebx.(1x)2)+az(ebz.(1z)2), where ax=az=30, andbx=110,bz=70. The free boundary is imposed by equivalent force in the free surface boundary .

The snapshots of solutions uxand uzand corresponding adapted grids are shown in Figures 9-11, respectively. In each figure, the illustrations (a) to (d) correspond to times0.298, 0.502, 0.658, and 0.886sec, respectively.It is obvious that, the points are properly adapted and most of the energy is confined in the top layer, the soft one.

## 4. 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. Figure 8.Schematic shape of a inclined two-layered media, solid-solid configuration. The soft layer is above a stiff layer. Figure 9.Snapshots of solutionuxat times: a)0.298, b)0.502, c) 0.658, d) 0.886. Figure 10.Snapshots of solutionuzat times: a)0.298, b) 0.502, c) 0.658, d) 0.886. Figure 11.Adapted grid points at times: a)0.298, b) 0.502, c) 0.658, d) 0.886.

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2013 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Hassan Yousefi and Asadollah Noorzad (February 13th 2013). Wavelet Based Simulation of Elastic Wave Propagation, Wave Propagation Theories and Applications, Yi Zheng, IntechOpen, DOI: 10.5772/52097. Available from:

### chapter statistics

1Crossref citations

### Related Content

Next chapter

#### Shear Wave Propagation in Soft Tissue and Ultrasound Vibrometry

By Yi Zheng, Xin Chen, Aiping Yao, Haoming Lin, Yuanyuan Shen, Ying Zhu, Minhua Lu, Tianfu Wang and Siping Chen

#### Acoustic Waves

Edited by Don Dissanayake

First chapter

#### The Eigen Theory of Waves in Piezoelectric Solids

By Shaohua Guo

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.