Time-accurate solutions of Kortewegâde Vries equation using wavelet ...
Applied Mathematics and Computation 162 (2005) 447–460 www.elsevier.com/locate/amc
Time-accurate solutions of Korteweg–de Vries equation using wavelet Galerkin method B.V. Rathish Kumar *, Mani Mehra Department of Mathematics and Scientific Computing, Indian Institute of Technology, Kanpur 208016, India
Abstract In this study we propose a space and time-accurate numerical method for Korteweg– de Vries equation. In deriving the computational scheme, Taylor generalized Euler time discretization is performed prior to wavelet based Galerkin spatial approximation. This leads to the implicit system which can also be solved by explicit algorithms. Korteweg– de Vries equation is also solved by a operator splitting method using wavelet-Taylor– Galerkin approach. Asymptotic stability of the schemes are verified. 2004 Elsevier Inc. All rights reserved. Keywords: Korteweg–de Vries equation; Wavelets; Taylor–Galerkin method; Operator splitting
1. Introduction The non-linear waves in a shallow water becomes of great interest to many researchers. Its importance comes from its capability for describing many physical phenomena in the fields of physics, mathematics and engineering. These non-linear waves were mathematically modelled by Korteweg–de Vries [1]. The derived equation is well known Korteweg–de Vries (KdV) equation [2]. Such a powerful equation could successfully simulate the red spots at Jupiter. Moreover, the giant ocean waves known as (Tsunami) are also described by the
*
Corresponding author. E-mail addresses: [email protected] (B.V. Rathish Kumar), [email protected] (M. Mehra).
0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.104
448
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
KdV equation. The giant internal waves in the interior of the ocean arising from the temperature difference which may destruct marine vessels could be also described by such a power KdV equation. We consider it in the form ut ¼ Lu þ Nf ðuÞ, where L is linear differential operator and Nf ðuÞ is a nonlinear function. The application of methods based on wavelets to the numerical solution of partial differential equations (PDEs) has recently been studied both from the theoretical and the computational point of view due to its attractive feature: orthogonality, arbitrary regularity, good localization. Wavelet bases seem to combine the advantages of both spectral and finite element basis. Exploration of the usage of Daubechies wavelets to solve PDEs has been undertaken by a number of investigators such as Beylkin et al. [3], Glowinski et al. [4], Latto and Tenenbaum [5], and Qian and Weiss [6]. The fundamental concept behind the Taylor–Galerkin approach is to incorporate more analytical information into the numerical scheme in the most direct and natural way, so that the technique may be regarded as an extension to PDEs of the Obrechkoff methods [7] for ordinary differential equations. As a matter of fact this concept is not new and similar procedures have already been considered in the context of finite difference method [8,9] and has also been considered in conjunction with a spatial representation of spectral type [10]. Later Donea [11,12] has used it in deriving a time-accurate finite element scheme. Primarily their approach consists of extending the Taylor series in the time increment to the second-order before spatial discretization. This procedure has not be implemented so far in the wavelet approach to PDEs. Wavelet basis are known for their spacial accuracy. The aim of the present paper is to formulate wavelet-Taylor–Galerkin method (W-TGM) for KdV and generalized KdV equation in one dimension. The approximate factorization technique [13] when applied to the linear system resulting from applying the numerical scheme will lead to an simple explicit scheme for solving the linear system. This implicit time stepping scheme also leads to circulant system which can be solved using FFTs. We compute the solution by explicit scheme and compare it with the solution obtained from usual matrix inversion. Next, we also apply the operator splitting method in conjunction with appropriate wavelet-Taylor–Galerkin schemes for non-linear conservation and linear dispersive phases for solving KdV equation is proposed. We will now describe the content of the paper more precisely. In Section 2 we summarize some basics of wavelet analysis. In Section 3 we consider the application of W-TGM to KdV equation. In this section we also discuss approximate factorization techniques and operator splitting approach for KdV equation. In Section 4 we have studied stability of the scheme in terms of ODE context. Numerical results are provided in Section 5. Conclusion are given in Section 6.
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
449
2. Wavelet preliminaries Ingrid Daubechies define the class of compactly supported wavelets [14]. P Briefly, let / be a solution of the scaling relation /ðxÞ ¼ k ak /ð2x kÞ. The ak are a collection of coefficients that categorize the specific wavelet basis. The expression / is called the scaling function. The associated wavelet function w is P k defined by the equation wðxÞ ¼ ð1Þ a1k /ð2x kÞ. The normalization k R / dx ¼ 1 of the scaling function leads to the condition X ak ¼ 2: k
The translates of / are required to be orthonormal i.e. Z /ðx kÞ/ðx mÞ dx ¼ dk;m : The scaling relation implies the condition N 1 X
ak ak2m ¼ d0;m ;
k¼0
where N is the order of wavelet. For coefficients verifying the above two conditions, the functions consisting of translates and dilations of the wavelet function, wð2j x kÞ, form a complete, orthogonal basis for square integrable functions on the real line, L2 ðRÞ. If only a finite number of the ak are non-zero then / will have compact support. Since Z X k /ðxÞwðx mÞ dx ¼ ð1Þ a1k ak2m ¼ 0 k
the translates of the scaling function and wavelet define orthogonal subspaces Vj ¼ f2j=2 /ð2j x kÞ; m ¼ . . . ; 1; 0; 1; . . .g; Wj ¼ f2j=2 wð2j x kÞ; m ¼ . . . ; 1; 0; 1; . . .g: The relation Vjþ1 ¼ Vj þ Wj implies the Mallat transform [14] V0 V1 Vjþ1 ; Vjþ1 ¼ V0 W0 W1 Wj :
450
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
Smooth scaling functions arise as a consequence of the degree of approximation of the translates. The result that the polynomials 1; x; . . . ; xp1 be expressed as linear combinations of the translates of /ðx kÞ is implied by the conditions X k ð1Þ k m ak ¼ 0 k
for m ¼ 0; 1; . . . ; p 1. The following are equivalent results. • • • •
. . . ; xp1 g are linear combinations of /ðx R kÞ. f1; x; P kf cjK /ð2j x kÞk 6 C2jp kf p k, where cjk ¼ f ðxÞ/ð2j x kÞ dx. R m ¼ 0 for m ¼ 0; 1; . . . ; p 1. R x wðxÞ dx f ðxÞwð2j xÞ dx 6 c2jp .
For the Daubechies scaling/wavelet function DN have p ¼ N =2. In Fig. 1 we see an example of a compactly supported scaling function and its associated fundamental wavelet function. By rescaling and translation we obtain a complete orthonormal system for L2 ðRÞ which has a sufficient smoothness to also be a basis for H 1 ðRÞ. This wavelet system then yields a basis for second-order elliptic boundary problems on intervals on the real line. The illustrated example has fundamental support ½0; 5 . For arbitrarily large even N there is Daubechies example of a fundamental scaling function defining a wavelet family with support in the interval ½0; N 1 [14]. For any j, l 2 Z we define the 1-periodic scaling function ~ ðxÞ ¼ / j;l
1 X
1 X
/j;l ðx þ nÞ ¼ 2j=2
n¼1
/ð2j ðx þ nÞ lÞ;
x2R
n¼1
2 Scaling function Wavelet function
1.5 1 0.5 0 –0.5 –1 –1.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 1. Daubechies scaling and wavelet functions for N ¼ 6 with support on ½0; 5 .
ð2:1Þ
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
451
and the 1-periodic wavelet ~ ðxÞ ¼ w j;l
1 X
wj;l ðx þ nÞ ¼ 2j=2
n¼1
1 X
wð2j ðx þ nÞ lÞ;
x 2 R:
ð2:2Þ
n¼1
The 1 periodicity can be verified as follows 1 X
~ ðx þ 1Þ ¼ / j;l
/j;l ðx þ n þ 1Þ ¼
n¼1
1 X
~ ðxÞ /j;l ðx þ mÞ ¼ / j;l
m¼1
~ ðx þ 1Þ ¼ w ~ ðxÞ. and similarly w j;l j;l 2.1. The wavelet Galerkin method For a PDE of the form F ðu; ut ; . . . ; ux ; uxx ; . . .Þ ¼ 0:
ð2:3Þ
The wavelet approximation of the form uj ðx; tÞ ¼
2j 1 X
~ ðxÞ; cj;k ðtÞ/ j;k
ð2:4Þ
k;l¼0
where cj;k is unknown coefficient of scaling function expansion. Since we assume periodic boundary conditions there is a periodic wrap around in ðx; yÞ and we let the period scale with the number of terms in the expansion. To determine the coefficient of the expansion (2.4) we substitute (2.4) into Eq. (2.3) and again project the resulting expression onto subspace Vj . The projection requires cj;k to satisfy the equations Z 1 ~ ðxÞF ðuj ; ujt ; ujx ; ujxx ; . . .Þ dx ¼ 0: / j;m 1
To evaluate this expression we must know the coefficients of the form Z 1 Y n d1 ;d2 ;...;dn ^ðl1 ; l2 ; . . . ; ln ; d1 ; d2 ; . . . ; dn Þ ¼ ^l1 ;l2 ;...;ln ¼ /dlii ðyÞ dy ¼
Z
1 i¼1
1
n Y
1 i¼1
/dlii ðyÞ dy:
We can alter a doubly subscripted connection coefficient in to a single subscripted one, and a triply subscripted connection coefficient in to a doubly subscripted one. We therefore define the two and three-term connection coefficient as Z 1 ^dlm1 d2 d3 ¼ /d1 ðyÞ/dl 2 ðyÞ/dm3 ðyÞ dy; 1
452
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
where di P 0. Since the scaling function used to define compact wavelets has a limited number of derivatives, the numerical evaluation of these expressions is often unstable or inaccurate. Special algorithm to evaluate the connection coefficients is devised by Latto et al. [15].
3. Wavelet-Taylor–Galerkin method The first step in the derivation of Taylor–Galerkin schemes [11] is to formulate a high-order time stepping scheme algorithm before the discretization of the spatial variable, which is therefore left momentarily continuous. There are however important differences between the present approach and the aforementioned studies. In the stability analysis the Taylor expansion is used to determine the truncation error of a numerical scheme which has been obtained by discretizing both space and time; here the expansion is used instead to discretize the differential equation in time only. In the first case, starting from a fully discretized equation, a partial differential equation (the modified equation) is generated that the given numerical scheme actually solves; in the present case, starting from the governing partial differential equation, a time stepping algorithm (the generalized time discretized equation) is obtained which is to be employed after adequate wavelet spatial discretization in the actual computations. As an illustration of the procedure above, we are considering KdV and generalized KdV equation. Korteweg–de Vries equation (KdV): We consider it in the form of 3 2 ut ¼ Lu þ Nf ðuÞ where L ¼ a oxo þ b oxo 3 , N ¼ a oxo and f ðuÞ ¼ u2 , then equation becomes ou o u2 o3 ¼a uþ ð3:1Þ þ b 3 u: ot ox 2 ox Time discretization: We have used second-order W-TGM because too many terms are introduced in the third-order time derivative term, especially for nonlinear problems. Let us leave the spatial variable x continuous and discretize (3.1) in time by the following forward Taylor series expansion, n
ðut Þ ¼
unþ1 un Dt n utt OðDt2 Þ 2 Dt
ð3:2Þ
which includes first-order and second-order time derivatives, while the former is provided directly by (3.1), the latter can be obtained by taking the time derivative of the governing PDEs. 2 u ut ¼ aox u þ aox ð3:3Þ þ bo3x u: 2
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
453
The time derivative of (3.3) is utt ¼ aox ut þ aox ðuut Þ þ bo3x ut
ð3:4Þ
and the substitution of (3.3) and (3.4) into the Taylor series expansion (3.2) gives Dt 1 ðaox þ aox un þ aun ox þ bo3x Þ ðunþ1 un Þ=Dt 2 2 n ðu Þ n ¼ aox u þ aox ð3:5Þ þ bo3x un : 2 Spatial discretization: To obtain fully discrete equation we apply wavelet Galerkin method (WGM) to Eq. (3.5) with approximation of the form uj ðx; tÞ ¼
2j 1 X
~ ðxÞ; cj;k ðtÞ/ j;k
ð3:6Þ
k¼0
where cj;k is unknown coefficient of scaling function expansion. The associated wavelet Galerkin equations are ! nþ1 nþ1 X X X cj;l cnj;l Dt cj;k cnj;k j 01 3j=2 j 3 03 a2 ^lk þ a2 ml;k þ bð2 Þ ^lk 2 Dt Dt k k k X XX X 3j=2 j 3 ¼ a2j cnj;k ^01 cnj;k cnj;m ^001 cnj;k ^03 lk þa2 lkm þbð2 Þ lk k
k
m
k
ð3:7Þ for l ¼ 0; 1; . . . ; 2j 1 and ml;k is given by the following expressions involving three level connection coefficient X 010 ml k ¼ cnj;m ð^001 lkm þ ^lkm Þ: m
Therefore calculating the coefficient cj;k reduces to solving a matrix equation: Acu ¼ F ; where cu denote the vector of scaling function coefficient corresponding to u. If simple Gaussian elimination is used to solve this system, then the cost of finding the vector cu is heavy: Oðn3 Þ operations, where n is the order of the matrix. Significantly better performance can be achieved by use of sparse matrix routines. However, because this system is circulant, using FFT, the solution can be performed in Oðn log2 nÞ operations. We can also use approximate factorization techniques as will describe in Section 3.1.
454
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
3.1. Approximate factorization techniques and explicit schemes The W-TGM requires to solve at each time step a system of linear algebraic equations of the form Acu ¼ F :
ð3:8Þ
F ¼ F ðcnu Þ is a known vector and M varies with time and must be recomputed, the exact factorization of the matrix will lead to a significant computational expense. Therefore approximate factorization techniques appear to be quite attractive for these applications. In the present context it is indeed essential to retain the consistent character of the wavelet-Taylor–Galerkin ÔmassÕ matrix A in the approximate factorization procedure. Let us consider the identity A ¼ L þ ðA LÞ;
ð3:9Þ
where L is the diagonal and positive matrix obtained from A by row-sum technique. X Lii ¼ Aii ; Lij ¼ 0; j 6¼ i: ð3:10Þ i
Since L is diagonal and has positive entries, therefore 1
1
A ¼ L2 ðI þ X ÞL2 ;
ð3:11Þ
where 1
1
X ¼ L2 ðA LÞL2 :
ð3:12Þ
Then, under the assumption kX k 6 1, the inverse of A in the form (3.11) can be expressed by the following series, 1
1
A1 ¼ L2 ðI X þ X 2 X 3 þ ÞL2 :
ð3:13Þ
Truncating the series after the term X gives the following two-term approximation of the inverse of A, 1 1 1 1 A ð2Þ ¼ 2L I AL ð3:14Þ 2 while retaining the term X 2 produces the three-term approximation 1 A1 ð3Þ ¼ 3L1 I AL1 þ AL1 AL1 3
ð3:15Þ
and so on. The successive approximations in the above approximate factorization technique can be generated by the following multi-pass algorithm. Consider the sequence of approximate solutions V ðgÞ , g ¼ 0; 1; . . . ; G, defined as
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
455
follows, start from V ð0Þ ¼ 0, then, for g ¼ 0; 1; . . . ; G 1, determine V ðgþ1Þ from V ðgÞ by means of the ÔdiagonalÕ linear system LV ðgþ1Þ ¼ F ðA LÞV ðgÞ :
ð3:16Þ
Finally assume V ¼ V ðGÞ The approximation (3.14) and (3.15) can be obtained by this simple procedure with G ¼ 2 and G ¼ 3, respectively. For G ¼ n we will call it as n-pass explicit scheme. 3.2. Wavelet-Taylor–Galerkin splitting method Here the basic idea is to test the strategy of operator splitting method for KdV equation by solving the non-linear conservation law ut ¼ Nf ðuÞ and the linear dispersive equation ut ¼ Lu sequentially with the appropriate W-TGM. We will call it as W-TGMS scheme. First phase: the non-linear conservation problem corresponds to u1t ¼ Nfu1 :
ð3:17Þ
The temporal discretization is achieved by second-order W-TGM ðunþ1 un1 Þ 1 1 ¼ un1t þ Dtun1t þ OðDt2 Þ; ð3:18Þ Dt 2 where un1 ¼ un . The spatial discretization of this semi-discrete equation is achieved by WGM. Second phase: the linear dispersive problem reads u2t ¼ Lu2 and a second-order temporal discretization is obtained by writing 1 1 DtL ðu2nþ1 un2 Þ=Dt ¼ Lun2 ; 2
ð3:19Þ
ð3:20Þ
and unþ2 ¼ unþ1 . A spatially discrete form of (3.20) is again where un2 ¼ unþ1 1 2 obtained by the WGM. Periodic boundary conditions are incorporated into both the first and second phase numerical schemes. 4. Theoretical stability of the linearized schemes We use the notion of asymptotic stability of a numerical method as it is defined in [16] for a discrete problem of the form dU ¼ LU ; dt where L is assumed to be a diagonalisable matrix.
456
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
Definition. The region of absolute stability of a numerical method is defined for the scalar model problem dU ¼ kU dt to be the set of all kDt such that kU n k is bounded as t ! 1. Finally we say that a numerical method is asymptotically stable for a particular problem if, for sufficiently small Dt > 0, the product of Dt times every eigenvalue of L lies within the region of absolute stability. Most crucial property of L is its spectrum, for this will determine the stability of the time discretization.
5. Numerical results and discussions In this section we present the results of numerical experiment in which we compute approximation to the solution of KdV equation [16] and generalized KdV equation [17]. All the results we present are obtained using Daubechies D6 scaling function and for j ¼ 6. KdV equation: The analytical solution of KdV equation (3.1) is given by " # rffiffiffiffiffiffi 1 Du 2 uðx; tÞ ¼ u0 þ Du sech ðx ct pÞ ; 2pkD 6 1 c ¼ 2p 1 þ u0 þ Du ; 3 "rffiffiffiffiffiffi, # pffiffiffiffiffiffiffiffiffi Du u0 ¼ 2kD 6Du tanh kD 24 on ð1; 1Þ. Chose kD ¼ 0:01, Du ¼ 0:2, a ¼ 2p, b ¼ 12 k2D ð2pÞ3 imposed periodicity on ð0; 2pÞ. The initial condition will agree with the true solution of the problem. Here W-TGM (W-TGM means implicit W-TGM), 3-pass explicit W-TGM and W-TGMS are compared. The relative L1 errors are reported in Table 1. In Fig. 2 solution to the problem using WGM and W-TGM for Dt ¼ 104 and j ¼ 6 are presented. For stability analysis we will consider linearized KdV equation 3 ou ¼ a oxo u þ b oxo 3 u þ acux where the linearization coefficient c stands for the ot value of u. We assume throughout that jcj 6 1. The kDt for this scheme is plotted in Fig. 3. Generalized KdV equation: In order to test our method for generalized KdV equation we use exact one-soliton solution for the KdV equation 1 ut þ ðu2 Þx þ uxxx ¼ 0: 2
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
457
Table 1 j
WGM
W-TGM
W-TGM (explicit)
W-TGMS
2
10
4 5 6
0.0310 0.0231 0.0179
0.0382 0.0271 0.0196
0.0374 0.0272 0.0196
0.0874 0.0647 0.0423
104
4 5 6
0.0304 0.0184 0.0175
0.0304 0.0184 0.0170
0.0304 0.0184 0.0170
0.0916 0.0656 0.0431
Dt
t=.01
0.2
t=.01
0.2
(b)
(a)
0.1
0.1
u
u
0.15
0.15
0.05
0.05
0
0
–0.05
0
0.2 0.4 0.6 0.8 1
1.2 1.4 1.6 1.8 2
–0.05
0
0.2 0.4 0.6 0.8
X\PI
1 1.2 1.4 1.6 1.8
2
X\PI
Fig. 2. (a) Solution of KdV equation by WGM, (b) by W-TGM.
0.025 0.8
(a)
0.02
0.6
(b)
0.015
0.4
0.01
0.2
0.005
0
0
–0.2
–0.005
–0.4
–0.01
–0.6
–0.015
–0.8
–0.02 –1.5
–1
–0.5
0
0.5
–0.025 –0.03 –0.02 –0.01
0
0.01
0.02
Fig. 3. For W-TGM Dt times the eigenvalues of L6 : (a) Dt ¼ 102 , (b) Dt ¼ 104 .
0.03
458
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460 t=.01
1 0.9 0.8 0.7
u
0.6 0.5 0.4 0.3 0.2 0.1 0 –2
–1.5
–1
–0.5
0
0.5
1
1.5
2
X Fig. 4. Solution of KdV equation by W-TGM.
0.5
x 10–3 5
(a)
0.4
4
0.3
3
0.2
2
0.1
1
0
0
–0.1
–1
–0.2
–2
–0.3
–3
–0.4
–4
–0.5 –16 –14 –12 –10 –8
–6
–4
–2
0
2
(b)
–5 –18 –16 –14 –12 –10 –8 –6 –4 –2
–3
x 10
Fig. 5. For W-TGM Dt times the eigenvalues of L6 : (a) Dt ¼ 102 , (b) Dt ¼ 104 .
The one-soliton solution is given by uðx; tÞ ¼ 3c sech2
rffiffiffiffiffi c ðx ctÞ : 4
0
2
x 10–7
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
459
We used ¼ 0:0013020833 and c ¼ 1=3. In this computation we use periodic boundary data in the interval ð2; 2Þ. For stability analysis we will consider linearized KdV equation ou þ cux þ uxxx ¼ 0 where the linearization coefficient ot c stands for the value of u. We assume throughout that jcj 6 1. The solution is plotted in Fig. 4 for Dt ¼ 104 and j ¼ 6. The kDt for this scheme is plotted in Fig. 5 for Dt ¼ 102 and Dt ¼ 104 . For Dt ¼ 102 time step our scheme is unstable.
6. Conclusion In wavelet-Taylor–Galerkin method the precedence of time discretization to space discretization in conjunction with wavelet bases for expressing spatial terms renders robustness to the proposed schemes and makes them space and time-accurate. We have also used operator splitting in conjunction with wavelet-Taylor–Galerkin method. It is easy to implement on a computer, and one can combine a variety of methods for each of the equations. Numerical results are presented for the solutions of KdV equations.
References [1] D.J. Korteweg, G. de Vries, On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves, Philos. Mag. 5 (39) (1895) 422–443. [2] R.M. Miura, The Korteweg–de Vries equation: a survey of results, SIAM Rev. 18 (1976) 412– 459. [3] G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transforms and numerical algorithms, Commun. Pure Appl. Math. XLIV (1991) 141–183. [4] R. Glowinski, W. Lawton, M. Ravachol, E. Tenenbaum, Wavelet solutions of linear and nonlinear elliptic, parabolic and hyperbolic problems in 1D, Comput. Methods Appl. Sci. Eng., SIAM, PA 4 (1990) 55–120. [5] A. Latto, E. Tenenbaum, Les ondelettes a support compact et la solution de 1Õequation de Burgers, C.R. Acad. Sci. Paris, Ser. I (1990) 903–909. [6] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, J. Comput. Phys. 106 (1992) 155–175. [7] J.D. Lambert, Computational Methods for Ordinary Differential Equations, Wiley, London, 1973. [8] P.D. Lax, B. Wendroff, Systems of conservation laws, Commun. Pure Appl. Math. 13 (1960) 217. [9] P.D. Lax, B. Wendroff, On the stability of difference schemes, Commun. Pure Appl. Math. 15 (1962) 363. [10] J. Gazdag, Numerical convective schemes based on accurate computation of space derivatives, J. Comput. Phys. 13 (1973) 100–113. [11] J. Donea, A Taylor–Galerkin method for convective transport problems, Int. J. Numer. Methods Eng. 20 (1984) 101–119. [12] J. Donea, L. Quartapelle, V. Selmin, An analysis of time discretization in finite element solution of hyperbolic problems, J. Comput. Phys. 70 (1987) 463.
460
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
[13] J. Donea, S. Giuliani, H. Laval, Time-accurate solution of advection-diffusion problems by finite elements, Comput. Methods Appl. Mech. Eng. 45 (1984) 123–146. [14] I. Daubechies, Orthonormal basis of compactly supported wavelets, Commun. Pure Appl. Math. 41 (1988) 906–966. [15] A. Latto, H.L. Resnikoff, E. Tenenbaum, The evaluation of connection coefficients of compactly supported wavelets, in: Proceedings of the French–USA Workshop on Wavelets and Turbulence, Princeton, Springer-Verlag, New York, 1991. [16] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, 1988. [17] H. Holden, K.H. Karlsen, N.H. Risebro, Operator splitting methods for generalized KdV equations, J. Comput. Phys. 153 (1999) 203–222.
Time-accurate solutions of Korteweg–de Vries equation using wavelet Galerkin method B.V. Rathish Kumar *, Mani Mehra Department of Mathematics and Scientific Computing, Indian Institute of Technology, Kanpur 208016, India
Abstract In this study we propose a space and time-accurate numerical method for Korteweg– de Vries equation. In deriving the computational scheme, Taylor generalized Euler time discretization is performed prior to wavelet based Galerkin spatial approximation. This leads to the implicit system which can also be solved by explicit algorithms. Korteweg– de Vries equation is also solved by a operator splitting method using wavelet-Taylor– Galerkin approach. Asymptotic stability of the schemes are verified. 2004 Elsevier Inc. All rights reserved. Keywords: Korteweg–de Vries equation; Wavelets; Taylor–Galerkin method; Operator splitting
1. Introduction The non-linear waves in a shallow water becomes of great interest to many researchers. Its importance comes from its capability for describing many physical phenomena in the fields of physics, mathematics and engineering. These non-linear waves were mathematically modelled by Korteweg–de Vries [1]. The derived equation is well known Korteweg–de Vries (KdV) equation [2]. Such a powerful equation could successfully simulate the red spots at Jupiter. Moreover, the giant ocean waves known as (Tsunami) are also described by the
*
Corresponding author. E-mail addresses: [email protected] (B.V. Rathish Kumar), [email protected] (M. Mehra).
0096-3003/$ - see front matter 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2003.12.104
448
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
KdV equation. The giant internal waves in the interior of the ocean arising from the temperature difference which may destruct marine vessels could be also described by such a power KdV equation. We consider it in the form ut ¼ Lu þ Nf ðuÞ, where L is linear differential operator and Nf ðuÞ is a nonlinear function. The application of methods based on wavelets to the numerical solution of partial differential equations (PDEs) has recently been studied both from the theoretical and the computational point of view due to its attractive feature: orthogonality, arbitrary regularity, good localization. Wavelet bases seem to combine the advantages of both spectral and finite element basis. Exploration of the usage of Daubechies wavelets to solve PDEs has been undertaken by a number of investigators such as Beylkin et al. [3], Glowinski et al. [4], Latto and Tenenbaum [5], and Qian and Weiss [6]. The fundamental concept behind the Taylor–Galerkin approach is to incorporate more analytical information into the numerical scheme in the most direct and natural way, so that the technique may be regarded as an extension to PDEs of the Obrechkoff methods [7] for ordinary differential equations. As a matter of fact this concept is not new and similar procedures have already been considered in the context of finite difference method [8,9] and has also been considered in conjunction with a spatial representation of spectral type [10]. Later Donea [11,12] has used it in deriving a time-accurate finite element scheme. Primarily their approach consists of extending the Taylor series in the time increment to the second-order before spatial discretization. This procedure has not be implemented so far in the wavelet approach to PDEs. Wavelet basis are known for their spacial accuracy. The aim of the present paper is to formulate wavelet-Taylor–Galerkin method (W-TGM) for KdV and generalized KdV equation in one dimension. The approximate factorization technique [13] when applied to the linear system resulting from applying the numerical scheme will lead to an simple explicit scheme for solving the linear system. This implicit time stepping scheme also leads to circulant system which can be solved using FFTs. We compute the solution by explicit scheme and compare it with the solution obtained from usual matrix inversion. Next, we also apply the operator splitting method in conjunction with appropriate wavelet-Taylor–Galerkin schemes for non-linear conservation and linear dispersive phases for solving KdV equation is proposed. We will now describe the content of the paper more precisely. In Section 2 we summarize some basics of wavelet analysis. In Section 3 we consider the application of W-TGM to KdV equation. In this section we also discuss approximate factorization techniques and operator splitting approach for KdV equation. In Section 4 we have studied stability of the scheme in terms of ODE context. Numerical results are provided in Section 5. Conclusion are given in Section 6.
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
449
2. Wavelet preliminaries Ingrid Daubechies define the class of compactly supported wavelets [14]. P Briefly, let / be a solution of the scaling relation /ðxÞ ¼ k ak /ð2x kÞ. The ak are a collection of coefficients that categorize the specific wavelet basis. The expression / is called the scaling function. The associated wavelet function w is P k defined by the equation wðxÞ ¼ ð1Þ a1k /ð2x kÞ. The normalization k R / dx ¼ 1 of the scaling function leads to the condition X ak ¼ 2: k
The translates of / are required to be orthonormal i.e. Z /ðx kÞ/ðx mÞ dx ¼ dk;m : The scaling relation implies the condition N 1 X
ak ak2m ¼ d0;m ;
k¼0
where N is the order of wavelet. For coefficients verifying the above two conditions, the functions consisting of translates and dilations of the wavelet function, wð2j x kÞ, form a complete, orthogonal basis for square integrable functions on the real line, L2 ðRÞ. If only a finite number of the ak are non-zero then / will have compact support. Since Z X k /ðxÞwðx mÞ dx ¼ ð1Þ a1k ak2m ¼ 0 k
the translates of the scaling function and wavelet define orthogonal subspaces Vj ¼ f2j=2 /ð2j x kÞ; m ¼ . . . ; 1; 0; 1; . . .g; Wj ¼ f2j=2 wð2j x kÞ; m ¼ . . . ; 1; 0; 1; . . .g: The relation Vjþ1 ¼ Vj þ Wj implies the Mallat transform [14] V0 V1 Vjþ1 ; Vjþ1 ¼ V0 W0 W1 Wj :
450
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
Smooth scaling functions arise as a consequence of the degree of approximation of the translates. The result that the polynomials 1; x; . . . ; xp1 be expressed as linear combinations of the translates of /ðx kÞ is implied by the conditions X k ð1Þ k m ak ¼ 0 k
for m ¼ 0; 1; . . . ; p 1. The following are equivalent results. • • • •
. . . ; xp1 g are linear combinations of /ðx R kÞ. f1; x; P kf cjK /ð2j x kÞk 6 C2jp kf p k, where cjk ¼ f ðxÞ/ð2j x kÞ dx. R m ¼ 0 for m ¼ 0; 1; . . . ; p 1. R x wðxÞ dx f ðxÞwð2j xÞ dx 6 c2jp .
For the Daubechies scaling/wavelet function DN have p ¼ N =2. In Fig. 1 we see an example of a compactly supported scaling function and its associated fundamental wavelet function. By rescaling and translation we obtain a complete orthonormal system for L2 ðRÞ which has a sufficient smoothness to also be a basis for H 1 ðRÞ. This wavelet system then yields a basis for second-order elliptic boundary problems on intervals on the real line. The illustrated example has fundamental support ½0; 5 . For arbitrarily large even N there is Daubechies example of a fundamental scaling function defining a wavelet family with support in the interval ½0; N 1 [14]. For any j, l 2 Z we define the 1-periodic scaling function ~ ðxÞ ¼ / j;l
1 X
1 X
/j;l ðx þ nÞ ¼ 2j=2
n¼1
/ð2j ðx þ nÞ lÞ;
x2R
n¼1
2 Scaling function Wavelet function
1.5 1 0.5 0 –0.5 –1 –1.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 1. Daubechies scaling and wavelet functions for N ¼ 6 with support on ½0; 5 .
ð2:1Þ
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
451
and the 1-periodic wavelet ~ ðxÞ ¼ w j;l
1 X
wj;l ðx þ nÞ ¼ 2j=2
n¼1
1 X
wð2j ðx þ nÞ lÞ;
x 2 R:
ð2:2Þ
n¼1
The 1 periodicity can be verified as follows 1 X
~ ðx þ 1Þ ¼ / j;l
/j;l ðx þ n þ 1Þ ¼
n¼1
1 X
~ ðxÞ /j;l ðx þ mÞ ¼ / j;l
m¼1
~ ðx þ 1Þ ¼ w ~ ðxÞ. and similarly w j;l j;l 2.1. The wavelet Galerkin method For a PDE of the form F ðu; ut ; . . . ; ux ; uxx ; . . .Þ ¼ 0:
ð2:3Þ
The wavelet approximation of the form uj ðx; tÞ ¼
2j 1 X
~ ðxÞ; cj;k ðtÞ/ j;k
ð2:4Þ
k;l¼0
where cj;k is unknown coefficient of scaling function expansion. Since we assume periodic boundary conditions there is a periodic wrap around in ðx; yÞ and we let the period scale with the number of terms in the expansion. To determine the coefficient of the expansion (2.4) we substitute (2.4) into Eq. (2.3) and again project the resulting expression onto subspace Vj . The projection requires cj;k to satisfy the equations Z 1 ~ ðxÞF ðuj ; ujt ; ujx ; ujxx ; . . .Þ dx ¼ 0: / j;m 1
To evaluate this expression we must know the coefficients of the form Z 1 Y n d1 ;d2 ;...;dn ^ðl1 ; l2 ; . . . ; ln ; d1 ; d2 ; . . . ; dn Þ ¼ ^l1 ;l2 ;...;ln ¼ /dlii ðyÞ dy ¼
Z
1 i¼1
1
n Y
1 i¼1
/dlii ðyÞ dy:
We can alter a doubly subscripted connection coefficient in to a single subscripted one, and a triply subscripted connection coefficient in to a doubly subscripted one. We therefore define the two and three-term connection coefficient as Z 1 ^dlm1 d2 d3 ¼ /d1 ðyÞ/dl 2 ðyÞ/dm3 ðyÞ dy; 1
452
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
where di P 0. Since the scaling function used to define compact wavelets has a limited number of derivatives, the numerical evaluation of these expressions is often unstable or inaccurate. Special algorithm to evaluate the connection coefficients is devised by Latto et al. [15].
3. Wavelet-Taylor–Galerkin method The first step in the derivation of Taylor–Galerkin schemes [11] is to formulate a high-order time stepping scheme algorithm before the discretization of the spatial variable, which is therefore left momentarily continuous. There are however important differences between the present approach and the aforementioned studies. In the stability analysis the Taylor expansion is used to determine the truncation error of a numerical scheme which has been obtained by discretizing both space and time; here the expansion is used instead to discretize the differential equation in time only. In the first case, starting from a fully discretized equation, a partial differential equation (the modified equation) is generated that the given numerical scheme actually solves; in the present case, starting from the governing partial differential equation, a time stepping algorithm (the generalized time discretized equation) is obtained which is to be employed after adequate wavelet spatial discretization in the actual computations. As an illustration of the procedure above, we are considering KdV and generalized KdV equation. Korteweg–de Vries equation (KdV): We consider it in the form of 3 2 ut ¼ Lu þ Nf ðuÞ where L ¼ a oxo þ b oxo 3 , N ¼ a oxo and f ðuÞ ¼ u2 , then equation becomes ou o u2 o3 ¼a uþ ð3:1Þ þ b 3 u: ot ox 2 ox Time discretization: We have used second-order W-TGM because too many terms are introduced in the third-order time derivative term, especially for nonlinear problems. Let us leave the spatial variable x continuous and discretize (3.1) in time by the following forward Taylor series expansion, n
ðut Þ ¼
unþ1 un Dt n utt OðDt2 Þ 2 Dt
ð3:2Þ
which includes first-order and second-order time derivatives, while the former is provided directly by (3.1), the latter can be obtained by taking the time derivative of the governing PDEs. 2 u ut ¼ aox u þ aox ð3:3Þ þ bo3x u: 2
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
453
The time derivative of (3.3) is utt ¼ aox ut þ aox ðuut Þ þ bo3x ut
ð3:4Þ
and the substitution of (3.3) and (3.4) into the Taylor series expansion (3.2) gives Dt 1 ðaox þ aox un þ aun ox þ bo3x Þ ðunþ1 un Þ=Dt 2 2 n ðu Þ n ¼ aox u þ aox ð3:5Þ þ bo3x un : 2 Spatial discretization: To obtain fully discrete equation we apply wavelet Galerkin method (WGM) to Eq. (3.5) with approximation of the form uj ðx; tÞ ¼
2j 1 X
~ ðxÞ; cj;k ðtÞ/ j;k
ð3:6Þ
k¼0
where cj;k is unknown coefficient of scaling function expansion. The associated wavelet Galerkin equations are ! nþ1 nþ1 X X X cj;l cnj;l Dt cj;k cnj;k j 01 3j=2 j 3 03 a2 ^lk þ a2 ml;k þ bð2 Þ ^lk 2 Dt Dt k k k X XX X 3j=2 j 3 ¼ a2j cnj;k ^01 cnj;k cnj;m ^001 cnj;k ^03 lk þa2 lkm þbð2 Þ lk k
k
m
k
ð3:7Þ for l ¼ 0; 1; . . . ; 2j 1 and ml;k is given by the following expressions involving three level connection coefficient X 010 ml k ¼ cnj;m ð^001 lkm þ ^lkm Þ: m
Therefore calculating the coefficient cj;k reduces to solving a matrix equation: Acu ¼ F ; where cu denote the vector of scaling function coefficient corresponding to u. If simple Gaussian elimination is used to solve this system, then the cost of finding the vector cu is heavy: Oðn3 Þ operations, where n is the order of the matrix. Significantly better performance can be achieved by use of sparse matrix routines. However, because this system is circulant, using FFT, the solution can be performed in Oðn log2 nÞ operations. We can also use approximate factorization techniques as will describe in Section 3.1.
454
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
3.1. Approximate factorization techniques and explicit schemes The W-TGM requires to solve at each time step a system of linear algebraic equations of the form Acu ¼ F :
ð3:8Þ
F ¼ F ðcnu Þ is a known vector and M varies with time and must be recomputed, the exact factorization of the matrix will lead to a significant computational expense. Therefore approximate factorization techniques appear to be quite attractive for these applications. In the present context it is indeed essential to retain the consistent character of the wavelet-Taylor–Galerkin ÔmassÕ matrix A in the approximate factorization procedure. Let us consider the identity A ¼ L þ ðA LÞ;
ð3:9Þ
where L is the diagonal and positive matrix obtained from A by row-sum technique. X Lii ¼ Aii ; Lij ¼ 0; j 6¼ i: ð3:10Þ i
Since L is diagonal and has positive entries, therefore 1
1
A ¼ L2 ðI þ X ÞL2 ;
ð3:11Þ
where 1
1
X ¼ L2 ðA LÞL2 :
ð3:12Þ
Then, under the assumption kX k 6 1, the inverse of A in the form (3.11) can be expressed by the following series, 1
1
A1 ¼ L2 ðI X þ X 2 X 3 þ ÞL2 :
ð3:13Þ
Truncating the series after the term X gives the following two-term approximation of the inverse of A, 1 1 1 1 A ð2Þ ¼ 2L I AL ð3:14Þ 2 while retaining the term X 2 produces the three-term approximation 1 A1 ð3Þ ¼ 3L1 I AL1 þ AL1 AL1 3
ð3:15Þ
and so on. The successive approximations in the above approximate factorization technique can be generated by the following multi-pass algorithm. Consider the sequence of approximate solutions V ðgÞ , g ¼ 0; 1; . . . ; G, defined as
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
455
follows, start from V ð0Þ ¼ 0, then, for g ¼ 0; 1; . . . ; G 1, determine V ðgþ1Þ from V ðgÞ by means of the ÔdiagonalÕ linear system LV ðgþ1Þ ¼ F ðA LÞV ðgÞ :
ð3:16Þ
Finally assume V ¼ V ðGÞ The approximation (3.14) and (3.15) can be obtained by this simple procedure with G ¼ 2 and G ¼ 3, respectively. For G ¼ n we will call it as n-pass explicit scheme. 3.2. Wavelet-Taylor–Galerkin splitting method Here the basic idea is to test the strategy of operator splitting method for KdV equation by solving the non-linear conservation law ut ¼ Nf ðuÞ and the linear dispersive equation ut ¼ Lu sequentially with the appropriate W-TGM. We will call it as W-TGMS scheme. First phase: the non-linear conservation problem corresponds to u1t ¼ Nfu1 :
ð3:17Þ
The temporal discretization is achieved by second-order W-TGM ðunþ1 un1 Þ 1 1 ¼ un1t þ Dtun1t þ OðDt2 Þ; ð3:18Þ Dt 2 where un1 ¼ un . The spatial discretization of this semi-discrete equation is achieved by WGM. Second phase: the linear dispersive problem reads u2t ¼ Lu2 and a second-order temporal discretization is obtained by writing 1 1 DtL ðu2nþ1 un2 Þ=Dt ¼ Lun2 ; 2
ð3:19Þ
ð3:20Þ
and unþ2 ¼ unþ1 . A spatially discrete form of (3.20) is again where un2 ¼ unþ1 1 2 obtained by the WGM. Periodic boundary conditions are incorporated into both the first and second phase numerical schemes. 4. Theoretical stability of the linearized schemes We use the notion of asymptotic stability of a numerical method as it is defined in [16] for a discrete problem of the form dU ¼ LU ; dt where L is assumed to be a diagonalisable matrix.
456
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
Definition. The region of absolute stability of a numerical method is defined for the scalar model problem dU ¼ kU dt to be the set of all kDt such that kU n k is bounded as t ! 1. Finally we say that a numerical method is asymptotically stable for a particular problem if, for sufficiently small Dt > 0, the product of Dt times every eigenvalue of L lies within the region of absolute stability. Most crucial property of L is its spectrum, for this will determine the stability of the time discretization.
5. Numerical results and discussions In this section we present the results of numerical experiment in which we compute approximation to the solution of KdV equation [16] and generalized KdV equation [17]. All the results we present are obtained using Daubechies D6 scaling function and for j ¼ 6. KdV equation: The analytical solution of KdV equation (3.1) is given by " # rffiffiffiffiffiffi 1 Du 2 uðx; tÞ ¼ u0 þ Du sech ðx ct pÞ ; 2pkD 6 1 c ¼ 2p 1 þ u0 þ Du ; 3 "rffiffiffiffiffiffi, # pffiffiffiffiffiffiffiffiffi Du u0 ¼ 2kD 6Du tanh kD 24 on ð1; 1Þ. Chose kD ¼ 0:01, Du ¼ 0:2, a ¼ 2p, b ¼ 12 k2D ð2pÞ3 imposed periodicity on ð0; 2pÞ. The initial condition will agree with the true solution of the problem. Here W-TGM (W-TGM means implicit W-TGM), 3-pass explicit W-TGM and W-TGMS are compared. The relative L1 errors are reported in Table 1. In Fig. 2 solution to the problem using WGM and W-TGM for Dt ¼ 104 and j ¼ 6 are presented. For stability analysis we will consider linearized KdV equation 3 ou ¼ a oxo u þ b oxo 3 u þ acux where the linearization coefficient c stands for the ot value of u. We assume throughout that jcj 6 1. The kDt for this scheme is plotted in Fig. 3. Generalized KdV equation: In order to test our method for generalized KdV equation we use exact one-soliton solution for the KdV equation 1 ut þ ðu2 Þx þ uxxx ¼ 0: 2
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
457
Table 1 j
WGM
W-TGM
W-TGM (explicit)
W-TGMS
2
10
4 5 6
0.0310 0.0231 0.0179
0.0382 0.0271 0.0196
0.0374 0.0272 0.0196
0.0874 0.0647 0.0423
104
4 5 6
0.0304 0.0184 0.0175
0.0304 0.0184 0.0170
0.0304 0.0184 0.0170
0.0916 0.0656 0.0431
Dt
t=.01
0.2
t=.01
0.2
(b)
(a)
0.1
0.1
u
u
0.15
0.15
0.05
0.05
0
0
–0.05
0
0.2 0.4 0.6 0.8 1
1.2 1.4 1.6 1.8 2
–0.05
0
0.2 0.4 0.6 0.8
X\PI
1 1.2 1.4 1.6 1.8
2
X\PI
Fig. 2. (a) Solution of KdV equation by WGM, (b) by W-TGM.
0.025 0.8
(a)
0.02
0.6
(b)
0.015
0.4
0.01
0.2
0.005
0
0
–0.2
–0.005
–0.4
–0.01
–0.6
–0.015
–0.8
–0.02 –1.5
–1
–0.5
0
0.5
–0.025 –0.03 –0.02 –0.01
0
0.01
0.02
Fig. 3. For W-TGM Dt times the eigenvalues of L6 : (a) Dt ¼ 102 , (b) Dt ¼ 104 .
0.03
458
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460 t=.01
1 0.9 0.8 0.7
u
0.6 0.5 0.4 0.3 0.2 0.1 0 –2
–1.5
–1
–0.5
0
0.5
1
1.5
2
X Fig. 4. Solution of KdV equation by W-TGM.
0.5
x 10–3 5
(a)
0.4
4
0.3
3
0.2
2
0.1
1
0
0
–0.1
–1
–0.2
–2
–0.3
–3
–0.4
–4
–0.5 –16 –14 –12 –10 –8
–6
–4
–2
0
2
(b)
–5 –18 –16 –14 –12 –10 –8 –6 –4 –2
–3
x 10
Fig. 5. For W-TGM Dt times the eigenvalues of L6 : (a) Dt ¼ 102 , (b) Dt ¼ 104 .
The one-soliton solution is given by uðx; tÞ ¼ 3c sech2
rffiffiffiffiffi c ðx ctÞ : 4
0
2
x 10–7
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
459
We used ¼ 0:0013020833 and c ¼ 1=3. In this computation we use periodic boundary data in the interval ð2; 2Þ. For stability analysis we will consider linearized KdV equation ou þ cux þ uxxx ¼ 0 where the linearization coefficient ot c stands for the value of u. We assume throughout that jcj 6 1. The solution is plotted in Fig. 4 for Dt ¼ 104 and j ¼ 6. The kDt for this scheme is plotted in Fig. 5 for Dt ¼ 102 and Dt ¼ 104 . For Dt ¼ 102 time step our scheme is unstable.
6. Conclusion In wavelet-Taylor–Galerkin method the precedence of time discretization to space discretization in conjunction with wavelet bases for expressing spatial terms renders robustness to the proposed schemes and makes them space and time-accurate. We have also used operator splitting in conjunction with wavelet-Taylor–Galerkin method. It is easy to implement on a computer, and one can combine a variety of methods for each of the equations. Numerical results are presented for the solutions of KdV equations.
References [1] D.J. Korteweg, G. de Vries, On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves, Philos. Mag. 5 (39) (1895) 422–443. [2] R.M. Miura, The Korteweg–de Vries equation: a survey of results, SIAM Rev. 18 (1976) 412– 459. [3] G. Beylkin, R. Coifman, V. Rokhlin, Fast wavelet transforms and numerical algorithms, Commun. Pure Appl. Math. XLIV (1991) 141–183. [4] R. Glowinski, W. Lawton, M. Ravachol, E. Tenenbaum, Wavelet solutions of linear and nonlinear elliptic, parabolic and hyperbolic problems in 1D, Comput. Methods Appl. Sci. Eng., SIAM, PA 4 (1990) 55–120. [5] A. Latto, E. Tenenbaum, Les ondelettes a support compact et la solution de 1Õequation de Burgers, C.R. Acad. Sci. Paris, Ser. I (1990) 903–909. [6] S. Qian, J. Weiss, Wavelets and the numerical solution of partial differential equations, J. Comput. Phys. 106 (1992) 155–175. [7] J.D. Lambert, Computational Methods for Ordinary Differential Equations, Wiley, London, 1973. [8] P.D. Lax, B. Wendroff, Systems of conservation laws, Commun. Pure Appl. Math. 13 (1960) 217. [9] P.D. Lax, B. Wendroff, On the stability of difference schemes, Commun. Pure Appl. Math. 15 (1962) 363. [10] J. Gazdag, Numerical convective schemes based on accurate computation of space derivatives, J. Comput. Phys. 13 (1973) 100–113. [11] J. Donea, A Taylor–Galerkin method for convective transport problems, Int. J. Numer. Methods Eng. 20 (1984) 101–119. [12] J. Donea, L. Quartapelle, V. Selmin, An analysis of time discretization in finite element solution of hyperbolic problems, J. Comput. Phys. 70 (1987) 463.
460
B.V. Rathish Kumar, M. Mehra / Appl. Math. Comput. 162 (2005) 447–460
[13] J. Donea, S. Giuliani, H. Laval, Time-accurate solution of advection-diffusion problems by finite elements, Comput. Methods Appl. Mech. Eng. 45 (1984) 123–146. [14] I. Daubechies, Orthonormal basis of compactly supported wavelets, Commun. Pure Appl. Math. 41 (1988) 906–966. [15] A. Latto, H.L. Resnikoff, E. Tenenbaum, The evaluation of connection coefficients of compactly supported wavelets, in: Proceedings of the French–USA Workshop on Wavelets and Turbulence, Princeton, Springer-Verlag, New York, 1991. [16] C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, 1988. [17] H. Holden, K.H. Karlsen, N.H. Risebro, Operator splitting methods for generalized KdV equations, J. Comput. Phys. 153 (1999) 203–222.
