- Email: [email protected]
Wavelet based schemes for linear advection–dispersion equation
Applied Mathematics and Computation 218 (2011) 3786–3798
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Wavelet based schemes for linear advection–dispersion equation q K. Sandeep a,⇑, Shikha Gaur b, D. Dutta c, H.S. Kushwaha c a b c
Department of Mechanical Engineering, Institute of Technology, BHU, Varanasi, India Department of Applied Mathematics, Institute of Technology, BHU, Varanasi, India Bhabha Atomic Research Centre, Trombay, Mumbai, India
a r t i c l e
i n f o
Keywords: Wavelet-Galerkin Advection–dispersion Adaptive solution Linear B-spline wavelets Multiscale decomposition
a b s t r a c t In this paper, two wavelet based adaptive solvers are developed for linear advection– dispersion equation. The localization properties and multilevel structure of the wavelets in the physical space are used for adaptive computational methods for solution of equation which exhibit both smooth and shock-like behaviour. The first framework is based on wavelet-Galerkin and the second is based on multiscale decomposition of finite element method. Coiflet wavelet filter is incorporated in both the methods. The main advantage of both the adaptive methods is the elimination of spurious oscillations at very high Peclet number. 2011 Elsevier Inc. All rights reserved.
1. Introduction Many interesting physical systems are characterized by the presence of localized structure or sharp transition, which might occur anywhere in the domain or change their locations in space with time. Popular methods such as finite element, so-called meshless and recently developed wavelet methods, to solve these problems efficiently, use adaptive grid techniques. Adaptive refinement techniques can also be profitably applied in solving partial differential equations useful in many applications, including simulation, animation, computer vision, etc. The currently existing adaptive grid techniques may be roughly classified as either subdivision schemes or basis refinement techniques. The major difference between these approaches is that subdivision schemes solve problems in the physical space by increasing the nodes while basis refinement techniques (including hierarchical basis in finite element method) solve problems in coefficient space. Though both the adaptive grid techniques are well understood, a lot of work has to be done for efficient implementation in complex domain, in particular to reduce computational time. Wavelet has high potential for fast, hierarchical and locally adaptive algorithms because of their compactly supported refinable basis functions [1–4]. Liandrat and Tchamitchian [5] proposed the first algorithm based on a spatial approximation exploiting the regularity properties of an orthonormal wavelet basis. Beylkin and Keiser [6] used wavelet expansion for adaptively updating numerical solution of nonlinear partial differential equations, which exhibit both smooth and shock-like behaviour. Due to the signal processing base of traditional wavelet, the research in PDE simulations [7,8] was limited to simple domain and boundary conditions. This limitation has been eliminated with the development of the lifting scheme [9] and stable completion [10,11]. By using lifting scheme, Vasilyev and Paolucci [12] developed wavelet collocation method to adapt computational refinements to local demands of the solution. Krysl et al. [13] developed conforming hierarchical adaptive refinement methods where hierarchical refinement treats refinement as the addition of finer level ‘‘detail function’’ to an
q
This study is supported by the DAE, BRNS, Mumbai, India (2007/36/74-BRNS/2688) and DST, India Grant no. SR/S4/MS:367/06.
⇑ Corresponding author.
E-mail address: [email protected] (K. Sandeep). 0096-3003/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.09.023
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3787
unchanged set of coarse-level functions. Amaratunga and Sudarshan [14] customized the second-generation wavelets to generate hierarchical basis for finite element method to solve PDEs both hierarchically and adaptively. Advection–dispersion equation exhibits discontinuity (shocks) after a finite time. Further, the numerical solution shows spurious oscillations when dispersion coefficient is small as compared to the velocity of flow, i.e., at high Peclet number. To capture the singular effects in the solution, the domain would require very fine resolution near singularities. The classical discretization based on uniform grid will be highly uneconomical. The existing numerical techniques use artificial dispersion to overcome the stability problem. A comparison of various methods is discussed by Johnson [15], and Zienkiewicz and Taylor [16]. Main interest of our work is to remove numerical instability by adaptive grid generation (very fine grid at the critical zone) and show the effect of filtering by using wavelets without sacrificing the accuracy of the results. In the paper, two wavelet based methods are presented and the results are compared with some recent finite difference methods. In the first method, linear advection–dispersion equation is solved by using wavelet-Galerkin method. For calculation of inner product, Newton–cotes method is used which can be replaced by recently developed highly efficient methods [17,18]. The basic idea behind the adaptive solution is simply based on the analysis of wavelet coefficients, which gives information about the region where sharp change starts or ends. At any time step only local matrix reflecting the local changes in the solution, is solved. The method uses efficient data structure of uniform grid and periodic basis function to evaluate the entries of the stiffness matrix. In the second method, the finest scale finite element solution space is projected onto the scaling and wavelet spaces resulting in the decomposition of high- and low-scale components. Repetition of such a projection results in multi-scale decomposition of the fine scale solution. In the proposed wavelet projection method, the fine scale solution can be obtained by any other numerical method also. Subsequently the properties of the wavelet functions are exploited to eliminate the nodes from the smooth region where the wavelet coefficients will not exceed a preset tolerance. This wavelet-based multi-scale transformation hierarchically filters out the less significant part of the solution, and thus provides an effective framework for the selection of significant part of the solution. In this process, the ‘big’ coefficient matrix at the finest level will be calculated once for complete domain whereas the ‘small’ adaptively compressed coefficient matrix for a priory known localized dynamic zone of high gradient, which will be considerably less expensive to solve, will be used for the solution in every step of the solution. Similar technique is used in the software QUADFLOW [19] using finite volume method. The paper presents simple, general methods with minimal mathematical framework. The present methods remove a number of implementation headache associated with adaptive grid techniques and is a general technique, independent of domain dimension. These methods have very important and highly practical consequence because they reduce the computational time significantly. Description of the different element of the algorithm in combination with different mathematical comments on the method, are provided. The resulting algorithms, while capturing full generality of methods, are surprisingly simple. A set of concrete, compelling examples based on our implementation is also the contribution of the paper. The rest of the paper is organized as follows. Advection–dispersion equation is presented in Section 2. Multiscaling using wavelet is briefly discussed in Section 3. In Section 4, Method-I, i.e. wavelet Galerkin method is discussed in detail. Multiscale decomposition of finite element, i.e. Method-II is discussed in Section 5. Results obtained using these two methods and their comparative study is presented in Section 6. Finally, Section 7 contains conclusion. 2. Advection–dispersion equation Generally pollutant concentration in atmosphere is governed by advection–dispersion equation. The one dimensional advection–dispersion equation can be written as:
@C @C @2C ¼ u þ D 2 kC; @t @x @x
ð1Þ
where C is concentration (mg/l), t duration (days), u the flow velocity (m/day), x is the distance along the direction of flow from the upstream boundary of modeled domain (m), D the dispersive coefficient (m2/day), and k the decay constant (day1). We are considering following boundary and initial condition:
at t P 0;
x ¼ 0;
C ¼ C 0 ðtÞ; @C ¼ 0; at t P 0; x ¼ L; @x at t ¼ 0; 0 < x 6 L; C ¼ 0;
ð2aÞ ð2bÞ ð2cÞ
where C0 is the concentration of constant magnitude. Applying weighted residual method in the advection–dispersion Eq. (1), we get
Z 0
L
! @C @C @2C þu D 2 þ kC dx ¼ 0: w @t @x @x
ð3Þ
3788
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Applying the natural boundary conditions in the weak form, we get
Z 0
L
@w @C @C þ kwC þ w dx ¼ 0: wu þ D @x @x @t
ð4Þ
The mathematical difficulties that arise from advection–dispersion problems are due to the fact that the variational formulation has a coercive constant which is extremely small when compared to the continuity constant. On the discrete level, the problem shows unstable behaviour. Adding stabilization terms may be viewed as a way of improving the coercive constant [20]. In other words, in case of high Peclet number Pe = (uDx)/D, i.e., when the dispersion coefficient D is small compared to characteristic velocity u of the flow, the numerical methods which do not properly take into account this aspect lead to the numerical instabilities, resulting in loss of both qualitative and quantitative accuracy. Assuming element size Dx to be very small and uniform in the whole domain is not a practical solution. Many different methods based on upwind scheme and artificial diffusion are developed to suppress spurious numerical oscillations. On the other hand, our wavelet based two different methods can be interpreted as augmenting the test space by a suitable supplementary space similar to those developed by Brezzi et al. [21]. 3. Multiscaling using wavelet transforms The concept of multiresolution analysis is to interpolate an unknown field at a coarse level by means of so-called scaling functions. Any improvement to the initial approximation consists in adding ‘details’ wherever required provided by new functions known as wavelets. Thus it is well suited for multiscale solutions. A multi-resolution analysis is a nested sequence
V 0 V 1 V 2 . . . . . . L2 ðRÞ
ð5aÞ
satisfying the following properties:
Vj Vj+1, f ðxÞ 2 V j () f ð2xÞ 2 V jþ1 , f ðxÞ 2 V 0 () f ðx þ 1Þ 2 V 0 , [V j is dense in L2(R);\ Vj = {0}.
Each subspace Vj is spanned by a set of scaling function {/j, k(x), "k 2 Z}. The complement of Vj in Vj+1 is defined as subspace Wj, i.e.,
V jþ1 ¼ V j W j
8j 2 Z;
ð5bÞ
The space Vj+1 can be decomposed in a consecutive manner
V jþ1 ¼ V 0 W 0 W 1 W 2 . . . W j :
ð5cÞ
The basis functions in Wj are called wavelet functions and are denoted by wj,k. These wavelet and scaling functions in different scales are used for wavelet-based multi-scaling. A function f 2 L2(R) is approximated by its projection Pjf onto the space Vj 1 X
Pj f ¼
ck /j;k :
ð6aÞ
k¼1
Let us denote the projection of f on Wj as Qjf. Then we have
Pj f ¼ Pj1 f þ Q j1 f :
ð6bÞ
In the multi-scaling, scaling coefficients cj,k are decomposed into the scaling coefficients cj1,k of the approximation Pj1f and wavelet coefficients dj1,k of Qj1f at the next coarser scale. 4. Method-I Wavelets and multiscale may act as the natural functional framework for stabilized advection–dispersion problems [20]. In this paper, compactly supported piecewise polynomials B-spline wavelets, which can be generated by convolution as well as recursive relation proposed by Cox and Carl de Boor [22], are used. The constant B-spline over the ith subinterval is defined as
/1k ðxÞ ¼
1 tk 6 x < t kþ1 0
otherwise
;
ð7aÞ
where subscript j denoting multiscale level is dropped for convenience. The mth order B-spline can be recursively calculated by using a formula of Cox- de Boor
/m k ðxÞ ¼
x tk t kþm x m1 /m1 ðxÞ þ / ðxÞ: t kþm t kþ1 kþ1 t kþm1 tk k
ð7bÞ
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3789
With equally spaced ‘‘knots tk’’ of unit length, (7) can be simplified greatly because the B-spline blending function is same for all the segments. The support of the basis function goes up one point for each increase in order m i.e., for any positive integer m, sup[/m] = [0, m]. In the uniform knot vector case, the functions are successive convolutions of the ‘‘box’’ function. The first order spline is defined as a box function /1 = v[0, 1)(x). The mth order spline can be obtained by the convolution
/m ðxÞ ¼ ð/m1 /1 ÞðxÞ ¼
Z
1
/m1 ðx tÞdt;
ð8Þ
0
where ⁄ denotes the convolution operator. Chui and Wang [23] and Unser et al. [24] developed independently B-spline wavelets for multiresolution analysis. The wavelet function for mth order B-spline
wm ¼
2m2 X
1 2
ð1Þi /2m ði þ 1ÞU2m m ð2x iÞ
m1
ð9aÞ
i¼0
with support [0, 2 m 1], is a basic wavelet that generates wavelet function, and
U2m m ð2xÞ ¼
m X
ð1Þi
i¼0
m i
/m ðx iÞ:
ð9bÞ
The 2nd order B-spline and its corresponding wavelet function are shown in Fig. 1. On each level of resolution j an unknown function f(x) can be approximated by using B-spline scaling function
fj ðxÞ ¼
X
aj;k /m ð2j x kÞ ¼
X
k
aj;k /m j;k ðxÞ:
ð10Þ
k
We can define the detail function dj (x) as
dj ðxÞ ¼
X
X
bj;k wm ð2j x kÞ ¼
k
bj;k wm j;k ðxÞ;
ð11Þ
k
where aj,k and bj,k are the scaling and wavelet coefficients at level ‘‘j’’ that defines special resolution of the approximations. The following relationship holds for the function at next higher resolution:
fjþ1 ðxÞ ¼ fj ðxÞ þ dj ðxÞ
ð12aÞ
or
fjþ1 ðxÞ ¼
X
X
ajþ1;k /m jþ1;k ðxÞ ¼
k
aj;k /m j;k ðxÞ þ
k
X
bj;k wm j;k ðxÞ:
ð12bÞ
k
If we start with the level j0 scaling function, then the approximation of level ‘‘j + 1’’ is
fjþ1 ðxÞ ¼
X k
aj0;k /m j0 ;k ðxÞ þ
j X X i¼j0
bi;k wm i;k ðxÞ:
ð13Þ
k
Similarly, concentration function can be expressed as the linear combination of scaling and wavelet functions
Cðx; tÞ ¼
X k
aj0;k ðtÞ/m j0 ;k ðxÞ þ
j X X i¼j0
bi;k ðtÞwm i;k ðxÞ:
k
Fig. 1. (a) Second order B-spline basis function, (b) corresponding wavelet.
ð14Þ
3790
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Substituting this equation in the weak form (4) and using the weight function w as mth order B-spline /m and its wavelet functions wm, we get
_ ¼ 0: ½KfCg þ ½MfCg
ð15Þ
Considering the level of resolution from jo to j, the stiffness and mass matrices of the dynamic Eq. (15) are
2
K/jo /jo
K/jo wjo
6K 6 wjo/jo 6 ½K ¼ 6 6 6 4
K/jo wj 7 7 7 7 7 7 5 Kwj wj
Kwjo wjo
Kwj /jo
3
: K/jo wj
Kwj wjo
ð16aÞ
and
2
M/jo /jo 6M 6 wjo/jo 6 ½M ¼ 6 6 6 4 Mwj /jo
M/jo wjo Mwjo wjo Mwj wjo
M/jo wj
3
M/jo wj 7 7 7 7 7; 7 5 Mw j w j
ð16bÞ
where the typical elements of stiffness and mass matrices at any resolution are
2RL 6 6 K/jo /jo ¼ 6 6 4
0
ð/jo;1
RL
@/jo;1 @x
@/ ð/jo;n @xjo;1 0
þD
þ
@/jo;1 @/jo;1 @x @x
þ k/jo;1 /jo;1 Þdx
@/ @/ D @xjo;n @xjo;1
þ k/jo;1 /jo;1 Þdx : :
RL 0
ð/jo;1
@/jo;n @x
þD
@/jo;1 @/jo;n @x @x
þ k/jo;1 /jo;n Þdx
RL
@/ ð/jo;n @xjo;n 0
þ
@/ @/ D @xjo;n @xjo;n
3 7 7 7 7 5
ð17aÞ
þ k/jo;1 /jo;n Þdx
and
2RL 6 6 M/jo wi ¼ 6 6 4
0
RL 0
/jo;1 wi;1 dx
/jo;n wi;1 dx
RL 0
/jo;1 wi;n dx
RL
00
/jo;n wi;n dx
3 7 7 7: 7 5
ð17bÞ
Similarly, the stiffness and mass matrices at other resolution can be obtained. The unknown coefficients
3 a_ jo 7 6 6b 7 6_ 7 6 jo 7 _ ¼ 6 bjo 7 fCg ¼ 6 7 and fCg 7 6 4 5 5 4 _ bj bj 2
ajo
3
2
ð18Þ
at time ts+1 can be obtained from the coefficients at ts by using the recursive relation of a-family of approximation
f½M þ Dt a½KgfCgsþ1 ¼ f½M ð1 aÞDt½KgfCgs :
ð19Þ
Algorithm Outline The algorithm that uses explicitly the scaling and wavelet basis stands on the fact that the large wavelet coefficient represent high gradient region on the physical domain. The problem can be solved by inverting the M at each scale once and for all at the beginning of the calculation. Then the complete solution can be obtained by matrix multiplication at each time advancing step. On the other hand, our main goal is to compute scaling and wavelet coefficients in a region as close as possible to the space where sharp transition occurs. The time marching version of our algorithm is designed in order to fulfill efficiently this requirement. It involves the following steps.
Fig. 2. (a) Vj, (b) Vj1 Wj1, (c) Vj2 Wj2 Wj1.
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3791
Fig. 3. Solution of the linear advection–dispersion equation using Method I.
Fig. 4. Evolution, decay and marching of wavelet coefficients.
Flowchart for construction of wavelet-Galerkin method. i at ith node. 1. Discretize scaling function /i and its derivative @/ @x i 2. Discretize wavelet wi and its derivative @w at ith node. @x i i 3. Translate periodic /i ; @/ ; wi and @w and use pointwise product and Newton cotes integration to calculate elements of @x @x stiffness matrix K and mass matrix M. 4. Select first few nodes to form ‘‘active’’ space. 5. Loop over the time domain a. Calculate scaling coefficient by solving dynamic Eq. (19) in the active space b. Use wavelet spectrum to estimate active space for the next time step. 6. End loop 7. Use Coiflet filter to eliminate spurious oscillation
3792
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 5. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
Fig. 6. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
5. Method-II In this variational multiscale stabilization method, the problem is discretized at the finest scale and then forward transformation is applied for multiscale projection. The standard Galerkin finite element method employing simple polynomialbased elements is used for ease of calculations of inner products. Similar approach has been used by Brezzi et al. [21], who have accounted fine scale in order to accurately calculate the coarse scale without ‘seeing’ the fine scale feature. Let the finite element equation at the finest resolution J be
AJ C J ¼ fJ :
ð20Þ j
j
Using B-spline wavelets [25] for forward transformation [P jQ ] = TJ, the finest scale is replaced by multiscale as shown in the Fig. 2, i.e.,
T J1
C J1 ¼ CJ : dJ1
ð21Þ
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3793
Fig. 7. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
Fig. 8. Comparison of Method I with analytical solution.
Substituting (21) in (20), we get
AJ T J1
C J1 ¼ fJ dJ1
ð22Þ
or
T TJ1 AJ T J1
C J1 ¼ T TJ1 fJ ; dJ1
ð23Þ
where CJ1 are the scaling coefficients, dJ1 are the wavelets coefficients at the lower resolution. The wavelet transform can be continued until a desired coarsest level is achieved. The modified stiffness matrix can be expressed as
T TJ1 AJ T J1 ¼ AJ1 :
ð24Þ
The objective of multiscale is to eliminate small detail coefficients and develop the compressed stiffness matrix from the finest scales solution which will be considerably less expensive to solve. It will also eliminate the process of re-meshing the FE domain and re-computing the stiffness matrix with change in location of sudden jump. Multiscale is performed by
3794
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 9. Comparison of Method I, Method II and analytical solution.
Fig. 10. Concentration using Method I before filtering and after fifth level filtering.
projecting the solution at the fine scale space Vj onto next coarser levels Wj1 and Vj1, recursively. It is well established fact that the coefficients of wavelet space are negligible in the smooth region. Retaining grid points associated with wavelets coefficients at high gradient region and eliminating them from a priori known smooth region will generate adaptive grid. In other words, at any time, the computational grid should include points near the sharp gradient zone. Computational grid should also consist of grid points associated with wavelets whose coefficients can possibly become significant during the next time step. Finite element stiffness matrix at the finest scale will be transformed to the coarsest scale only once, i.e. there is no need of re-meshing and re-computing the stiffness matrix with change of sharp gradient in space and time. Compression of stiffness matrix by eliminating the small wavelet coefficients from the smooth region will be required in all subsequent steps of adaptive solution. Fig. 2 depicts wavelet transformation. Coiflet filter is included to eliminate high frequency and small magnitude oscillations which are typical numerical instability phenomena near jump at high Peclet number.
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3795
Fig. 11. Concentration using Method II for coarse grid before filtering and after seventh level filtering.
Fig. 12. Concentration using Method II for fine grid before filtering and after seventh level filtering.
6. Results and discussions To demonstrate the algorithm in our experiment we begin with the simple example of linear advection–dispersion, which can be viewed as a simple case of complex three-dimensional pollution problem. Mostly advection–dispersion solvers use upwind scheme to reduce numerical instability. In the present methods, wavelet based adaptive grid and filtering are used for accuracy and numerical stability. Simple time marching algorithm (a-family) is employed only on the small set of active nodes. The decay constant k = 0 is assumed in this analysis. The change in concentration with time for Peclet number Pe = 50 at locations x = 200 m and 500 m are shown in Fig. 3 and corresponding wavelet coefficients are shown in Fig. 4. It can be observed that changes in wavelet coefficients are very small compared to the scaling coefficients and start well before the sudden jump. This helps in automatic grid refinement with time at the location of sudden jump. The results of Method-I are compared with the analytical method and method based on upwind scheme [26] in Figs. 5–7 for Pe = 10, 100,
3796
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 13. Detail coefficients from first to eighth level for Method II.
200, respectively. It can be observed that the solution by Method-I is in good agreement with the analytical solution even at high Peclet number. However, spurious oscillations at very high Peclet number Pe = 500 are clearly visible and shown in Fig. 8. Newton–cotes method is used for calculation of inner product which makes Method-I inefficient and it cannot be extended in two-dimensional advection–dispersion problems. Therefore, it is necessary to use numerical techniques proposed by Beylkin [17] and Dahmen and Micchelli [18] and will be used in future work. Method-II utilizes the stiffness matrix developed by finite element method. The stiffness matrix is projected into wavelet space and obtained multiscale structure. A comparison of Method-II with Method-I and analytical solution is shown in Fig. 9. It can be observed that Method-II is not as stable as Method-I. However, Method-II is more efficient than Method-I. MATLAB programs developed for calculation of stiffness matrix for Method- I takes 4.35 s while Method-II takes only 1.13 s. Fig. 10 shows the response of Method-I at Pe = 1 and the result is improved by using Coiflet filters. It can be observed that the filtering reduces the spurious oscillations, significantly. Similar filtering at Pe = 1 in Method-II is applied and the results before filtering and after filtering are shown in Fig. 11. Adaptive grid refinement improves the results and it is shown in Fig. 12. A typical pattern of wavelet coefficients which are eliminated during filtering in the preceding case at various levels are shown in Fig. 13. The elimination of wavelet coefficient up to level VII improves the result which is presented in Fig. 12. However, deterioration of the result is
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3797
Fig. 14. Concentration using Method II before filtering and after eighth level filtering.
clearly visible in Fig. 14 due to unjustifiable elimination of wavelet coefficients of level VIII which have large magnitude and low frequency. 7. Conclusion In case of hyperbolic PDEs which show jump discontinuity, the standard numerical schemes such as finite element and finite difference will in general exhibit spurious oscillations. To overcome this difficulty, two multiscale wavelet based methods are proposed and wavelet filtering is used. Results are compared with analytical and other numerical solutions. It has been observed that after filtering, Method-II is as good as Method-I. Generation of stiffness matrix in Method-I is more expensive than in Method-II. Therefore, a customized numerical technique in Method-I should be employed to improve its efficiency. Moreover, the original spur for reduction of degree of freedom, which will be very useful for higher dimensional problems, can be easily achieved in Method-II. References [1] Y. Meyer, Wavelets with compact support, in: U. Chicago: Zygmund Lectures, 1987. [2] S.G. Mallat, A theory for multiresoltion signal decomposition: The wavelet representation, Communications on Pure and Applied Mathematics 41 (1988) 674–693. [3] I. Daubechies, Orthonormal bases of compactly supported wavelets, Communications on Pure and Applied Mathematics 41 (1988) 909–996. [4] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, MA, 1996. [5] J. Liandrat and Ph. Tchamitchian, Resolution of the 1D regularized Burgers equation using a spatial wavelet approximation, Report No: NASA CR – 187480, NASA Langley Research Centre, Hamptonva, (1990). [6] G. Beylkin, J.M. Keiser, On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases, Journal of Computational Physics 132 (1997) 233–259. [7] K. Amaratunga, J.R. Williams, Wavelet based Green’s function approach to 2D PDEs, Engineering Computations 10 (4) (1993) 349–367. [8] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, Journal of Computational Physics 106 (1993) 155–175. [9] W. Swelden, The lifting scheme: A construction of second generation wavelets, SIAM Journal of Mathematical Analysis 29 (2) (1998) 511–546. [10] J. Carnicer, W. Dahmen, J. Pena, Local decompositions of refinable spaces, Applied Computational and Harmonic Analysis 3 (1996) 125–153. [11] W. Dahmen, R. Stevenson, Element-by-element construction of wavelets satisfying stability and moment conditions, SIAM on Numerical Analysis 37 (1999) 319–352. [12] O.V. Vasilyev, S. Paolucci, A dynamically adaptive multilevel wavelet collocation method for solving partial differential equations in finite domain, Journal of Computational Physics 125 (1996) 498–512. [13] P. Krysl, E. Grinspun, P. Schroder, Natural hierarchical refinement for finite element methods, International Journal for Numerical Methods in Engineering 56 (8) (2002) 1109–1124. [14] K. Amaratunga, R. Sudarshan, Multiresolution modeling with operator-customized wavelet derived from finite elements, Computer Methods in Applied Mechanics and Engineering 195 (2006) 2509–2532. [15] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, 1987. [16] O.C. Zienkiewicz, R.L. Taylor, Finite Element Method Volume 3: Fluid Dynamics, Fifth ed., Butterworth-Heinemann, 2000. [17] G. Beylkin, On the representation of operators in bases of compactly supported wavelets, SIAM Journal on Numerical Analysis 6 (6) (1992) 1716–1740. [18] W. Dahmen, C.A. Micchelli, Using the refinement equation for evaluating integrals of wavelets, SIAM J. Numerical Analysis 30 (1993) 507–537. [19] A. Cohen, S.M. Kaber, S. Muller, M. Postel, Fully adaptive multiresolution finite volume schemes for conservation laws, Mathematics of Computation 24 (2001) 183–225.
3798 [20] [21] [22] [23]
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
S. Bertoluzza, C. Canuto, A. Tabacco, Negative norm R stabilization of convection–diffusion problems, Applied Mathematics Letter 13 (2000) 121–127. F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, b ¼ g, Computer Methods in Applied Mechanics and Engineering 145 (1997) 329–340. D.F. Rogers, J.A. Adams, Mathematical Elements for Computer Graphics, McGraw-Hill Publishing Company, Singapore, 1990. C.K. Chui, J.Z. Wang, On compacty supported spline wavelets and a duality principle, Transaction of the American Mathematical Society 330 (2) (1992) 903–915. [24] M. Unser, A. Aldroubi, M. Eden, A family of polynomial spline wavelet transforms, NCRR rep.153/90, Nat. Inst. Health, (1990). [25] E.J. Stollnitz, T.D. Derose, D.H. Salesin, Wavelets for Computer Graphics, Morgan Kaufmann Pub. Inc., California, 1996. 209–213. [26] Li. Guoyuan, C.R. Jackson, Simple, accurate, and efficient revisions to MacCormack and Saulyev schemes: High Peclet numbers, Applied Mathematics and Computation 186 (2007) 610–622.
Contents lists available at SciVerse ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Wavelet based schemes for linear advection–dispersion equation q K. Sandeep a,⇑, Shikha Gaur b, D. Dutta c, H.S. Kushwaha c a b c
Department of Mechanical Engineering, Institute of Technology, BHU, Varanasi, India Department of Applied Mathematics, Institute of Technology, BHU, Varanasi, India Bhabha Atomic Research Centre, Trombay, Mumbai, India
a r t i c l e
i n f o
Keywords: Wavelet-Galerkin Advection–dispersion Adaptive solution Linear B-spline wavelets Multiscale decomposition
a b s t r a c t In this paper, two wavelet based adaptive solvers are developed for linear advection– dispersion equation. The localization properties and multilevel structure of the wavelets in the physical space are used for adaptive computational methods for solution of equation which exhibit both smooth and shock-like behaviour. The first framework is based on wavelet-Galerkin and the second is based on multiscale decomposition of finite element method. Coiflet wavelet filter is incorporated in both the methods. The main advantage of both the adaptive methods is the elimination of spurious oscillations at very high Peclet number. 2011 Elsevier Inc. All rights reserved.
1. Introduction Many interesting physical systems are characterized by the presence of localized structure or sharp transition, which might occur anywhere in the domain or change their locations in space with time. Popular methods such as finite element, so-called meshless and recently developed wavelet methods, to solve these problems efficiently, use adaptive grid techniques. Adaptive refinement techniques can also be profitably applied in solving partial differential equations useful in many applications, including simulation, animation, computer vision, etc. The currently existing adaptive grid techniques may be roughly classified as either subdivision schemes or basis refinement techniques. The major difference between these approaches is that subdivision schemes solve problems in the physical space by increasing the nodes while basis refinement techniques (including hierarchical basis in finite element method) solve problems in coefficient space. Though both the adaptive grid techniques are well understood, a lot of work has to be done for efficient implementation in complex domain, in particular to reduce computational time. Wavelet has high potential for fast, hierarchical and locally adaptive algorithms because of their compactly supported refinable basis functions [1–4]. Liandrat and Tchamitchian [5] proposed the first algorithm based on a spatial approximation exploiting the regularity properties of an orthonormal wavelet basis. Beylkin and Keiser [6] used wavelet expansion for adaptively updating numerical solution of nonlinear partial differential equations, which exhibit both smooth and shock-like behaviour. Due to the signal processing base of traditional wavelet, the research in PDE simulations [7,8] was limited to simple domain and boundary conditions. This limitation has been eliminated with the development of the lifting scheme [9] and stable completion [10,11]. By using lifting scheme, Vasilyev and Paolucci [12] developed wavelet collocation method to adapt computational refinements to local demands of the solution. Krysl et al. [13] developed conforming hierarchical adaptive refinement methods where hierarchical refinement treats refinement as the addition of finer level ‘‘detail function’’ to an
q
This study is supported by the DAE, BRNS, Mumbai, India (2007/36/74-BRNS/2688) and DST, India Grant no. SR/S4/MS:367/06.
⇑ Corresponding author.
E-mail address: [email protected] (K. Sandeep). 0096-3003/$ - see front matter 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2011.09.023
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3787
unchanged set of coarse-level functions. Amaratunga and Sudarshan [14] customized the second-generation wavelets to generate hierarchical basis for finite element method to solve PDEs both hierarchically and adaptively. Advection–dispersion equation exhibits discontinuity (shocks) after a finite time. Further, the numerical solution shows spurious oscillations when dispersion coefficient is small as compared to the velocity of flow, i.e., at high Peclet number. To capture the singular effects in the solution, the domain would require very fine resolution near singularities. The classical discretization based on uniform grid will be highly uneconomical. The existing numerical techniques use artificial dispersion to overcome the stability problem. A comparison of various methods is discussed by Johnson [15], and Zienkiewicz and Taylor [16]. Main interest of our work is to remove numerical instability by adaptive grid generation (very fine grid at the critical zone) and show the effect of filtering by using wavelets without sacrificing the accuracy of the results. In the paper, two wavelet based methods are presented and the results are compared with some recent finite difference methods. In the first method, linear advection–dispersion equation is solved by using wavelet-Galerkin method. For calculation of inner product, Newton–cotes method is used which can be replaced by recently developed highly efficient methods [17,18]. The basic idea behind the adaptive solution is simply based on the analysis of wavelet coefficients, which gives information about the region where sharp change starts or ends. At any time step only local matrix reflecting the local changes in the solution, is solved. The method uses efficient data structure of uniform grid and periodic basis function to evaluate the entries of the stiffness matrix. In the second method, the finest scale finite element solution space is projected onto the scaling and wavelet spaces resulting in the decomposition of high- and low-scale components. Repetition of such a projection results in multi-scale decomposition of the fine scale solution. In the proposed wavelet projection method, the fine scale solution can be obtained by any other numerical method also. Subsequently the properties of the wavelet functions are exploited to eliminate the nodes from the smooth region where the wavelet coefficients will not exceed a preset tolerance. This wavelet-based multi-scale transformation hierarchically filters out the less significant part of the solution, and thus provides an effective framework for the selection of significant part of the solution. In this process, the ‘big’ coefficient matrix at the finest level will be calculated once for complete domain whereas the ‘small’ adaptively compressed coefficient matrix for a priory known localized dynamic zone of high gradient, which will be considerably less expensive to solve, will be used for the solution in every step of the solution. Similar technique is used in the software QUADFLOW [19] using finite volume method. The paper presents simple, general methods with minimal mathematical framework. The present methods remove a number of implementation headache associated with adaptive grid techniques and is a general technique, independent of domain dimension. These methods have very important and highly practical consequence because they reduce the computational time significantly. Description of the different element of the algorithm in combination with different mathematical comments on the method, are provided. The resulting algorithms, while capturing full generality of methods, are surprisingly simple. A set of concrete, compelling examples based on our implementation is also the contribution of the paper. The rest of the paper is organized as follows. Advection–dispersion equation is presented in Section 2. Multiscaling using wavelet is briefly discussed in Section 3. In Section 4, Method-I, i.e. wavelet Galerkin method is discussed in detail. Multiscale decomposition of finite element, i.e. Method-II is discussed in Section 5. Results obtained using these two methods and their comparative study is presented in Section 6. Finally, Section 7 contains conclusion. 2. Advection–dispersion equation Generally pollutant concentration in atmosphere is governed by advection–dispersion equation. The one dimensional advection–dispersion equation can be written as:
@C @C @2C ¼ u þ D 2 kC; @t @x @x
ð1Þ
where C is concentration (mg/l), t duration (days), u the flow velocity (m/day), x is the distance along the direction of flow from the upstream boundary of modeled domain (m), D the dispersive coefficient (m2/day), and k the decay constant (day1). We are considering following boundary and initial condition:
at t P 0;
x ¼ 0;
C ¼ C 0 ðtÞ; @C ¼ 0; at t P 0; x ¼ L; @x at t ¼ 0; 0 < x 6 L; C ¼ 0;
ð2aÞ ð2bÞ ð2cÞ
where C0 is the concentration of constant magnitude. Applying weighted residual method in the advection–dispersion Eq. (1), we get
Z 0
L
! @C @C @2C þu D 2 þ kC dx ¼ 0: w @t @x @x
ð3Þ
3788
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Applying the natural boundary conditions in the weak form, we get
Z 0
L
@w @C @C þ kwC þ w dx ¼ 0: wu þ D @x @x @t
ð4Þ
The mathematical difficulties that arise from advection–dispersion problems are due to the fact that the variational formulation has a coercive constant which is extremely small when compared to the continuity constant. On the discrete level, the problem shows unstable behaviour. Adding stabilization terms may be viewed as a way of improving the coercive constant [20]. In other words, in case of high Peclet number Pe = (uDx)/D, i.e., when the dispersion coefficient D is small compared to characteristic velocity u of the flow, the numerical methods which do not properly take into account this aspect lead to the numerical instabilities, resulting in loss of both qualitative and quantitative accuracy. Assuming element size Dx to be very small and uniform in the whole domain is not a practical solution. Many different methods based on upwind scheme and artificial diffusion are developed to suppress spurious numerical oscillations. On the other hand, our wavelet based two different methods can be interpreted as augmenting the test space by a suitable supplementary space similar to those developed by Brezzi et al. [21]. 3. Multiscaling using wavelet transforms The concept of multiresolution analysis is to interpolate an unknown field at a coarse level by means of so-called scaling functions. Any improvement to the initial approximation consists in adding ‘details’ wherever required provided by new functions known as wavelets. Thus it is well suited for multiscale solutions. A multi-resolution analysis is a nested sequence
V 0 V 1 V 2 . . . . . . L2 ðRÞ
ð5aÞ
satisfying the following properties:
Vj Vj+1, f ðxÞ 2 V j () f ð2xÞ 2 V jþ1 , f ðxÞ 2 V 0 () f ðx þ 1Þ 2 V 0 , [V j is dense in L2(R);\ Vj = {0}.
Each subspace Vj is spanned by a set of scaling function {/j, k(x), "k 2 Z}. The complement of Vj in Vj+1 is defined as subspace Wj, i.e.,
V jþ1 ¼ V j W j
8j 2 Z;
ð5bÞ
The space Vj+1 can be decomposed in a consecutive manner
V jþ1 ¼ V 0 W 0 W 1 W 2 . . . W j :
ð5cÞ
The basis functions in Wj are called wavelet functions and are denoted by wj,k. These wavelet and scaling functions in different scales are used for wavelet-based multi-scaling. A function f 2 L2(R) is approximated by its projection Pjf onto the space Vj 1 X
Pj f ¼
ck /j;k :
ð6aÞ
k¼1
Let us denote the projection of f on Wj as Qjf. Then we have
Pj f ¼ Pj1 f þ Q j1 f :
ð6bÞ
In the multi-scaling, scaling coefficients cj,k are decomposed into the scaling coefficients cj1,k of the approximation Pj1f and wavelet coefficients dj1,k of Qj1f at the next coarser scale. 4. Method-I Wavelets and multiscale may act as the natural functional framework for stabilized advection–dispersion problems [20]. In this paper, compactly supported piecewise polynomials B-spline wavelets, which can be generated by convolution as well as recursive relation proposed by Cox and Carl de Boor [22], are used. The constant B-spline over the ith subinterval is defined as
/1k ðxÞ ¼
1 tk 6 x < t kþ1 0
otherwise
;
ð7aÞ
where subscript j denoting multiscale level is dropped for convenience. The mth order B-spline can be recursively calculated by using a formula of Cox- de Boor
/m k ðxÞ ¼
x tk t kþm x m1 /m1 ðxÞ þ / ðxÞ: t kþm t kþ1 kþ1 t kþm1 tk k
ð7bÞ
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3789
With equally spaced ‘‘knots tk’’ of unit length, (7) can be simplified greatly because the B-spline blending function is same for all the segments. The support of the basis function goes up one point for each increase in order m i.e., for any positive integer m, sup[/m] = [0, m]. In the uniform knot vector case, the functions are successive convolutions of the ‘‘box’’ function. The first order spline is defined as a box function /1 = v[0, 1)(x). The mth order spline can be obtained by the convolution
/m ðxÞ ¼ ð/m1 /1 ÞðxÞ ¼
Z
1
/m1 ðx tÞdt;
ð8Þ
0
where ⁄ denotes the convolution operator. Chui and Wang [23] and Unser et al. [24] developed independently B-spline wavelets for multiresolution analysis. The wavelet function for mth order B-spline
wm ¼
2m2 X
1 2
ð1Þi /2m ði þ 1ÞU2m m ð2x iÞ
m1
ð9aÞ
i¼0
with support [0, 2 m 1], is a basic wavelet that generates wavelet function, and
U2m m ð2xÞ ¼
m X
ð1Þi
i¼0
m i
/m ðx iÞ:
ð9bÞ
The 2nd order B-spline and its corresponding wavelet function are shown in Fig. 1. On each level of resolution j an unknown function f(x) can be approximated by using B-spline scaling function
fj ðxÞ ¼
X
aj;k /m ð2j x kÞ ¼
X
k
aj;k /m j;k ðxÞ:
ð10Þ
k
We can define the detail function dj (x) as
dj ðxÞ ¼
X
X
bj;k wm ð2j x kÞ ¼
k
bj;k wm j;k ðxÞ;
ð11Þ
k
where aj,k and bj,k are the scaling and wavelet coefficients at level ‘‘j’’ that defines special resolution of the approximations. The following relationship holds for the function at next higher resolution:
fjþ1 ðxÞ ¼ fj ðxÞ þ dj ðxÞ
ð12aÞ
or
fjþ1 ðxÞ ¼
X
X
ajþ1;k /m jþ1;k ðxÞ ¼
k
aj;k /m j;k ðxÞ þ
k
X
bj;k wm j;k ðxÞ:
ð12bÞ
k
If we start with the level j0 scaling function, then the approximation of level ‘‘j + 1’’ is
fjþ1 ðxÞ ¼
X k
aj0;k /m j0 ;k ðxÞ þ
j X X i¼j0
bi;k wm i;k ðxÞ:
ð13Þ
k
Similarly, concentration function can be expressed as the linear combination of scaling and wavelet functions
Cðx; tÞ ¼
X k
aj0;k ðtÞ/m j0 ;k ðxÞ þ
j X X i¼j0
bi;k ðtÞwm i;k ðxÞ:
k
Fig. 1. (a) Second order B-spline basis function, (b) corresponding wavelet.
ð14Þ
3790
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Substituting this equation in the weak form (4) and using the weight function w as mth order B-spline /m and its wavelet functions wm, we get
_ ¼ 0: ½KfCg þ ½MfCg
ð15Þ
Considering the level of resolution from jo to j, the stiffness and mass matrices of the dynamic Eq. (15) are
2
K/jo /jo
K/jo wjo
6K 6 wjo/jo 6 ½K ¼ 6 6 6 4
K/jo wj 7 7 7 7 7 7 5 Kwj wj
Kwjo wjo
Kwj /jo
3
: K/jo wj
Kwj wjo
ð16aÞ
and
2
M/jo /jo 6M 6 wjo/jo 6 ½M ¼ 6 6 6 4 Mwj /jo
M/jo wjo Mwjo wjo Mwj wjo
M/jo wj
3
M/jo wj 7 7 7 7 7; 7 5 Mw j w j
ð16bÞ
where the typical elements of stiffness and mass matrices at any resolution are
2RL 6 6 K/jo /jo ¼ 6 6 4
0
ð/jo;1
RL
@/jo;1 @x
@/ ð/jo;n @xjo;1 0
þD
þ
@/jo;1 @/jo;1 @x @x
þ k/jo;1 /jo;1 Þdx
@/ @/ D @xjo;n @xjo;1
þ k/jo;1 /jo;1 Þdx : :
RL 0
ð/jo;1
@/jo;n @x
þD
@/jo;1 @/jo;n @x @x
þ k/jo;1 /jo;n Þdx
RL
@/ ð/jo;n @xjo;n 0
þ
@/ @/ D @xjo;n @xjo;n
3 7 7 7 7 5
ð17aÞ
þ k/jo;1 /jo;n Þdx
and
2RL 6 6 M/jo wi ¼ 6 6 4
0
RL 0
/jo;1 wi;1 dx
/jo;n wi;1 dx
RL 0
/jo;1 wi;n dx
RL
00
/jo;n wi;n dx
3 7 7 7: 7 5
ð17bÞ
Similarly, the stiffness and mass matrices at other resolution can be obtained. The unknown coefficients
3 a_ jo 7 6 6b 7 6_ 7 6 jo 7 _ ¼ 6 bjo 7 fCg ¼ 6 7 and fCg 7 6 4 5 5 4 _ bj bj 2
ajo
3
2
ð18Þ
at time ts+1 can be obtained from the coefficients at ts by using the recursive relation of a-family of approximation
f½M þ Dt a½KgfCgsþ1 ¼ f½M ð1 aÞDt½KgfCgs :
ð19Þ
Algorithm Outline The algorithm that uses explicitly the scaling and wavelet basis stands on the fact that the large wavelet coefficient represent high gradient region on the physical domain. The problem can be solved by inverting the M at each scale once and for all at the beginning of the calculation. Then the complete solution can be obtained by matrix multiplication at each time advancing step. On the other hand, our main goal is to compute scaling and wavelet coefficients in a region as close as possible to the space where sharp transition occurs. The time marching version of our algorithm is designed in order to fulfill efficiently this requirement. It involves the following steps.
Fig. 2. (a) Vj, (b) Vj1 Wj1, (c) Vj2 Wj2 Wj1.
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3791
Fig. 3. Solution of the linear advection–dispersion equation using Method I.
Fig. 4. Evolution, decay and marching of wavelet coefficients.
Flowchart for construction of wavelet-Galerkin method. i at ith node. 1. Discretize scaling function /i and its derivative @/ @x i 2. Discretize wavelet wi and its derivative @w at ith node. @x i i 3. Translate periodic /i ; @/ ; wi and @w and use pointwise product and Newton cotes integration to calculate elements of @x @x stiffness matrix K and mass matrix M. 4. Select first few nodes to form ‘‘active’’ space. 5. Loop over the time domain a. Calculate scaling coefficient by solving dynamic Eq. (19) in the active space b. Use wavelet spectrum to estimate active space for the next time step. 6. End loop 7. Use Coiflet filter to eliminate spurious oscillation
3792
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 5. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
Fig. 6. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
5. Method-II In this variational multiscale stabilization method, the problem is discretized at the finest scale and then forward transformation is applied for multiscale projection. The standard Galerkin finite element method employing simple polynomialbased elements is used for ease of calculations of inner products. Similar approach has been used by Brezzi et al. [21], who have accounted fine scale in order to accurately calculate the coarse scale without ‘seeing’ the fine scale feature. Let the finite element equation at the finest resolution J be
AJ C J ¼ fJ :
ð20Þ j
j
Using B-spline wavelets [25] for forward transformation [P jQ ] = TJ, the finest scale is replaced by multiscale as shown in the Fig. 2, i.e.,
T J1
C J1 ¼ CJ : dJ1
ð21Þ
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3793
Fig. 7. Comparison of Method I with analytical solution and Saulyevc scheme solution [26].
Fig. 8. Comparison of Method I with analytical solution.
Substituting (21) in (20), we get
AJ T J1
C J1 ¼ fJ dJ1
ð22Þ
or
T TJ1 AJ T J1
C J1 ¼ T TJ1 fJ ; dJ1
ð23Þ
where CJ1 are the scaling coefficients, dJ1 are the wavelets coefficients at the lower resolution. The wavelet transform can be continued until a desired coarsest level is achieved. The modified stiffness matrix can be expressed as
T TJ1 AJ T J1 ¼ AJ1 :
ð24Þ
The objective of multiscale is to eliminate small detail coefficients and develop the compressed stiffness matrix from the finest scales solution which will be considerably less expensive to solve. It will also eliminate the process of re-meshing the FE domain and re-computing the stiffness matrix with change in location of sudden jump. Multiscale is performed by
3794
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 9. Comparison of Method I, Method II and analytical solution.
Fig. 10. Concentration using Method I before filtering and after fifth level filtering.
projecting the solution at the fine scale space Vj onto next coarser levels Wj1 and Vj1, recursively. It is well established fact that the coefficients of wavelet space are negligible in the smooth region. Retaining grid points associated with wavelets coefficients at high gradient region and eliminating them from a priori known smooth region will generate adaptive grid. In other words, at any time, the computational grid should include points near the sharp gradient zone. Computational grid should also consist of grid points associated with wavelets whose coefficients can possibly become significant during the next time step. Finite element stiffness matrix at the finest scale will be transformed to the coarsest scale only once, i.e. there is no need of re-meshing and re-computing the stiffness matrix with change of sharp gradient in space and time. Compression of stiffness matrix by eliminating the small wavelet coefficients from the smooth region will be required in all subsequent steps of adaptive solution. Fig. 2 depicts wavelet transformation. Coiflet filter is included to eliminate high frequency and small magnitude oscillations which are typical numerical instability phenomena near jump at high Peclet number.
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3795
Fig. 11. Concentration using Method II for coarse grid before filtering and after seventh level filtering.
Fig. 12. Concentration using Method II for fine grid before filtering and after seventh level filtering.
6. Results and discussions To demonstrate the algorithm in our experiment we begin with the simple example of linear advection–dispersion, which can be viewed as a simple case of complex three-dimensional pollution problem. Mostly advection–dispersion solvers use upwind scheme to reduce numerical instability. In the present methods, wavelet based adaptive grid and filtering are used for accuracy and numerical stability. Simple time marching algorithm (a-family) is employed only on the small set of active nodes. The decay constant k = 0 is assumed in this analysis. The change in concentration with time for Peclet number Pe = 50 at locations x = 200 m and 500 m are shown in Fig. 3 and corresponding wavelet coefficients are shown in Fig. 4. It can be observed that changes in wavelet coefficients are very small compared to the scaling coefficients and start well before the sudden jump. This helps in automatic grid refinement with time at the location of sudden jump. The results of Method-I are compared with the analytical method and method based on upwind scheme [26] in Figs. 5–7 for Pe = 10, 100,
3796
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
Fig. 13. Detail coefficients from first to eighth level for Method II.
200, respectively. It can be observed that the solution by Method-I is in good agreement with the analytical solution even at high Peclet number. However, spurious oscillations at very high Peclet number Pe = 500 are clearly visible and shown in Fig. 8. Newton–cotes method is used for calculation of inner product which makes Method-I inefficient and it cannot be extended in two-dimensional advection–dispersion problems. Therefore, it is necessary to use numerical techniques proposed by Beylkin [17] and Dahmen and Micchelli [18] and will be used in future work. Method-II utilizes the stiffness matrix developed by finite element method. The stiffness matrix is projected into wavelet space and obtained multiscale structure. A comparison of Method-II with Method-I and analytical solution is shown in Fig. 9. It can be observed that Method-II is not as stable as Method-I. However, Method-II is more efficient than Method-I. MATLAB programs developed for calculation of stiffness matrix for Method- I takes 4.35 s while Method-II takes only 1.13 s. Fig. 10 shows the response of Method-I at Pe = 1 and the result is improved by using Coiflet filters. It can be observed that the filtering reduces the spurious oscillations, significantly. Similar filtering at Pe = 1 in Method-II is applied and the results before filtering and after filtering are shown in Fig. 11. Adaptive grid refinement improves the results and it is shown in Fig. 12. A typical pattern of wavelet coefficients which are eliminated during filtering in the preceding case at various levels are shown in Fig. 13. The elimination of wavelet coefficient up to level VII improves the result which is presented in Fig. 12. However, deterioration of the result is
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
3797
Fig. 14. Concentration using Method II before filtering and after eighth level filtering.
clearly visible in Fig. 14 due to unjustifiable elimination of wavelet coefficients of level VIII which have large magnitude and low frequency. 7. Conclusion In case of hyperbolic PDEs which show jump discontinuity, the standard numerical schemes such as finite element and finite difference will in general exhibit spurious oscillations. To overcome this difficulty, two multiscale wavelet based methods are proposed and wavelet filtering is used. Results are compared with analytical and other numerical solutions. It has been observed that after filtering, Method-II is as good as Method-I. Generation of stiffness matrix in Method-I is more expensive than in Method-II. Therefore, a customized numerical technique in Method-I should be employed to improve its efficiency. Moreover, the original spur for reduction of degree of freedom, which will be very useful for higher dimensional problems, can be easily achieved in Method-II. References [1] Y. Meyer, Wavelets with compact support, in: U. Chicago: Zygmund Lectures, 1987. [2] S.G. Mallat, A theory for multiresoltion signal decomposition: The wavelet representation, Communications on Pure and Applied Mathematics 41 (1988) 674–693. [3] I. Daubechies, Orthonormal bases of compactly supported wavelets, Communications on Pure and Applied Mathematics 41 (1988) 909–996. [4] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, MA, 1996. [5] J. Liandrat and Ph. Tchamitchian, Resolution of the 1D regularized Burgers equation using a spatial wavelet approximation, Report No: NASA CR – 187480, NASA Langley Research Centre, Hamptonva, (1990). [6] G. Beylkin, J.M. Keiser, On the adaptive numerical solution of nonlinear partial differential equations in wavelet bases, Journal of Computational Physics 132 (1997) 233–259. [7] K. Amaratunga, J.R. Williams, Wavelet based Green’s function approach to 2D PDEs, Engineering Computations 10 (4) (1993) 349–367. [8] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, Journal of Computational Physics 106 (1993) 155–175. [9] W. Swelden, The lifting scheme: A construction of second generation wavelets, SIAM Journal of Mathematical Analysis 29 (2) (1998) 511–546. [10] J. Carnicer, W. Dahmen, J. Pena, Local decompositions of refinable spaces, Applied Computational and Harmonic Analysis 3 (1996) 125–153. [11] W. Dahmen, R. Stevenson, Element-by-element construction of wavelets satisfying stability and moment conditions, SIAM on Numerical Analysis 37 (1999) 319–352. [12] O.V. Vasilyev, S. Paolucci, A dynamically adaptive multilevel wavelet collocation method for solving partial differential equations in finite domain, Journal of Computational Physics 125 (1996) 498–512. [13] P. Krysl, E. Grinspun, P. Schroder, Natural hierarchical refinement for finite element methods, International Journal for Numerical Methods in Engineering 56 (8) (2002) 1109–1124. [14] K. Amaratunga, R. Sudarshan, Multiresolution modeling with operator-customized wavelet derived from finite elements, Computer Methods in Applied Mechanics and Engineering 195 (2006) 2509–2532. [15] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Cambridge University Press, 1987. [16] O.C. Zienkiewicz, R.L. Taylor, Finite Element Method Volume 3: Fluid Dynamics, Fifth ed., Butterworth-Heinemann, 2000. [17] G. Beylkin, On the representation of operators in bases of compactly supported wavelets, SIAM Journal on Numerical Analysis 6 (6) (1992) 1716–1740. [18] W. Dahmen, C.A. Micchelli, Using the refinement equation for evaluating integrals of wavelets, SIAM J. Numerical Analysis 30 (1993) 507–537. [19] A. Cohen, S.M. Kaber, S. Muller, M. Postel, Fully adaptive multiresolution finite volume schemes for conservation laws, Mathematics of Computation 24 (2001) 183–225.
3798 [20] [21] [22] [23]
K. Sandeep et al. / Applied Mathematics and Computation 218 (2011) 3786–3798
S. Bertoluzza, C. Canuto, A. Tabacco, Negative norm R stabilization of convection–diffusion problems, Applied Mathematics Letter 13 (2000) 121–127. F. Brezzi, L.P. Franca, T.J.R. Hughes, A. Russo, b ¼ g, Computer Methods in Applied Mechanics and Engineering 145 (1997) 329–340. D.F. Rogers, J.A. Adams, Mathematical Elements for Computer Graphics, McGraw-Hill Publishing Company, Singapore, 1990. C.K. Chui, J.Z. Wang, On compacty supported spline wavelets and a duality principle, Transaction of the American Mathematical Society 330 (2) (1992) 903–915. [24] M. Unser, A. Aldroubi, M. Eden, A family of polynomial spline wavelet transforms, NCRR rep.153/90, Nat. Inst. Health, (1990). [25] E.J. Stollnitz, T.D. Derose, D.H. Salesin, Wavelets for Computer Graphics, Morgan Kaufmann Pub. Inc., California, 1996. 209–213. [26] Li. Guoyuan, C.R. Jackson, Simple, accurate, and efficient revisions to MacCormack and Saulyev schemes: High Peclet numbers, Applied Mathematics and Computation 186 (2007) 610–622.