Wavelet Based Simulation of Elastic Wave Propagation

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, $$L^2(ℝ)$$, are decomposed as a sequence of detail subspaces, denoted by $$\left\{w_k\right\}$$, and an approximation subspace, indicated with$$v_j$$. The approximation of $$f(t)$$ at resolution level $$j$$, $$f(t)$$, is in $$v_j$$ and the details $$d_k(t)$$ are in $$w_k$$ (detail sub-spaces of level $$k$$). The corresponding scale of resolution level $$j$$ is usually chosen to be of order $$2^{-j}$$ [69, 70].In orthogonal wavelet systems, the multiresolution analysis of $$L^2(ℝ)$$ is nested sequences of the subspaces $$\left\{v_j\right\}$$ such that:

1. $$\mathrm{...}\subset v_{-1}\subset v_0\subset v_1\mathrm{...}\subset L^2(ℝ)$$

2. $$v_{-\infty }=\left\{0\right\},v_{+\infty }=L^2$$

3. $$f(t)\in v_j\iff f(2t)\in v_{j+1}$$

4. $$f(t)\in v_0\Rightarrow f(t-k)\in v_0$$

5. Exists a function$$\varphi (t)$$, called the scaling function such that set $${\left\{\varphi (t-k)\right\}}_{k\in Z}$$is a basis of $$v_0$$.

The sub-space $$v_j$$ denotes the space spanned by family $$\{{\varphi }_{j,k}(t)\}$$, i.e., $$v_j=\overline{\underset{k}{span}\{{\varphi }_{j,k}(t)\}}$$ where $${\varphi }_{j,k}(t)=2^{j/2}\varphi (2^{j}t-k)$$. The $${\varphi }_{j,k}(t)$$ is a scaled and shifted version of the $$\varphi (t)$$; thereby the function $$\varphi (t)$$ is known as the father wavelet. The scale functions ($${\varphi }_{j,k}(t)$$) are localized in both spatial (or time) and frequency (scale) spaces. The function $$\varphi (t)$$ is usually designed so that:$${\displaystyle \int \varphi (t)}dt=1$$&$${{\displaystyle \int \left|\varphi (t)\right|}}^{2}dt=1$$. The second equation implies that the scaling function ($$\varphi (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, $${\varphi }_{j,k}(t)$$ is usually normalized with the coefficient $$2^{j/2}$$ to preserve the energy conservation concept; namely, $${{\displaystyle \int \left|{\varphi }_{j,k}(t)\right|}}^{2}dt=1$$. Since $$v_j\subset v_{j+1}$$, there exist a detail space $$w_j$$ that are complementary of $$v_j$$ in $$v_{j+1}$$, i.e., $$v_j\oplus w_j=v_{j+1}$$. The subspace $$w_j$$ itself is spanned by a dilated and shifted wavelet function family, i.e. $$\left\{{\psi }_{j,k}(t)\right\}$$, where $${\psi }_{j,k}(t)=2^{j/2}\psi (2^{j}t-k)$$; the function $$\psi (t)$$ is usually referred as the mother wavelet. The wavelet function, $${\psi }_{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.:$${\displaystyle \int {\psi }_{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., $${{\displaystyle \int \left|{\psi }_{j,k}(t)\right|}}^{2}dt=1$$. The approximate and detail subspaces satisfy orthogonally conditions as follows: $$v_j\perp w_j$$& ($$w_j\perp w_{j^{\prime }}forj\ne j^{\prime }$$). These relations lead to:$$⟨{\psi }_{j,k},{\varphi }_{j,l}⟩=0$$&$$⟨{\psi }_{j,k},{\psi }_{j^{\prime },l}⟩={\delta }_{kl}{\delta }_{j{j}^{\prime }}$$&$$⟨{\phi }_{j,k},{\phi }_{j,l}⟩={\delta }_{kl}$$where $$<f(t).g(t)>={\displaystyle \int f^{*}(t)g(t)dt}$$ (the inner product).Due to the fact that $$v_0\subset v_1$$and $$w_0\subset v_1$$, then any function in $$v_0$$ or$$w_0$$ can be expanded in terms of the basis function of $$v_1$$, i.e.:$$\varphi (t)=\sqrt{2}{\displaystyle \sum _k{h}_k\varphi (2t-k)}$$&$$\psi (t)=\sqrt{2}{\displaystyle \sum _k{h^{\prime }}_k\varphi (2t-k)}$$.These important equations are known as: dilation equations, refinement equations or two-scale relationships [69-72]. The $$h_k$$ and $${h^{\prime }}_k$$ are called filter coefficients, and can be obtained, in general, by the following relationships:$$h_k=<{\varphi }_{1,k},\varphi >$$ and $${h^{\prime }}_k=<{\varphi }_{1,k},\psi >$$.The orthogonality condition of $$\varphi (x-k)$$ and $$\psi (x-k)$$ leads to relationship:$${h^{\prime }}_k={(-1)}^k{h}_{N-1-k}$$where $$N$$ is length of the scaling coefficient filters, $$\left\{h_k\right\},k=1,\dots ,N$$.

As mentioned before, the multiresolution decomposition of $$L^2(ℝ)$$ leads to a set of subspaces with different resolution levels; i.e.,$$L^2=v_j\oplus w_j\oplus w_{j+1}\oplus \mathrm{...}$$. In this regard, by using one decomposition level, a function $$f(x)\in v_{J_{\max }}$$ (a space with sampling step $$1/2^{J_{\max }}$$) can be expanded as:$$f(x)=\left({\displaystyle \sum _{l=-\infty }^{+\infty }c(J_{\max }-1,k){\varphi }_{J_{\max }-1,,l}(x)}\right)+\left({\displaystyle \sum _{n=-\infty }^{\infty }d(J_{\max }-1,n){\psi }_{J_{\max -1},n}(x)}\right)$$

By following thestep by step decomposition of approximation space, if the coarsest resolution level is $$J_{\min }\le J_{\max -1}$$, then the function $$f(x)\in v_{J_{\max }}$$can be represented as:

$$f(x)=\left({\displaystyle \sum _{l=-\infty }^{+\infty }c(J_{\min },k){\varphi }_{J_{\min },l}(x)}\right)+\left({\displaystyle \sum _{j=J_{\min }}^{J_{\max }-1}{\displaystyle \sum _{n=-\infty }^{\infty }d(j,n){\psi }_{j,n}(x)}}\right)$$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:$$\left\{\left\{d(J_{\max }-1,n)\right\},\left\{d(J_{\max }-2,n)\right\},\dots ,\left\{d(J_{\min },n)\right\},\left\{c(J_{\min },n)\right\}\right\}$$; 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),{\varphi }_{j,k}(x)>$$&$$d(j,k)=<f(x),{\psi }_{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 ($$f_1$$) and corresponding wavelet-based decompositions. There, the symbol $$a_0$$ refers to the approximation space ($$a_0={\displaystyle {\sum }_{l=-\infty }^{+\infty }c(0,k){\varphi }_{0,l}(x)}$$) and $$d_0$$to $$d_9$$ denote the detail spaces ($$d_j={\displaystyle {\sum }_{n=-\infty }^{+\infty }d(j,n){\psi }_{j,n}(x)}$$); where, the finest and coarsest resolution levels are $$J_{\max }=10$$ and $$J_{\min }=0$$, respectively. The superposition of all projected data, $$f_2$$(i.e., $$f_2=a_0+{\displaystyle {\sum }_{j=0}^9{d}_j}$$) and the difference $$f_2-f_1$$ are presented as well. It is clear that $$a_0$$ approximates the overall smooth behavior; the projections $$d_0$$-$$d_9$$include localfluctuations in different resolutions. There, the frequency content of $$d_j$$is 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: $${\displaystyle \int \varphi (t)}dt=1$$ which leads to $${\displaystyle \sum _{k=1}^N{h}_k}=1$$; $$N$$ denotes filter length.

2. Orthogonality condition: $${\displaystyle \int \varphi (t).\varphi (t+l)}dt={\delta }_{0,l}$$ or equivalently $${\displaystyle \sum _{k=1}^N{h}_k}{h}_{k+2l}={\delta }_{0,l}$$; the parameter $${\delta }_{0,l}$$ denotes the Kronecker delta. For filter set $$\left\{h_k\right\}:k=1,2,\dots ,N$$where $$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 $$\left\{\phi (x-k)\right\}$$can exactly reconstruct polynomials of order upto but not greater than $$p$$ [71, 72]. The polynomial can be represented as:$$f(x)={\alpha }_0+{\alpha }_{1}x+\dots +{\alpha }_{p-1}{x}^{p-1}$$; on the other hands: $$f(x)={\displaystyle \sum _{k=-\infty }^{\infty }{c}_k\phi (x-k)}$$.

By taking the inner product of the wavelet function ($$\psi (x)$$) with the above equation, other conditions can be obtained; since:$$<f(x),\psi (x)>={\displaystyle \sum _{k=-\infty }^{\infty }{c}_k<\phi (x-k),\psi (x)>}\equiv 0$$, or:$$<f(x),\psi (x)>={\alpha }_0{\displaystyle \int \psi (x)dx}+{\alpha }_1{\displaystyle \int x.}\psi (x)dx+\dots +{\alpha }_{p-1}{\displaystyle \int x^{p-1}}.\psi (x)dx=0$$.

As $${\alpha }_i$$coefficients are arbitrary, then it is necessary that each of the above integration to be equal to zero:$${\displaystyle \int x^l}.\psi (x)dx=0,l=0,1,\dots ,p-1$$; theseequations lead to $$p$$ equations where $$p-1$$ 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 $${\left[d^l\widehat{\psi }(w)/d{w}^l\right]}_{w=0}=0,l=0,1,\dots ,p-1$$. It can be shown that these conditions lead to conditions:$${\displaystyle \sum _k^N{(-1)}^k{h}_k}{k}^l=0,l=0,1,\dots ,p-1$$.

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

In Figure 2 the Daubechies scaling and wavelet functions of order 12 in spatial ($$\phi (x)$$&$$\psi (x)$$) and frequency ($$\left|\widehat{\phi }(w)\right|$$&$$\left|\widehat{\psi }(w)\right|$$) domains are illustrated. It is evident that functions $$\phi (x)$$&$$\psi (x)$$, and $$\left|\widehat{\phi }(w)\right|$$&$$\left|\widehat{\psi }(w)\right|$$ have localized feature. In this figure $${\phi }^{\prime }(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 $$T$$denotes an operator of the following form:$$T:L^2(ℝ)\to L^2(ℝ)$$. The aim is to represent the operator T in the wavelet spaces; this can be done by projection the operator in the wavelet spaces.

The projection of the operator in the approximation space of resolution level$$j$$ ($$v_j$$) can be represented as:$$P_j:L^2(ℝ)\to v_j$$where,$$(P_{j}f)(x)={\displaystyle \sum _{}⟨f,{\varphi }_{j,k}⟩{\varphi }_{j,k}(x)}$$. In the same way, the projection of the operator in the detail subspace$$w_j$$, of resolution$$j$$, is:$$Q_j:L^2(ℝ)\to w_j;Q_j=P_{j+1}-P_j$$where,$$(Q_{j}f)(x)={\displaystyle \sum ⟨f,{\psi }_{j,k}⟩{\psi }_{j,k}(x)}$$. The $$Q_j$$ definition is directly resulted from the multiresolution property, i.e., $$v_{j-1}\subset v_j$$ and $$v_j=v_{j-1}\oplus w_{j-1}$$.

For representing the operator $$T$$ in the multiresolution form, firstly, a signal $$x\in v_{J_{\max }}$$ is considered, where $$J_{\max }$$ denotes the finest resolution level, where$$dx=1/2^{J_{\max }}$$. The data $$x$$ can then be projected into the scaling (approximation) and detail spaces of resolution $$j=J_{\max }-1$$ by one step wavelet transform, i.e.: $$x=P_j(x)+Q_j(x)$$.

Considering a linear operator (function) $$T$$ and multiresolution feature, the function $$T(x)$$ can be presented as follows:$$T(x)=T_{J_{\max }}=T(P_j+Q_j)=T(P_j)+T(Q_j)=T{P}_j+T{Q}_j$$.

However $$T(P_j)$$&$$T(Q_j)$$ are no longer orthogonal to each other; so each of them can be re-projected to $$v_j$$ and $$w_j$$ as follows:$$T(P_j)=P_{j}T{P}_j+Q_{j}T{P}_j\&T(Q_j)=P_{j}T{Q}_j+Q_{j}T{Q}_j$$.

By substituting these relationships inthe equation$$T(x)=T_{J_{\max }}$$, we have:

$$T(x)=(Q_{j}T{Q}_j+Q_{j}T{P}_j+P_{j}T{Q}_j)+P_{j}T{P}_j$$E6

Each term of the above equation belongs to either $$v_j$$ or $$w_j$$ as follows:

$$A_j=Q_{j}f{Q}_j\in w_j;$$$$B_j=Q_{j}f{P}_j\in w_j;$$$${\Gamma }_j=P_{j}f{Q}_j\in v_j;$$

$$T_j=P_{j}f{P}_j\in v_j$$E7

In the above equations, $$B_j$$ and $${\Gamma }_j$$ represent interrelationship effects of subspaces$$v_j$$and$$w_j$$. Using these symbols, the operator $$T$$can be rewritten as:$$T_{J_{\max }}=T_{j+1}=(A_j+B_j+{\Gamma }_j)+T_j$$By continuously repeating the above mentioned procedure for operators $$T_j$$, finally, the$$T_{J_{\max }}$$ can be expressed in the multiresolution representation as follows:

$$T_{J_{\max }}={\displaystyle \sum _{i=J_{\min }}^{J_{\max }}(A_i+B_i+{\Gamma }_i)}+T_{J_{\min }}$$E8

where $$J_{\min }$$ denotes the coarsest resolution level (i.e., $$dx=1/2^{J_{\min }}$$). This representation is the telescopic form of the operator $$T$$.

The schematic shape of the operator $$T$$in telescopic (multiresolution) form is presented in Figure 3; this form of representation is known as the Non-Standard form (NS form). In this figure, it is assumed that: $$J_{\min }=1$$, $$J_{\max }=4$$. There, the coefficients$$d^i$$and$$s^i$$ are the scale and detail coefficients, respectively; these coefficients are obtained from the common discrete wavelet transform of data $$x$$. The $${\widehat{d}}^i\&{\widehat{s}}^i$$are the NS form of the wavelet coefficients, and should be converted to standard form by a proper algorithm (will be discussed).

The projection of the operator $$T$$in the wavelet space results to set $${\left\{A_j,B_j,{\Gamma }_j\right\}}_{j\in Z}$$, where$$j$$denotes resolution levels. This form is called NS from, since both of the scale ($$s^j$$) and detail ($$d^j$$) coefficients are simultaneously appeared in the formulation, see Figure 3.

Figure 3.

Schematic shape of a NS form of the operator$$T$$.

The matrix elements of projected operators $$A_j$$, $$B_j$$, $${\Gamma }_j$$, and $$T_j$$ are $${\alpha }_{{}^{\text{'}}}^j$$, $${\beta }_{{}^{\text{'}}}^j$$, $${\gamma }_{{}^{\text{'}}}^j$$and $$s_{{}^{\text{'}}}^j$$, respectively; for the derivative operator of order $$n$$, $$d^n/d{x}^n$$,the element definitions are:

$${\alpha }_{il}^j={\left(2^j\right)}^n{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\psi (2^{j}x-i){\psi }^{(n)}(2^{j}x-l)\left(2^j\right)dx={\left(2^j\right)}^n{\alpha }_{i-l}}$$E9

$${\beta }_{il}^j={\left(2^j\right)}^n{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\psi (2^{j}x-i){\varphi }^{(n)}(2^{j}x-l)\left(2^j\right)dx={\left(2^j\right)}^n{\beta }_{i-l}}$$E10

$${\gamma }_{il}^j={\left(2^j\right)}^n{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\varphi (2^{j}x-i){\psi }^{(n)}(2^{j}x-l)\left(2^j\right)dx={\left(2^j\right)}^n{\gamma }_{i-l}}$$E11

$$s_{il}^j={\left(2^j\right)}^n{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\varphi (2^{j}x-i){\varphi }^{(n)}(2^{j}x-l)\left(2^j\right)dx={\left(2^j\right)}^n{s}_{i-l}}$$E12

Where:

$${\beta }_l={\displaystyle {\int }_{-\infty }^{+\infty }\psi (x-l)\frac{d^n}{d{x}^n}\varphi (x)dx}$$

$${\alpha }_l={\displaystyle {\int }_{-\infty }^{+\infty }\psi (x-l)\frac{d^n}{d{x}^n}\psi (x)dx},$$E13

$${\gamma }_l={\displaystyle {\int }_{-\infty }^{+\infty }\varphi (x-l)\frac{d^n}{d{x}^n}\psi (x)dx},$$

$$s_l={\displaystyle {\int }_{-\infty }^{+\infty }\varphi (x-l)\frac{d^n}{d{x}^n}\varphi (x)dx}$$E14

The coefficients$${\alpha }_{il}^j$$, $${\beta }_{il}^j$$, $${\gamma }_{il}^j$$and $$s_{il}^j$$ are not independent; the coefficients$${\alpha }_{il}^j$$, $${\beta }_{il}^j$$, $${\gamma }_{il}^j$$can be expressed in terms of$$s_{il}^j$$. 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 $$A_j$$, $$B_j$$, $${\Gamma }_j$$ elements.The NS form of the operator $$d/dx$$obtained by the Daubechies wavelet of order 12 ($$Db12$$) is presented in Figure 4; there $$J_{\min }=7$$&$$J_{\max }=10$$. It is clear that the projected operator is banded in the wavelet space.

To convert $${\left\{{\widehat{d}}_j,{\widehat{s}}_j\right\}}_{J_{\min }\le j\le J_{\max }-1}$$to the standard form $$\left\{{\left\{d_j\right\}}_{J_{\min }\le j\le J_{\max }-1},s_{J_{\min }}\right\}$$, the vector $${\widehat{s}}_j$$ is expanded for $$J_{\min }\le j\le J_{\max }-1$$ by the following algorithm [73]:

1. $${\overline{d}}_{J_{\max }-1}={\overline{s}}_{J_{\max }-1}=0$$
2. $$j=J_{\max }-1,J_{\max }-2,\dots ,J_{\min }$$E15

(2.1.) If $$j\ne J_{\max }-1$$then evaluate$${\overline{d}}_{j-1}$$&$${\overline{s}}_{j-1}$$from equation$${\overline{s}}_j+{\widehat{s}}_j={\overline{d}}_{j-1}+{\overline{s}}_{j-1}$$, where $${\overline{d}}_{j-1}=Q_{j-1}({\overline{s}}_j+{\widehat{s}}_j)$$and $${\overline{s}}_{j-1}=P_{j-1}({\overline{s}}_j+{\widehat{s}}_j)$$.

(2.2.) evaluate$$d_{j-1}={\overline{d}}_{j-1}+{\widehat{d}}_{j-1}$$,

1. At level $$j=J_{\min }$$, we have $$s_{J_{\min }}={\overline{s}}_{J_{\min }}+{\widehat{s}}_{J_{\min }}$$.

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)\text{in}\Omega \in {ℝ}^d$$

Figure 4.

The NS form of operator $$d/dx$$ obtained by $$Db12$$; it is assumed: $$J_{\min }=7$$&$$J_{\max }=10$$.

where: $$L$$ and$$N$$ represent the linear and non-linear terms, respectively; $$u=u(x,t)$$; $$x\in {ℝ}^d$$, $$d=1,2,3$$;$$t\in [0,T]$$. The initial condition is:$$u(x,0)=u_0(x)\text{in}\Omega$$, and the linear boundary condition is:$$Bu(x,0)=0o\text{n}\partial \Omega \in {ℝ}^{d-1},t\in [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)=e^{t.L}.u(x,0)+{\displaystyle {\int }_0^t{e}^{(t-\tau )L}N(u(x,\tau ))d\tau }$$.

For numerical simulations, the $$u(x,t)$$should be discretized in time; the discretized value at time $$t_n=t_0+n\Delta t$$($$\Delta t$$ is the time step) will be denoted by $$u_n\equiv u(x,t_n)$$. In the same way the discrete form of$$N(u(x,t))$$ at$$t=t_n$$ is $$N_n\equiv N(u(x,t_n))$$.

If the linear operator is a constant, i.e., $$L=q$$, the discretized form of the above equation is [74]:$$u_{n+1}=e^{q.l.\Delta t}{u}_{n+1-l}+\Delta t(\gamma .N_{n+1}+{\displaystyle \sum _{m=0}^{M-1}{\beta }_m.N_{n-m}})$$, where $$M+1$$ is the number of time levels considered in the discretization and $$l\le M$$; the coefficients $$\gamma$$and$${\beta }_m$$are functions of$$q.\Delta t$$. It is clear that the explicit solution is obtained when $$\gamma =0$$; for other choices the scheme is implicit.For case $$l=1$$&$$\gamma =0$$(the explicit method) the coefficients $$\gamma$$and$${\beta }_m$$are presented in Table (1).In this table, for linear operator $$L$$ the coefficient $$Q_k$$ is [74]:

$$Q_k=Q_k(L.\Delta t);Q_j(L.\Delta t)=\frac{e^{L.\Delta t}-E_j(L.\Delta t)}{{(L.\Delta t)}^j};E_j(L.\Delta t)={\displaystyle \sum _{k=0}^{j-1}\frac{{(L.\Delta t)}^k}{k!}};j=0,1,\dots$$E16

For $$j=0,1,2$$ the above mentioned relations yield:$$Q_0(L\Delta t)=e^{L\Delta t};$$$$Q_1(L\Delta t)=(e^{L\Delta t}-I){(L\Delta t)}^{-1};$$$$Q_2(L\Delta t)=(e^{L\Delta t}-I-L\Delta t){(L\Delta t)}^{-2}$$; where $$I$$ is the identity matrix.

$$order$$	$${\beta }_2$$	$${\beta }_1$$	$${\beta }_0$$	$$M$$
1	$$0$$	$$0$$	$$Q_1$$	1
2	$$0$$	$$-Q_2$$	$$Q_1+Q_2$$	2
3	$$Q_2/2+Q_3$$	$$-2(Q_2+Q_3)$$	$$Q_1+3Q_2/2+Q_3$$	3

Table 1.

Coefficient values for case $$\gamma =0\&l=1$$ (the explicit scheme), where $$Q_k=Q_k\left(q\Delta t\right)$$.

2.4. Simulation of 2D SH propagating fronts

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

$$\rho \frac{{\partial }^2{u}_y}{\partial t^2}=\frac{\partial }{\partial x}\left(\mu \frac{\partial u_y}{\partial x}\right)+\frac{\partial }{\partial z}\left(\mu \frac{\partial u_y}{\partial z}\right)+f_y$$E17

Where $$u_y=u_y(x,z)$$ is the out-of plane displacement; $$\mu$$ and $$\rho$$ are shear modules and density, respectively. By defining a linear operator $$L_y$$, the above mentioned equation can be rewritten as:$$\frac{{\partial }^2{u}_y}{\partial t^2}=L_y.u_y+\frac{f_y}{\rho }$$where$$L_y=\frac{1}{\rho }\frac{\partial }{\partial x}\left(\mu \frac{\partial }{\partial x}\right)+\frac{1}{\rho }\frac{\partial }{\partial z}\left(\mu \frac{\partial }{\partial z}\right)$$.

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

$$\frac{\partial u_y}{\partial t}=v_y \& \frac{\partial v_y}{\partial t}=L_y.u_y+\frac{f_y}{\rho }$$E18

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

$$U=LU+F whereU=\left(\begin{array}{c}{u}_y\\ v_y\end{array}\right);L=\left(\begin{array}{cc}0& I\\ L_y& 0\end{array}\right);F=\frac{1}{\rho }\left(\begin{array}{c}0\\ f_y\end{array}\right)$$E19

The simplest explicit semi-group time integration scheme is obtained for case$$\gamma =0\&M=1$$; in this case the discretized form of the wave equation is:$$U_{n+1}=e^{\Delta tL}{U}_n+\Delta t.{\beta }_0.F_n$$.

For utilizing the semi-group method the non-linear term $$e^{\Delta tL}$$ is approximated by corresponding Taylor expansion [48]:$$e^{\Delta tL}=I+\Delta t.L+\frac{\Delta t^2}{2!}{L}^2+\frac{\Delta t^3}{3!}{L}^3+\frac{\Delta t^4}{4!}{L}^4+\mathrm{...}$$.

The coefficient $${\beta }_0$$ can be evaluated as:$${\beta }_0=Q_1(L\Delta t)=(e^{L\Delta t}-I){(L\Delta t)}^{-1}$$. Similarly, the $${\beta }_0$$ can be approximated by its Taylor expansion, i.e.:$${\beta }_0=I+\frac{\Delta t}{2}.L+\frac{\Delta t^2}{6}{L}^2+\frac{\Delta t^3}{24}{L}^3+\mathrm{...}$$.

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).{\dot{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:

$$\frac{{\partial }^2{u}_y}{\partial t^2}+Q\frac{\partial u_y}{\partial t}=\frac{1}{\rho }\left(\frac{\partial }{\partial x}\left(\mu \frac{\partial u_y}{\partial x}\right)+\frac{\partial }{\partial z}\left(\mu \frac{\partial u_y}{\partial z}\right)\right)+\frac{f_y}{\rho }$$E20

And the modified vector form of the equation is:$$U=\left(\begin{array}{c}{u}_y\\ v_y\end{array}\right);L=\left(\begin{array}{cc}0& I\\ L_y& -Q\end{array}\right);F=\frac{1}{\rho }\left(\begin{array}{c}0\\ f_y\end{array}\right)$$.

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)=a_x.\left(e^{b_x.X^2}+e^{b_x.{\left(X-n_x\right)}^2}\right)+a_z.\left(e^{b_z.Z^2}+e^{b_z.{\left(Z-n_z\right)}^2}\right)$$ where:$$X=x/dx$$; $$Z=z/dz$$;$$dx=dz=1/128$$; $$a_x=a_z=10000$$; $$b_x=b_z=-0.02$$; $$n_x=n_z=128$$; $$x\in \left[0,1\right]$$; $$z\in \left[0,1\right]$$. In simulations it assumed: $$J_{\max }=7$$ (the finest resolution level); $$J_{\min }=4$$. The Daubechies wavelet of order 12 is considered in calculations. The assumed mechanical properties are: $$\mu =1.8\times {10}^4\text{}kPa$$&$$\rho =2\text{}\text{}ton/m^3$$.

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 $$u_y(x,z,t)$$) at different time steps are presented in Figure 5; there the light gray circle represents the tunnel. The displacement $$u_y(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\times 143$$ and $$200\times 200$$ uniform points. Temporal integrations are done by the 4^th 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 ($$u_y(x,z,t)$$) at different times.

Figure 6.

The displacement $$u_y(x,z,t)$$ at time $$t=0.0048$$, obtained by the wavelet-based method.

Figure 7.

The displacement $$u_y(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\times 143$$ and $$200\times 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 ($$\epsilon$$); 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 $$2M-1$$ (with support $$Supp(\varphi )=[-2M+1,2M-1]$$), is obtained by auto-correlations of Daubechies scaling function of order $$M$$ (with$$M$$ 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 $V_j$ and $W_j$, respectively. These sets are locations of the wavelet transform coefficients: the $c(j,k)$and $d(j,k)$locations belong to $V_j$ and $W_j$, respectively. These locations are:

$$\begin{array}{l}{V}_j=\{x_{j,k}\in [0,1]:x_{j,k}=k/2^j\};j\in \mathbb{Z},k\in \{0,1,\cdots ,2^j\}\\ W_j=\{x_{j+1,2k+1}\in (0,1):x_{j+1,2k+1}=(2k+1)/2^{j+1}\};j\in \mathbb{Z},k\in \{0,1,\cdots ,2^j-1\}\end{array}$$E21

Regarding interpolation property of D-D scaling functions, the approximation coefficients ($c(J_{\min },k)$) at points $x_{J_{\min },k}\in V_{J_{\min }}$ are equal to sampled values of a considered function $$f(x)$$at these points, i.e., $c(J_{\min },k)=f(x_{J_{\min },k})$. The detail coefficients measured at points $$x_{j+1,2k+1}\in W_j$$ (of resolution $j$) is the difference of the function at points $$x_{j+1,2k+1}$$ (i.e., $f(x_{j+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 $V_j$); the estimated values are denoted by $P{f}_{j+1}(x_{j+1,2k+1})$. In the D-D wavelets, a simple and physical concept exist for such estimation; the estimation at $$x_{j+1,2k+1}$$ is attained by the local Lagrange interpolation by the known surrounding grid points $\left\{x_{j+1,2k}=x_{j,k}\right\}\in V_j$ (namely, the even-numbered grid points in $V_{j+1}$). For the D-D wavelet of order$$2M-1$$, $2M$most neighbor points, including in $V_j$, are selected in the vicinity of $x_{j+1,2k+1}$ for interpolation; for points far enough from boundary points, the selected points are: $\left\{x_{j+1,2k-2n}\right\}n\in \{-M+1,-M+2,\cdots ,M\}$. Using such set, the estimation at the point $x_{j+1,2k+1}$ is denoted by$P{f}_{j+1}(x_{j+1,2n+1})$, and the detail coefficients are: $d_{j,n}=f(x_{j+1,2n+1})-P{f}_{j+1}(x_{j+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=t_n$$), if the solution of PDE is$$f(x,t)$$, then the procedure for wavelet-based adaptive solution is:

1. Determining the grids, adapted by adaptive wavelet transform, using $$f(x,t_{n-1})$$ (step$$n-1$$). The values of points without$$f(x,t_{n-1})$$, 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=t_n;

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 degree$$2m-1$$, needs$$m^2.n$$operations [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 $$y_j=y(x_j)$$ where $$j=1,2,\mathrm{...},n$$ and $$x_1<x_2<\mathrm{...}<x_n$$ are assumed in order to determine a function$$f(x)$$, that $$y_j=f(x_j)+{\epsilon }_j$$. $${\epsilon }_j$$ are random, uncorrelated errors with zero mean and variance $${\sigma }_j^2$$. 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:

$${\displaystyle \sum _{j=1}^n{W}_j{\left|y_j-f(x_j)\right|}^2}+\frac{(1-p)}{p}{\displaystyle \int {\left|(d^{m}f(x)/d{x}^m)\right|}^{2}dx},0\le p\le 1$$E22

where, $$\lambda =(1-p)/p$$($$0\le \lambda \le +\infty$$) is a Lagrangian parameter, n is the number of observations,$$W_j$$ is weight factor at point$$x_j$$and $$m$$is the derivative order.

It can be shown that spline of degree$$k=2m-1$$, having $$2m-2$$ continuous derivatives, is an optimal solution; where, $$n\ge 2m$$. In this chapter the cubic smoothing spline is chosen to have a minimum curvature property; hence, $$m=2$$ ($$2m-1=3$$) and $$f\in C^2[x_1,x_n]$$[80-82].

According to this formula, the natural cubic spline interpolation is obtained by $$p=1$$ and the least-squares straight line fit by$$p=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$$\left\{W_j\right\}=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. $$W_j\ge 1$$. The values of $$\left\{W_i\right\}$$ and $$p$$ can be constant or variable in$$\{(x_i,y_i)\}$$sequence [85]. Here, the constant weights and smoothing parameter are studied. The $$\left\{W_j\right\}$$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/2^8$$spatial 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 $$u_x(x,z,t=0)$$ which is: $$u_z(x,z,t=0)=\exp (-500({(x-0.35)}^2+{(z-0.25)}^2))$$.In the numerical simulation it is assumed that: $$p=0.85$$ and$$\epsilon ={10}^{-5}$$. As mentioned before, the absorbing boundary condition is considered explicitly for simulation of infinite boundaries.This modification, performed for P-SV wave equations, is:

$$\begin{array}{c}\left((\lambda +2\mu )u_{x,xx}+\mu u_{x,zz}\right)+\left((\lambda +\mu )u_{z,xz}\right)=\rho ({\ddot{u}}_x+Q(x,z).{\dot{u}}_x)\\ \left((\lambda +2\mu )u_{z,zz}+\mu u_{z,xx}\right)+\left((\lambda +\mu )u_{x,xz}\right)=\rho ({\ddot{u}}_z+Q(x,z).{\dot{u}}_z)\end{array}$$E23

In the above equation, it is assumed that: $$Q=a_x(e^{b_x.x^2}+e^{b_x.{(1-x)}^2})+a_z(e^{b_z.{(1-z)}^2})$$, where $$a_x=a_z=30$$, and$$b_x=-110,b_z=-70$$. The free boundary is imposed by equivalent force in the free surface boundary [48].

The snapshots of solutions $$u_x$$ and $$u_z$$ 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 solution$$u_x$$at times: a)0.298, b)0.502, c) 0.658, d) 0.886.

Figure 10.

Snapshots of solution $$u_z$$ 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.