InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Mechanical Engineering » "Wave Propagation Theories and Applications", book edited by Yi Zheng, ISBN 978-953-51-0979-2, Published: February 13, 2013 under CC BY 3.0 license. © The Author(s).

Chapter 15

Wavelet Based Simulation of Elastic Wave Propagation

By Hassan Yousefi and Asadollah Noorzad
DOI: 10.5772/52097

Article top


The El Centro acceleration and corresponding multiresolution representation.
Figure 1. The El Centro acceleration and corresponding multiresolution representation.
The Daubechies scaling and wavelet functions of order 12 in spatial and frequency domains, as well as first derivative of the scaling function.
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.
Schematic shape of a NS form of the operatorT.
Figure 3. Schematic shape of a NS form of the operatorT.
The NS form of operator d/dx obtained by Db12; it is assumed: Jmin=7&Jmax=10.
Figure 4. The NS form of operator d/dx obtained by Db12; it is assumed: Jmin=7&Jmax=10.
The snapshots of displacement (uy(x,z,t)) at different times.
Figure 5. The snapshots of displacement (uy(x,z,t)) at different times.
The displacement uy(x,z,t) at time t=0.0048, obtained by the wavelet-based method.
Figure 6. The displacement uy(x,z,t) at time t=0.0048, obtained by the wavelet-based method.
The displacement uy(x,z,t) at t=0.0048, obtained by the finite difference method; the right and the left figures correspond to grids of size 143×143 and 200×200, respectively.
Figure 7. The displacement uy(x,z,t) at t=0.0048, obtained by the finite difference method; the right and the left figures correspond to grids of size 143×143 and 200×200, respectively.
Schematic shape of a inclined two-layered media, solid-solid configuration. The soft layer is above a stiff layer.
Figure 8. Schematic shape of a inclined two-layered media, solid-solid configuration. The soft layer is above a stiff layer.
Snapshots of solutionuxat times: a)0.298, b)0.502, c) 0.658, d) 0.886.
Figure 9. Snapshots of solutionuxat times: a)0.298, b)0.502, c) 0.658, d) 0.886.
Snapshots of solution uz at times: a)0.298, b) 0.502, c) 0.658, d) 0.886.
Figure 10. Snapshots of solution uz at times: a)0.298, b) 0.502, c) 0.658, d) 0.886.
Adapted grid points at 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.

Wavelet Based Simulation of Elastic Wave Propagation

Hassan Yousefi 1 and Asadollah Noorzad1

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 [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 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 [3]. 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 [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.

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 [14], 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 [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-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-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-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.

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 vj and the details dk(t) are in wk (detail sub-spaces of level k). The corresponding scale of resolution level j is 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 vj denotes 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/2 to preserve the energy conservation concept; namely, |ϕj,k(t)|2dt=1. Since vjvj+1, there exist a detail space wj that are complementary of vj in vj+1, i.e., vjwj=vj+1. The subspace wj itself 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 v0 orw0 can 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 equations or two-scale relationships [69-72]. The hk and hk are 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 N is 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:


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 a0 refers to the approximation space (a0=l=+c(0,k)ϕ0,l(x)) and d0to d9 denote the detail spaces (dj=n=+d(j,n)ψj,n(x)); where, the finest and coarsest resolution levels are Jmax=10 and Jmin=0, respectively. The superposition of all projected data, f2(i.e., f2=a0+j=09dj) and the difference f2f1 are presented as well. It is clear that a0 approximates 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 [71]:

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

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


Figure 1.

The El Centro acceleration and corresponding multiresolution representation.

Essential conditions 1 & 2 provide N/2+1 independent 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 p equations where p1 of them are independent [69, 71, 72]. These equations mean that the first p moments 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 N is even (for filter coefficients of lengthN) 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 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 [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].


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 Qj definition is directly resulted from the multiresolution property, i.e., vj1vj and vj=vj1wj1.

For representing the operator T in the multiresolution form, firstly, a signal xvJmax is considered, where Jmax denotes the finest resolution level, wheredx=1/2Jmax. The data x can then be projected into the scaling (approximation) and detail spaces of resolution j=Jmax1 by one step wavelet transform, i.e.: x=Pj(x)+Qj(x).

Considering a linear operator (function) T and 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 vj and wj as follows:T(Pj)=PjTPj+QjTPj&T(Qj)=PjTQj+QjTQj.

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


Each term of the above equation belongs to either vj or wj as follows:


In the above equations, Bj and Γj represent 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, theTJmax can be expressed in the multiresolution representation as follows:


where Jmin denotes 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 coefficientsdiandsi are 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.


Figure 3.

Schematic shape of a NS form of the operatorT.

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




The coefficientsαilj, βilj, γiljand silj are 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, Γj elements.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^j is expanded for JminjJmax1 by the following algorithm [73]:

  1. Setd¯Jmax1=s¯Jmax1=0 (the initialization step),
  2. For

(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


Figure 4.

The NS form of operator d/dx obtained by Db12; it is assumed: Jmin=7&Jmax=10.

where: L andN represent 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(Δt is the time step) will be denoted by unu(x,tn). In the same way the discrete form ofN(u(x,t)) att=tn is NnN(u(x,tn)).

If the linear operator is a constant, i.e., L=q, the discretized form of the above equation is [74]:un+1=eq.l.Δtun+1l+Δt(γ.Nn+1+m=0M1βm.Nnm), where M+1 is 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 L the coefficient Qk is [74]:

For j=0,1,2 the 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 I is the identity matrix.

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:


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ρwhereLy=1ρx(μx)+1ρz(μz).

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

uyt=vy &ρ

This system can be rewritten in vector notation as follows [48]:

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

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ΔtL is approximated by corresponding Taylor expansion [48]:eΔtL=I+Δt.L+Δt22!L2+Δt33!L3+Δt44!L4+....

The coefficient β0 can be evaluated as:β0=Q1(LΔt)=(eLΔtI)(LΔt)1. Similarly, the β0 can 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 [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 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:


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) [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.

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×143 and 200×200 uniform points. Temporal integrations are done by the 4th Runge-Kutta method. Corresponding displacements at t=0.0048 are 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 displacement uy(x,z,t) at time t=0.0048, obtained by the wavelet-based method.


Figure 7.

The displacement uy(x,z,t) at t=0.0048, obtained by the finite difference method; the right and the left figures correspond to grids of size 143×143 and 200×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 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 wavelet of order 2M1 (with support Supp(ϕ)=[2M+1,2M1]), is obtained by auto-correlations of Daubechies scaling function of order M (withM vanishing moments).

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 Vj and Wj, respectively. These sets are locations of the wavelet transform coefficients: the c(j,k)and d(j,k)locations belong to Vj and Wj, respectively. These locations are:


Regarding interpolation property of D-D scaling functions, the approximation coefficients (c(Jmin,k)) at points xJmin,kVJmin are 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+1 is 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+1 for 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+1 is 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 [78], 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 [36]. 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 [79], and a wavelet transform (employing pyramidal algorithm) uses n operations. 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 n observations yj=y(xj) where j=1,2,...,n and x1<x2<...<xn are assumed in order to determine a functionf(x), that yj=f(xj)+εj. εj are 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:


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

It can be shown that spline of degreek=2m1, having 2m2 continuous 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=1 and the least-squares straight line fit byp=0. In p<1 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{Wj}=1, the approximated proper values of p are 0.75- 0.95. The lower values of p are applicable for non-uniformly weighed data, i.e. Wj1. The values of {Wi} and p can be constant or variable in{(xi,yi)}sequence [85]. Here, the constant weights and smoothing parameter are studied. The {Wj}is assumed as 1 in 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 [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].

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.85 andε=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:


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 [48].

The snapshots of solutions ux and uz and 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 solution uz at 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.


1 - W. Sweldens, P. Schröder, Wavelet Based Simulation of Elastic Wave Propagationin: Wavelets in Computer Graphics. ACM SIGGRAPH Course Notes, ACM, 19961587
2 - K. Amaratunga, J. R. Williams, S. Qian, J. Weiss, Wavelet, Wavelet Based Simulation of Elastic Wave PropagationIESL Technical Report 92-059205Intelligent Engineering Systems Laboratory, MIT 1992
3 - S. Goedecker, Wavelets, Application. Their, the. for, of. Solution, Differential. Partial, in. Equations, Max. Physics, Max-Planck Institute for Solid State Research, Stuttgart, Germany;
4 - L. M. Jameson, T. Miyama, Wavelet Based Simulation of Elastic Wave PropagationMonthly Weather Review2000128515361548
5 - W. Cai, J. Wang, Wavelet Based Simulation of Elastic Wave PropagationSIAM Journal on Numerical Analysis1996333937970
6 - G. Beylkin, R. Coifman, V. Rokhin, Wavelet Based Simulation of Elastic Wave PropagationCommunications on Pure and Applied Mathematics199164141
7 - G. Beylkin, Wavelet Based Simulation of Elastic Wave PropagationSIAM Journal on Numerical Analysis19926617161740
8 - G. Beylkin, J. M. Keiser, An Adaptive Pseudo-Wavelet Approach for Solving Nonlinear Partial Differential Equations. In: Dahmen W, Kurdila A, Oswald P (eds.), Multiscale Wavelet Methods for Partial Differential Equations (Wavelet Analysis and Its Applications, 6San Diego: Academic Press; 1997137197
9 - S. G. A. Mallat, for. Theory, Signal. Multiresolution, The. Decomposition, Representation. I. E. E. Wavelet, IEEE Transactions on Pattern Analysis and Machine Intelligence 198927647693
10 - Pan GW.Wavelet Based Simulation of Elastic Wave PropagationNew Jersey: John Wiley & Sons; 2003
11 - Xu JC, Shann WC.Wavelet-Galerkin Methods for Two-point Boundary Value Problems. Numerische Mathematik1994372703
12 - L. U. Dianfeng, T. Ohyoshi, L. Zhu, Treatment of Boundary Condition in the Application of Wavelet-Galerkin Method to a SH Wave Problem. 1996Akita Univ. (Japan).
13 - V. Mishra, Wavelet. Sabina, Solutions. Galerkin, Ordinary. of, Equations. Differential, International Journal of Mathematical Analysis201159407424
14 - K. Amaratunga, J. R. Williams, Wavelet, Wavelet Based Simulation of Elastic Wave PropagationArchives of Computational Methods in Engineering199743243285
15 - J. R. Williams, K. Amaratunga, Based. Wavelet, Function. Green’s, to. . D. P. D. Approach, Es, Engineering Computations199310349
16 - S. Dahlke, I. Weinreich, Methods. Wavelet-Galerkin, Adapted. an, Wavelet. Biorthogonal, Basis, Applied and Computational Harmonic Analysis 199413237267
17 - J. R. Williams, K. A. Amaratunga, Wavelet. Multiscale, with. Solver, Complexity. O(n, Journal of computational physics 199512230
18 - S. Qian, J. Weiss, Wavelets and the Numerical Solution of Boundary Value Problems. Applied Mathematics Letter 1993614752
19 - M. Holmstrom, J. Walden, Wavelet Based Simulation of Elastic Wave PropagationJournal of Scientific Computing 19981311949
20 - M. Mehra, B. V. R. Kumar, Time, Wavelet Based Simulation of Elastic Wave PropagationCommunications in Numerical Methods in Engineering200521313
21 - B. V. R. Kumar, M. Mehra, Galerkin. Wavelet-Taylor, for. Method, Burgers. the, B. I. Equation, BITNumerical Mathematics 200545543
22 - M. Mehra, B. V. R. Kumar, Wavelet Based Simulation of Elastic Wave PropagationApplied Mathematics and Computation20071891292
23 - O. V. Vasilyev, N. K. R. Kevlahan, Wavelet Based Simulation of Elastic Wave PropagationJournal of Computational Physics20052062412431
24 - Alam JM, Kevlahan NKR, Vasilyev OV.Wavelet Based Simulation of Elastic Wave PropagationJournal of Computational Physics20062142829857
25 - S. Bertoluzza, L. Castro, Adaptive Wavelet Collocation for Elasticity: First Results. Technical Report 1276, Pub. IAN-CNR de Pavia. 2002
26 - M. Griebel, F. Koster, Adaptive Wavelet Solvers for the Unsteady Incompressible Navier-Stokes Equations. In: Malek J, Nečas J, Rokyta M (eds.) Advances in Mathematical Fluid Mechanics. Berlin: Springer; 200067118
27 - A. Latto, H. L. Resnikoff, E. Tenenbaum, The evaluation of connection coefficients of compactly supported wavelets. In: Proceedings of the French-USA Workshop on Wavelets and Turbulence, 1991, Princeton, New York: Springer-Verlag; 1992
28 - Jang GW, Kim JE, Kim YY.Wavelet Based Simulation of Elastic Wave PropagationInternational Journal for Numerical Methods in Engineering200459225
29 - Kim JE, Jang GW, Kim YY.Wavelet Based Simulation of Elastic Wave PropagationInternational Journal of Solids and Structures2003406473
30 - J. C. Santos, P. Cruz, MA Oliveira. P. J. Alves, F. D. Magalhães, A. Mendes, Wavelet Based Simulation of Elastic Wave PropagationComputer Methods in Applied Mechanics and Engineering20041933405425
31 - P. Cruz, A. Mendes, F. D. Magalhães, Wavelet, Wavelet Based Simulation of Elastic Wave PropagationAIChE Journal2002484774785
32 - P. Cruz, A. Mendes, F. D. Magalhães, Wavelet Based Simulation of Elastic Wave PropagationChemical Engineering Science2001561033053309
33 - L. M. A. Jameson, Very. Wavelet-Optimized, Order. High, Grid. Adaptive, Numerical. Order, S. I. A. Method, SIAM Journal on Scientific Computing 199819619802013
34 - M. Holmstrom, Wavelet Based Simulation of Elastic Wave PropagationSIAM Journal on Scientific Computing1999212405420
35 - S. Operto, J. Virieux, B. Hustedt, F. Malfanti, Wavelet Based Simulation of Elastic Wave PropagationGeophysical Journal International2002148476
36 - H. Yousefi, A. Noorzad, J. Farjoodi, . D. Simulating, Propagation. Waves, Elastic. in, Media. Solid, Wavelet. Using, Adaptive. Based, Method, Journal of Scientific Computing201042404
37 - H. Yousefi, A. Noorzad, J. Farjoodi, M. Vahidi, Multiresolution, Multiresolution-Based Adaptive Simulation of Wave Equation. Applied Mathematics & Information Sciences 2012S) 4758
38 - Zh. L. Pei, L. Y. Fu, G. X. Yu, L. X. A. Zhang, Adaptive. Wavelet-Optimized, Method. Grid, Finite. for, Simulation. Difference, Wave. of, Propagation, Bulletin of the Seismological Society of America2009991302313
39 - MA Cruz. P. Alves, A. Mendes, F. D. Magalhães, F. T. Pinho, P. J. Oliveira, Wavelet Based Simulation of Elastic Wave PropagationComputer Methods in Applied Mechanics and Engineering20021913909
40 - P. Cruz, MA Mendes. A. Alves, F. D. Magalhães, A. Mendes, Wavelet Based Simulation of Elastic Wave PropagationChemical Engineering Science2003581777
41 - A. Harten, Wavelet Based Simulation of Elastic Wave PropagationJournal of Computational Physics 19941152319338
42 - A. Cohen, S. M. Kaber, S. Muller, M. Postel, Wavelet Based Simulation of Elastic Wave PropagationMathematics of Computation 200372183
43 - S. Muller, Y. Stiriba, Wavelet Based Simulation of Elastic Wave PropagationJournal of Scientific Computing 2007303493531
44 - G. W. Wei, Y. Gu, Conjugate Filter Approach For Solving Burgers’ Equation. Journal of Computational and Applied Mathematics 2002
45 - Y. Gu, G. W. Wei, Conjugate Filter Approach for Shock Capturing. Communications in Numerical Methods in Engineering 2003
46 - D. C. Diez, M. Gunzburger, A. Kunoth, An Adaptive Wavelet Viscosity Method for Hyperbolic Conservation Laws. Numerical Methods for Partial Differential Equations 200824613881404
47 - T. K. Hong, B. L. N. A. Kennett, Method. Wavelet-Based, Simulation. for, Two. of-Dimensional, Wave. Elastic, Propagation, Geophysical Journal International 2002150610
48 - Hong TK, Kennett BLN.On a Wavelet-Based Method for the Numerical Simulation of Wave Propagation. Journal of Computational Physics 2002183577
49 - Hong TK, Kennett BLN.Scattering Attenuation of 2D Elastic Waves: Theory and Numerical Modeling Using a Wavelet-Based Method. Bulletin of the Seismological Society of America 2003932922938
50 - Hong TK, Kennett BLN.Modelling of Seismic Waves in Heterogeneous Media Using a Wavelet-Based Method: Application to Fault and Subduction Zones. Geophysical Journal International 2003154483
51 - Hong TK, Kennett BLN.Scattering of Elastic Waves in Media with a Random Distribution of Fluid-Filled Cavities: Theory and Numerical Modeling. Geophysical Journal International 2004159961
52 - Y. Wu, G. A. Mc Mechan, Wave Extrapolation in the Spatial Wavelet Domain with Application to Poststack Reverse-Time Migration. Geophysics 1998632589600
53 - S. Gopalakrishnan, M. Mitra, Wavelet Methods for Dynamical Problems with Application to Metallic, Composite, and Nano Composite Structures. New York: CRC Press; 2010
54 - P. Moczo, J. Kristek, L. Halada, The Finite-Difference Method for Seismologists, an Introduction. Comenius University Bratislava; 2004
55 - P. Moczo, J. Kristeka, M. Galisb, P. Pazaka, M. Balazovjech, Finite. The-Difference, Modeling. Finite-Element, Seismic. of, Propagation. Wave, Motion. Earthquake, Acta Physica Slovaca 2007572177406
56 - Prentice JSC.Truncation and Roundoff Errors in Three-Point Approximations of First and Second Derivatives. Applied Mathematics and Computation 20112174576
57 - Yamaleev NK.Minimization of the Truncation Error by Grid Adaptation. Journal of Computational Physics 2001
58 - R. Vichnevetsky, Propagation through Numerical Mesh Refinement for Hyperbolic Equations. Mathematics and Computers in Simulation 1981
59 - R. Vichnevetsky, Propagation and Spurious Reflection in Finite Element Approximations of Hyperbolic Equations. Computers & Mathematics with Applications 1985
60 - R. Vichnevetsky, Wave Propagation Analysis of Difference Schemes for Hyperbolic Equations: A Review. International Journal for Numerical Methods in Fluids 1987a
61 - R. Vichnevetsky, Wave Propagation and Reflection in Irregular Grids for Hyperbolic Equations. Applied Numerical Mathematics 1987b
62 - R. Grotjahn, J. J. Obrien, Some Inaccuracies in Finite Differencing Hyperbolic Equations. Monthly Weather Review 1976
63 - Trefethen LN.Group Velocity of Finite Difference Schemes. SIAM Review 1982
64 - Bazant ZP.Spurious Reflection of Elastic Waves in Nonuniform Finite Element Grids. Computer Methods in Applied Mechanics and Engineering 19781691
65 - Z. P. Bazant, Z. Celep, Spurious Reflection of Elastic Waves in Nonuniform Meshes of Constant and Linear Strain Finite Elements. Computers & Structures 1982
66 - D. Gottlieb, J. S. Hesthaven, Spectral Methods for Hyperbolic Problems. Journal of Computational and Applied Mathematics 2001
67 - B. Engquist, H. O. Kreiss, Difference and Finite Element Methods for Hyperbolic Differential Equations. Computer Methods in Applied Mechanics and Engineering 1979
68 - Chin RCY.Dispersion and Gibbs Phenomenon Associated with Difference Approximations to Initial Boundary-Value Problems for Hyperbolic Equations. Journal of Computational Physics 1975
69 - S. A. Mallet, Tour. Wavelet, Signal. of, Processing, San Diego: Academic Press; 1998
70 - B. Jawerth, W. Sweldens, An Overview of Wavelet Based Multiresolution Analysis. SIAM Review 1994363377412
71 - J. R. Williams, K. Amaratunga, Introduction to Wavelets in Engineering. IESL Technical Report 92-079207Intelligent Engineering Systems Lab oratory, MIT, 1992
72 - I. Daubechies, Orthonormal Bases of Compactly Supported Wavelets. Communications on Pure and Applied Mathematics 198841909
73 - D. Gines, G. Beylkin, J. L. U. Dunn, of. Factorization, Forms. Non-Standard, Multiresolution. Direct, Solvers, Applied and Computational Harmonic Analysis 199852156201
74 - G. Beylkin, J. M. Keiser, L. A. Vozovoi, Class. New, Time. of, Schemes. Discretization, the. for, of. Solution, P. D. Nonlinear, Es, Journal of Computational Physics 1998147362
75 - G. Beylkin, J. M. Keiser, On the Adaptive Numerical Solution of Nonlinear Partial Differential Equations in Wavelet Bases. Journal of Computational Physics 1997132233
76 - J. Sochacki, R. Kubichek, J. George, W. R. Fletcher, S. Smithson, Absorbing Boundary Conditions and Surface Waves. Geophysics 19875216071
77 - Boore DM.Finite Difference Methods for Seismic Wave Propagation in Heterogeneous Materials. In: Bolt BA (ed.) Methods in Computational Physics, 11Seismology: SurfaceWaves and Earth Oscillations. New York: Academic Press; 1972137
78 - Donoho DL. Interpolating Wavelet Transforms. Technical Report 408. Department of Statistics, StanfordUniversity. 1992
79 - Hutchinson MF, de Hoog FR.Smoothing Noisy Data with Spline Functions. Numerische Mathematik 1985
80 - Reinsch CH.Smoothing by Spline Functions. Numerische Mathematik 1967
81 - Reinsch CH.Smoothing by Spline Functions. II. Numerische Mathematik 1971
82 - M. Unser, A. Splines, Fit. Perfect, Signal. for, Processing. I. E. E. Image, IEEE Signal Proc Mag 1999
83 - P. Craven, G. Wahba, Smoothing Noisy Data with Spline Functions: Estimating the Correct Degree of Smoothing by the Method of Generalized Cross Validation. Numerische Mathematik 1979
84 - Lee TCM.Smoothing Parameter Selection for Smoothing Splines: A Simulation Study. Computational Statistics & Data Analysis 2003
85 - Lee TCM.Improved Smoothing Spline Regression by Combining Estimates of Different Smoothness. Statistics & Probability Letters 2004672133140
86 - Ragozin DL.Error Bounds for Derivative Estimates Based on Spline Smoothing of Exact or Noisy Data. Journal of Approximation Theory 1983
87 - Petrov YP, Sizikov VS. Well-Posed, Ill-Posed, and Intermediate Problems with Applications.VSP; 2005
88 - Hansen PC.Rank-Deficient and Discrete Ill-Posed Problems, Philadelphia: SIAM; 1998