Fast Wavelet-Taylor Method for PDEs
Fast Wavelet-Taylor Method for PDEs
We introduce the concept of fast wavelet-Taylor Galerkin methods for the numerical solution of partial
differential equations. In wavelet-Taylor Galerkin method discretization in time is performed before the
wavelet based spatial approximation by introducing accurate generalizations of the standard Euler, and
leap-frog time-stepping scheme with the help of Taylor series expansions in the time step. We will present
two different time-accurate wavelet schemes to solve the PDEs. First, numerical schemes taking advantage
of the wavelet bases capabilities to compress the operators and sparse representation of functions which are
smooth, except for in localized regions, up to any given accuracy are presented. Here numerical
experiments deal with advection equation with the spiky solution in one dimension, two dimensions, and
nonlinear equation with a shock in solution in two dimensions. Second, our schemes deal with more regular
class of problems where wavelets are not efficient procedure for data compression but we can use the good
approximation properties of wavelet. Here time-accurate schemes lead to consistent mass matrix in an
explicit time stepping, which can be solved by approximate factorization techniques. Numerical experi-
ment deals with more regular class of problems like heat equation as well as coupled linear system in two
dimensions. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 22: 274 –295, 2006
Keywords: Taylor-Galerkin method; wavelets; time marching scheme; hyperbolic equations; parabolic
equations
1. INTRODUCTION
Wavelets method is a new numerical tool for solving partial differential equations (PDEs).
Wavelet analysis assumed significance due to successful applications in signal and image
processing during the eighties. Compactly supported wavelets which are differentiable were
introduced by Daubechies in her celebrated article [1], which has had applications in a number
of areas. Since these functions combine orthogonality with localization and scaling properties,
it has been a natural idea to attempt to use these functions for the numerical approximation to
solutions to partial differential equations (PDEs). There have been a number of articles in this
Correspondence to: B. V. Rathish Kumar, Department of Mathematics, Indian Institute of Technology, Kampur 208016,
UP India (e-mail: bvrk@[Link])
© 2005 Wiley Periodicals, Inc.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 275
direction including studies of Burgers’ equation [2, 3] as well as two-dimensional flow [4 – 6].
Glowinski et al. [7] considered wavelet-based variational methods to solve one-dimension linear
and nonlinear ordinary differential equations. Whereas multiresolution analysis and their asso-
ciated scale function bases are used as alternative bases in Galerkin methods, such methods have
thus convergence properties similar to the ones of spectral methods, and simultaneously partial
derivative operators discretize similarly as in finite-difference methods.
In the literature, many tentatives have been performed, often based on Galerkin or Petrov-
Galerkin methods, which use the compression properties of wavelet bases, which contain
specific wavelet methods for PDEs. Some of them take advantage of the wavelet compression
of the solution [2]; others use instead the wavelet compression of the operator [8]. The aim of
the present article is to introduce the wavelet-Taylor Galerkin method, which has the benefit of
both these two properties. 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
[9] 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 [10 –12] and
has also been considered in conjunction with a spatial representation of spectral type [13]. Later
Donea [14, 15] has used it in deriving a time accurate finite element, but they are not taking the
advantage of sparsity of matrices that are coming in evolutionary problems after applying
high-order wavelet time-accurate schemes and wavelet compression property of smooth data.
Solutions to PDEs often behave differently in different areas. In fluid dynamics we have
shocks, boundary layers, and turbulence. In acoustics an example of PDEs is a low-frequency
wave, with a localized high-frequency burst. For these examples the solution can be smooth in
most of the solution domain, with small areas where the solution changes quickly. Where the
discretization in space (correlatively in time) ought to handle a huge number of degrees of
freedom. The aim of the present article is to examine the feasibility of applying fast wavelet
Taylor-Galerkin method (W-TGM) to such type of PDEs, which require the advantage of
wavelet compression and more regular class of PDEs in higher dimensions. We employ the
generalized time discretization based on Taylor series expansion and wavelet basis for spatial
approximation. This renders our method to be both time and space accurate. We present
wavelet-Taylor Galerkin method for advection equation in one dimension and advection
equation, nonlinear equation, heat equation, and coupled linear system in two dimensions are
considered.
2. WAVELET PRELIMINARIES
Compactly supported wavelets have several properties that are quite useful for representing
solutions of PDEs. The orthogonality, compact support, and exact representation of polynomials
of a fixed degree allow the efficient and stable calculation of regions with strong gradients or
oscillations. Ingrid Daubechies define the class of compactly supported wavelets [1]. Briefly, let
be a solution of the scaling relation
The ak are a collection of coefficients that categorize the specific wavelet basis. The expression
is called the scaling function. The associated wavelet function is defined by the equation
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
276 KUMAR AND MEHRA
共x兲 ⫽ 冘 共⫺1兲 a
k
k
共2x ⫺ k兲.
1⫺k
The normalization 兰 dx ⫽ 1 of the scaling function leads to the condition ¥k ak ⫽ 2. The
translates of are required to be orthonormal, i.e., 兰 (x ⫺ k) (x ⫺ m)dx ⫽ ␦k,m. The scaling
relation implies the condition ¥k⫽0N⫺1
akak⫺2m ⫽ ␦0,m, 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, (2jx ⫺ 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 nonzero, then will have compact support. Since
the translates of the scaling function and wavelet define orthogonal subspaces
The relation
V j⫹1 ⫽ V j ⫹ W j
V 0 傺 V 1 傺 · · · 傺 V j⫹1 ,
V j⫹1 ⫽ V 0 丣 W 0 丣 W 1 丣 · · · 丣 W j .
冘 共⫺1兲
k
⫺k m
k ak ⫽ 0
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
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 277
FIG. 1. Daubechies scaling and wavelet functions for N ⫽ 6 with support on [0, 5].
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
solution methods 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] [1]. As pointed out by Y. Meyer (1990) the complete toll box built
in L2(R) can be used in the periodic case L2([0, 1]) by introducing a standard periodization
technique. This technique consists at each scale in folding, around the integer values, the
wavelet and the scaling functions centered in [0, 1]. It writes ˜共x兲 ⫽ ¥⬁ 共x
j,k j,k j,l n⫽⫺⬁ j,l
Theorem. Let P ⫽ D/2 be the number of vanishing moments for a wavelet j,k and let f 僆
CP(R). Then the wavelet coefficients decay as 兩dj,k兩 ⱕ CP2⫺j[P⫹(1/ 2)]max僆Ij,k兩f (P)()兩.
The first step in the derivation of Taylor-Galerkin schemes [14] is to formulate a high-order
time-stepping scheme algorithm before the discretization of the spatial variable, which is
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
278 KUMAR AND MEHRA
therefore left momentarily continuous. The above procedure shows some resemblance with the
use of the Taylor series in the heuristic analysis of numerical stability by Hirt [17]. 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 that 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 equa-
tion) is obtained that is to be employed after adequate wavelet spatial discretization in the actual
computations. As an illustration of the above procedure, first we are considering W-TGM on
advection equation and 2D nonlinear equation, which is popularly referred to as a counterpart
of one-dimensional (1D) Burgers’ equation, where we are taking the advantage of wavelet
compression.
⭸ t u ⫽ a⭸ x u, (3.1)
␦t2 n ␦t3 n
u n⫹1 ⫽ u n ⫹ ␦ tu nt ⫹ u ⫹ u ··· (3.2)
2 tt 6 ttt
␦ t 2 n⫹1 ␦ t 3 n⫹1
u n ⫽ u n⫹1 ⫺ ␦ tu n⫹1
t ⫹ u ⫺ u ···. (3.3)
2 tt 6 ttt
u n⫹1 ⫺ u n 1 n ␦t n ␦t2 n
⫽ 共u t ⫹ u n⫹1 兲 ⫹ 共u ⫺ u n⫹1
兲 ⫹ 共u ⫹ u n⫹1
ttt 兲, (3.4)
␦t 2 t
4 tt tt
12 ttt
replacing the time derivatives by spatial derivatives, the associated wavelet Taylor-Galerkin
equations based on CN time gives W-TGM scheme
Au n⫹1 ⫽ Bu n , (3.5)
where A ⫽ I ⫺ (a␦t/2)⭸x ⫹ (a2␦t 2/4)⭸2x ⫺ (a2␦t 2/6)⭸2x and B ⫽ I ⫹ (a␦t/2)⭸x ⫹ (a2␦t 2/4)⭸2x ⫺
(a2␦t 2/6)⭸2x . Now wavelet Galerkin discretization turns the problem into a finite dimensional
space.
d un⫹1 ⫽ Ꮽ ⫺1 Ꮾd un ⫽ Ᏸd un . (3.6)
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 279
FIG. 2. Number of coefficients in the successive powers Dn (a) CN times stepping in wavelets and in
finite differences, versus x ⫽ 2n, n ⫽ 15, N ⫽ 1024, ␦t ⫽ 10⫺3, ⑀M ⫽ 10⫺8 (b) Taylor-Galerkin approach
in wavelets and in finite differences.
In this finite dimensional space un is to be replaced by the vector dnu along a wavelet finite basis,
and A and B are replaced by respectively Ꮽ and Ꮾ (finite) matrices. Due to second- and
third-order term in Taylor series our scheme leads to implicit method that needs inversion. Now
to solve equation (3.5) in wavelet basis we will compute Ꮽ⫺1 and Ꮽ⫺1Ꮾ once and store in
compressed form. We can now give a computational procedure for computing (3.5) using
wavelet compression.
Algorithm 1
1: trunc(Ꮽ⫺1, ⑀M) 3 (Ꮽ⫺1)⑀M
2: compute initial guess in wavelet basis 3 d0u
3: trunc(d0u, ⑀V) 3 (d0u)⑀V
For n ⫽ 0, 1, . . . , n1 ⫺ 1
4: (Ꮽ⫺1)⑀MᏮ(dnu)⑀V 3 dn⫹1 u
5: trunc(dn⫹1
u , ⑀V) 3 (dn⫹1
u )⑀V
where trunc(du, ⑀V) ⫽ {dk, 兩djk兩 ⬎ ⑀V} and trunc(Ꮽ, ⑀M) ⫽ {[Am,n] [Am,n], ⬎ ⑀M}.
j
A further property of the wavelet representation of operators is that the successive powers Ᏸn
of the time iteration matrix become sparser and sparser with increasing n. This property is very
specific to wavelets, as the opposite occurs with finite difference, where Ᏸn becomes a more and
more dense matrix as shown in Fig. 2. It is seen from Fig. 2 that in wavelet Taylor-Galerkin
approach compression in the matrix Dn is larger than wavelet Galerkin approach. From this
property we can obtain iterative speed of the wavelet Taylor-Galerkin scheme.
Algorithm 2
⑀M
1: Initialize (Ꮽ⫺1
0 ) and (d0u)⑀V
⑀M ⫺1 ⑀M
2: (Ᏸ0) 3 (Ꮽ0 ) Ꮾ
For n ⫽ 0, 1, . . . , n1 ⫺ 1
3: (Ᏸn)⑀M(dnu)⑀V 3 (dn⫹1u )⑀V
4: Ᏸn 3 Ᏸn⫹1. Then the approximate solution of PDE is at t ⫽ 2n␦t is d u共2 兲.
2 n
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
280 KUMAR AND MEHRA
Since differential operators are local operators, it seems that not much can be gained by
compression. But in wavelet basis it is possible to efficiently invert the differential operator and
then approximate (in a compressed form) the dense evolution operators. There is no need to
change from classical to wavelet coordinates till some time steps. In classical coordinate, the
evolution operator changes from very sparse to dense. In the wavelet representation we may start
the squaring in the classical coordinates and change to the wavelet basis at the point where the
wavelet representation is sparser. Thus we have the following algorithm.
Algorithm 3
1: For n ⫽ 0, 2, . . . , p
2: ( A)⫺1Bun 3 un⫹1
3: Initialize (Ꮽ⫺1)⑀V and (dpu)⑀V
For n ⫽ p ⫹ 1, p ⫹ 2, . . . , p ⫹ n1 ⫺ 1
4: (Ꮽ⫺1)⑀MᏮ(dnu)⑀V 3 (dn⫹1 u )
n⫹1 ⑀V
5: trunc(dn⫹1u , ⑀ V) 3 (du )
It is essential for the success of these algorithms that the computation of the matrix vector
product fully exploits the compressed form of both matrix and vectors. This can be done using
the algorithm in [16] or fast multiplication based on a general sparse format for both matrix and
vector. We can generalize W-TGM based on any other low-order time stepping, as an
illustration of above procedure we will do on the heat equation in two dimensions.
3.1.1. Theoretical Stability of the Linearized Schemes We use the notion of asymptotic
stability of a numerical method as it is defined in [9] for a discrete problem of the form du/dt ⫽
Lu, where L is assumed to be diagonalization matrix. Because for many evolution equations it
is necessary to adapt the time steps to the spatial resolution in order to maintain the stability and
precision of the numerical scheme. The region of absolute stability of a numerical method is
defined for the scalar model problem du/dt ⫽ u to be set of all ␦t such that 储un储 is bounded
as t 3 ⬁. Finally we say that a numerical method is asymptotic stable for a particular problem
if, for small ␦t ⬎ 0, the product of ␦t times every eigenvalues of L lies within the region of
absolute stability.
the domain except near x ⫽ 0.5, where we have a spike and the thresholded wavelet expansion
of the solution u(x, t) is shown in Fig. 3. This wavelet expansion will have few coefficients,
except for in the neighborhood of the spike at x ⫽ 0.5 ⫺ t. We stepped forward to t ⫽ 0.3, where
the solution for j ⫽ 8 and j ⫽ 9 is shown in Fig. 4 using the wavelet Galerkin method (WGM).
Here oscillation are produced throughout the domain due to the quick change of solution near
spike for j ⫽ 8 and oscillation are confined near the neighborhood of spike for j ⫽ 9. In terms
of finite difference methods, we want to have many points in areas where the solution has strong
variation and few points in area where the solution is smooth. If we use a Galerkin method, this
corresponds to the representation of the solution having fewer basis functions in the smooth
areas. Note that by thresholding a wavelet representation we have a way to automatically find
a sparse representation of smooth part. M. Holmstrom and J. Walden have applied adaptive
wavelet methods on such type of PDEs [18, 19]. But by our approach of W-TGM scheme we
are taking the advantage of time-accurate scheme as well as wavelet capabilities of compression
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 281
FIG. 3. (a) Initial function u0(x) (b) Truncated approximate solution by D6 with the threshold ⑀ ⫽ 0.001,
which leads to 13 retained wavelet coefficients, out of original 28.
to produce fast algorithm based on fast matrix vector product in terms of sparsity. We are getting
oscillation free solution as shown in Fig. 5 with and without truncation, where WGM is unable
to produce oscillation free solution near the spike with the same degree of freedom.
Error in approximated solution for WGM and W-TGM without truncation is shown in Fig.
6. Where near the spike error in W-TGM, scheme is small. From this we can conclude that near
the sharp gradients we can take the advantage of time accuracy and compression properties of
wavelet in W-TGM scheme. For higher-order wavelet improvement of error in W-TGM can be
seen in Fig. 6, where WGM shows no improvement. These nice properties are observed for all
the experiments in two dimensions also.
For some ⑀M ⬎ 0, ⑀v ⬎ 0, the density of Ᏸn is small but the truncation error 储Ᏸn ⫺ D⑀nM储 ⬍
C⑀M and 储un ⫺ u⑀nv储 ⬍ C⑀v is bounded for large n. Figure 7 gives an indication of the error
committed by truncation of the solution. In all the schemes ␦t is chosen such that it should
satisfy stability condition. The ␦t for the W-TGM scheme is plotted in Fig. 8.
As in the 1D case, we first examine a linear advection equation. Again it will test the method’s
FIG. 6.
KUMAR AND MEHRA
⑀v ⫽ 10⫺6.
282
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 283
u t ⫽ u x ⫹ u y. (3.7)
Here we are using the third-order accurate W-TGM scheme based on (3.4) and for the third term
we pose (1/2)(unt ⫹ un⫹1
t ) ⫽ (un⫹1 ⫺ un)/␦t. Let ut ⫽ ⵜu and time derivative is utt ⫽ ⵜ2u, putting
these values in equation (3.4):
Au n⫹1 ⫽ Bu n , (3.8)
where A ⫽ I ⫺ (a␦t/2)ⵜu ⫹ (a2␦t 2/4)ⵜ2 ⫺ (a2␦t 3/6)ⵜ2 and B ⫽ I ⫹ (a␦t/2)ⵜu ⫹ (a2␦t 2/4)ⵜ2 ⫺
(a2␦t 3/6)ⵜ2. Unfortunately for a N ⫻ N space discretization, this will result in a system of
algebraic equation that is sparse but is of the order N 2 ⫻ N 2. To overcome this problem, we are
taking full advantage of tensorial wavelet bases that leads to the factorization approximation of
2D operator.
The approximation space used here is VJ. Then our W-TGM scheme express themselves as
products of N ⫻ N matrices:
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
284 KUMAR AND MEHRA
d un⫹1 ⫽ Bd un B T ⫹ ␦ tA ⫺1 d f A ⫺1 .
T
(3.10)
Therefore all the algorithms in two dimensions will also take the advantage of wavelet
compression like 1D case.
3.2.1. Numerical Simulation of 2D Advection Equation The2 initial function is the smooth
function with a localized spike u0 (x, y ⫽ exp⫺␣(( x⫺1/ 2) ⫹( y⫺1/ 2) ) ⫺ 0.2 sin(2x)sin(2y) and
2
periodic boundary conditions. In Fig. 9 the initial function and solution at t ⫽ 0.3 without and
with truncation is shown.
u t ⫹ u共u x ⫹ u y 兲 ⫽ 共u xx ⫹ u yy 兲. (3.11)
Here we are using third-order accurate W-TGM scheme based on (3.4) and for the third term we
pose (1/2)(unt ⫹ un⫹1
t ) ⫽ (un⫹1 ⫺ un)/␦t. Let f(x, y) ⫽ u(ux ⫹ uy); then the original equation is
ut ⫽ ⌬u ⫺ f(x, y) and time derivative is utt ⫽ 2⌬2u ⫺ ⌬f, putting these values in equation
(3.4)
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 285
FIG. 9. (a) The initial function; (b) solution at t ⫽ 0.3 when j ⫽ 8, ⑀M ⫽ 10⫺5, ⑀v ⫽ 0; (c) solution at
t ⫽ 0.125 for ⑀M ⫽ 10⫺5, ⑀v ⫽ 10⫺4.
where A ⫽ I ⫺ (␦t/2)⌬u ⫹ (2␦t 2/4)⌬2 ⫺ (2␦t 3/6)⌬2 and B ⫽ I ⫹ (␦t/2)⌬u ⫹ (2␦t 2/4)⌬2 ⫺
(2␦t 3/6)⌬2. Here we are also using factorization approximation of 2D operator. Then our
W-TGM scheme express themselves as products of N ⫻ N matrices:
d un⫹1 ⫽ Bd un B T ⫹ ␦ tA ⫺1 d f A ⫺1 .
T
(3.13)
FIG. 10. (a) The initial function; (b) solution at t ⫽ 0.125 when j ⫽ 7, ⑀M ⫽ 10⫺5, ⑀v ⫽ 0; (c) solution
at t ⫽ 0.125 for ⑀M ⫽ 10⫺5, ⑀v ⫽ 10⫺4.
AV ⫽ F, (3.14)
where V ⫽ un⫹1 ⫺ un is the vector of nodal unknowns. F ⫽ F(un) is a known vector and A 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.15)
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 287
where L is the diagonal and positive matrix obtained from A by row-sum technique.
L ii ⫽ 冘A ,
i
ii L ij ⫽ 0, j ⫽ i. (3.16)
where
Then, under the assumption 储X储 ⱕ 1, the inverse of A in the form (3.17) can be expressed by the
following series,
Truncating the series after the term X gives the following two-term approximation of the inverse
of A,
A ⫺1 共2兲 ⫽ 2L ⫺1 I ⫺ 冉 1
2 冊
AL ⫺1 , (3.20)
冉
A ⫺1 共3兲 ⫽ 3L ⫺1 I ⫺ AL ⫺1 ⫹
1
3
AL ⫺1 AL ⫺1 冊 (3.21)
and so on. The successive approximations in the above approximate factorization technique can
be generated by the following multipass algorithm. Consider the sequence of approximate
solutions V ( g), g ⫽ 0, 1, . . . , G, defined as 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
Finally assume V ⫽ V (G). The approximation (3.20) and (3.21) 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.
Remark. The W-TGM based on and CN time stepping requires to solve at each time step
system of the form AV ⫽ F, where V ⫽ un⫹1 is the vector of nodal unknowns.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
288 KUMAR AND MEHRA
3.4.1. Scheme Based on Euler Time Stepping To this aim the following Taylor-series
expansion in the time increment ⌬t is used:
1 1
u n⫹1 ⫽ u n ⫹ ␦ tu nt ⫹ 共 ␦ t兲 2 u ntt ⫹ 共 ␦ t兲 3 u nttt ⫹ · · · . (3.23)
2 6
The time derivatives in the right-hand side of (3.23) are then replaced by their expressions
evaluated formally from the governing equation:
⭸u
⫽ F共x, y兲, (3.24)
⭸t
where F( x, y) ⫽ ⌬u ⫹ f(x, y) and ⌬ ⫽ (⭸2/⭸x2) ⫹ (⭸2/⭸ y2) is a 2D Laplacian operator and u
is the unknown scalar function defined on a bounded square domain with the periodic boundary
condition. Equation (3.24) is first discretized in time according to Taylor series (3.23), which
includes second- and third-order time derivatives.
combining equation (3.23) and (3.25) and replacing unt by (un⫹1 ⫺ un)/⌬t, the equation (3.24)
is replaced by the following expansion:
冉 1⫺
2␦ t 2 2
6
⌬ 冊冉 u n⫹1 ⫺ u n
␦t 冊⫽ ⌬u ⫹
n
2 ␦ t 2 n ␦ t
2
⌬u ⫹
2
⌬f ⫹ f. (3.26)
The fully discretized equations are subsequently obtained by means of the WGM.
u ttt ⫽ 3 ⌬ 3 u ⫹ 2 ⌬ 2 f, (3.27)
u n⫹1 ⫺ u n 2␦ t 2 n 3␦ t 2 3 n 2␦ t 2 2 ␦ t
⫽ ⌬u n ⫹ ⌬u ⫹ ⌬u ⫹ ⌬ f⫹ ⌬f ⫹ f. (3.28)
␦t 2 6 6 2
for ⫽ 1/2 ⫹ ␣␦t, where ␣ is a small positive constant. This method damps all component of
the solution and is formally second-order accurate in time. Here we are deriving W-TGM based
on time stepping. For the case ␣ ⫽ 0 we will get W-TGM based on CN time stepping. To
obtain an improved order of accuracy in ␦t we shall apply a Taylor-Galerkin method on the
following Taylor series:
1 1
u n⫹1 ⫽ u n ⫹ ␦ tu nt ⫹ 共 ␦ t兲 2 u ntt ⫹ 共 ␦ t兲 3 u nttt ⫹ O共 ␦ t 4 兲
2 6
1 1
u n ⫽ u n⫹1 ⫺ ␦ tu n⫹1
t ⫹ 共 ␦ t兲 2 u n⫹1
tt ⫺ 共 ␦ t兲 3 u n⫹1
ttt ⫹ O共 ␦ t 兲.
4
(3.30)
2 6
u n⫹1 ⫺ u n ␦t
关共1 ⫺ 兲 ⫹ 兴 ⫽ 关共1 ⫺ 兲u nt ⫹ u n⫹1 兲] ⫹ 关共1 ⫺ 兲u ntt ⫺ u n⫹1
tt 兴
␦t t
2
␦t2
⫹ 关共1 ⫺ 兲u nttt ⫹ u n⫹1
ttt 兴. (3.31)
6
u n⫹1 ⫺ u n
关共1 ⫺ 兲u nt ⫹ u n⫹1 兴⫽ . (3.32)
t
␦t
冉 1⫺
6 冊
2 ␦ t 2 2 u n⫹1 ⫺ u n
⌬
␦t
⫽ 关共1 ⫺ 兲 ⌬u n ⫹ ⌬u n⫹1 兲]
␦t ␦ t
⫹ 关共1 ⫺ 兲 2 ⌬ 2 u n ⫺ 2 ⌬ 2 u n⫹1 兴 ⫹ 共1 ⫺ 2 兲 ⌬f ⫹ f. (3.33)
2 2
The fully discretized equations are subsequently obtained by means of the WGM.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
290 KUMAR AND MEHRA
u n⫹1 ⫺ u n ␦t
关共1 ⫺ 兲 ⫹ 兴 ⫽ 关共1 ⫺ 兲 ⌬u n ⫹ ⌬u n⫹1 兲] ⫹ 关共1 ⫺ 兲 2 ⌬ 2 u n ⫺ 2 ⌬ 2 u n⫹1 兴
␦t 2
␦t2 2␦ t 2 2 ␦ t
⫹ 关共1 ⫺ 兲 3 ⌬ 3 u n ⫹ 3 ⌬ 3 u n⫹1 兴 ⫹ ⌬ f ⫹ 共1 ⫺ 2 兲 ⌬f ⫹ f. (3.34)
6 6 2
The fully discretized equations are subsequently obtained by means of the WGM.
3.4.3. W-TGM Based on Leap-Frog Time Stepping The Taylor-Galerkin based on leap-
frog time stepping can be obtained by the following Taylor series expression:
1 1
u n⫹1 ⫽ u n ⫹ ␦ tu nt ⫹ 共 ␦ t兲 2 u ntt ⫹ 共 ␦ t兲 3 u nttt ⫹ · · ·
2 6
1 1
u n⫺1 ⫽ u n ⫺ ␦ tu nt ⫹ 共 ␦ t兲 2 u ntt ⫺ 共 ␦ t兲 3 u nttt ⫹ · · · (3.35)
2 6
from which we deduce an improved approximation to the temporal derivative at time t ⫽ n␦t
in the form
u n⫹1 ⫺ u n⫺1 ␦ t 2 n
u nt ⫽ ⫺ u . (3.36)
2␦t 6 ttt
冉 1⫺
2␦ t 2 2
6
⌬ 冊冉 u n⫹1 ⫺ u n⫺1
2␦t 冊
⫽ ⌬u n ⫹ f. (3.37)
The fully discretized equations are subsequently obtained by means of the WGM.
u n⫹1 ⫺ u n⫺1 3␦ t 2 3 n 2␦ t 2 2
⫽ ⌬u n ⫹ ⌬u ⫹ ⌬ f ⫹ f. (3.38)
2␦t 6 6
The fully discretized equations are subsequently obtained by means of the WGM.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 291
3.4.4. Numerical Simulation of Heat Equation Let’s consider equation (3.24) on a doubly
periodic square, with vanishing (u(x, y, 0) ⫽ 0) initial condition and the following forcing term:
f 共x, y兲 ⫽ sin共2x兲cos共2y兲.
1
u共x, y, t兲 ⫽ 共1 ⫺ e ⫺8 t 兲sin共2x兲cos共2y兲.
2
8 2
We present numerical results obtained with W-TGM, W-TGMS, and 3-pass explicit W-TGM
for the 2D heat equation. Figure 11(a) illustrates the exact solution and Fig. 11(b) shows the
wavelet solution using Daubechies D6 scaling function on a mesh with n ⫽ 25 unknown in each
direction at time t ⫽ 0.005. All the results we present are obtained using Daubechies D6 scaling
function.
In Table I errors in L⬁ norm of the solution obtained by W-TGM, W-TGMS, and 3-pass
explicit W-TGM based on Euler scheme are presented. From Table I we can observe that the
error is not affected by decreasing ␦t beyond 10⫺4 and this is also true with other schemes.
TABLE I. Accuracy of results given by the L⬁ norm using Euler time-stepping scheme.
j ␦t W-TGM W-TGMS W-TGM (explicit)
TABLE II. Accuracy of results given by the L⬁ norm using time-stepping scheme.
j ␦t W-TGM W-TGMS W-TGM (explicit)
⫺4
3 10 0.0367 0.0367 0.0399
4 10⫺4 0.0049 0.0049 0.0194
5 10⫺4 0.0024 0.0024 0.0024
In Tables II and III errors in L⬁ norm of the solution obtained by W-TGM, W-TGMS, and
3-pass explicit W-TGM based on and CN scheme are presented.
In Table IV errors in L⬁ norm of the solution based on leap-frog scheme are presented.
One could argue that the presence of consistent mass matrix in our W-TGMs represent a
serious disadvantage from the point of view of the computational efficiency of the method.
Therefore the 3-pass explicit scheme is particularly useful for the solution of the algebraic
system resulting from W-TGM. As shown by the last columns of Tables I–IV, the 3-pass explicit
W-TGM are practically never less accurate than their implicit parents W-TGM. 储u储L⬁ ⫽
maxk⫽0,1,. . . ,2j⫺1兩u(k/(2j)兩.
⭸u
⭸t 冋
⫺2 1 ⭸u
⫽ 1 ⫺2
⭸x
⫹ 0 ⫺␣ 册
⫺1 0 ⭸u
⭸y
, 冋 册 (3.39)
⭸u 2 ⭸u 1 ⭸u 2 ⭸u 2
⫽ ⫺2 ⫺␣ . (3.40)
⭸t ⭸x ⭸x ⭸y
Here we consider the simple case of a second-order scheme obtained by dropping the last term
of the series in (3.23):
u 1n⫹1 ⫺ u 1n ␦t
⫽ u 1 tn ⫹ u 1 ttn ⫹ O共 ␦ t 2 兲
␦t 2
u 2n⫹1 ⫺ u 2n ␦t
⫽ u 2 tn ⫹ u 2 ttn ⫹ O共 ␦ t 2 兲. (3.41)
␦t 2
TABLE III. Accuracy of results given by the L⬁ norm using Crank-Nicolson time-stepping scheme.
j ␦t W-TGM W-TGMS W-TGM (explicit)
TABLE IV. Accuracy of results given by the L⬁ norm using leap-frog time-stepping scheme.
j ␦t W-TGM W-TGMS W-TGM (explicit)
⫺4
3 10 0.0375 0.0375 0.0375
4 10⫺4 0.005 0.005 0.005
5 10⫺4 0.0025 0.0037 0.0025
The substitution of u1t, u2t and u1tt, u2tt in Taylor expansion (3.41) yields the semi-discrete
equations:
u 1n⫹1 ⫺ u 1n
␦t
⫽ ⫺2
⭸u1n ⭸u2n ⭸u1n ␦t ⭸2u1n
⭸x
⫹
⭸x
⫺
⭸y
⫹
2 冉
⭸2u2n ⭸2u1n
5 2 ⫺4 2 ⫹ 2 ⫹4
⭸x ⭸x ⭸y
⭸2u1n
⭸x⭸y
⫺ 共1 ⫹ ␣兲
⭸2u2n
⭸x⭸y 冊
u 2n⫹1 ⫺ u 2n ⭸u 1n ⭸u 2n ⭸u 2n ␦ t
⫽ ⫺2 ⫺␣ ⫹
␦t ⭸x ⭸x ⭸y 2
⫻ 冉 ⫺4
⭸2u1n
⭸x 2 ⫹5
⭸2u2n
⭸x 2 ⫺ 共1 ⫹ ␣兲
⭸2u1n
⭸x⭸y
⭸2u1n
⫹ ␣2 2 ⫹ 4␣
⭸y
⭸2u2n
⭸x⭸y
. 冊 (3.42)
The fully discretized equations are subsequently obtained by means of the WGM.
u 1 ⫽ sin共x ⫺ t兲 ⫹ sin共y ⫺ t兲
FIG. 12. Elevation of u at time t ⫽ 0.005 (a) analytical solution (b) by second order accurate W-TGM
scheme.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
294 KUMAR AND MEHRA
Numerical results of W-TGM scheme for the coupled linear equations are shown in Fig. 12.
4. CONCLUSION
Taylor generalized time-stepping schemes in conjunctions with wavelet basis for spatial ap-
proximation are employed in deriving a new class of space-time accurate methods known as
W-TGMs for PDEs. The concepts are introduced through linear advection equation and 2D
nonlinear equation, a counterpart of a 1D Burgers’ equation, where we can show the power of
wavelet compression as well as time-accurate scheme and more regular class of PDEs like heat
equation, coupled linear equations. Factorization of the 2D operator, taking the advantage of
tensorial wavelet bases and multipass explicit schemes, based on approximate factorization, for
linear systems provide computationally efficient and accurate solution for higher dimension
PDEs where matrix inversion is demanding. Wavelet-Taylor Galerkin schemes with spatial
correction (W-TGMS) and with mixed spatial and temporal corrections (W-TGM) are exam-
ined and found to be of same order of accuracy in 2D problems. The results are quite
encouraging and we are extending our study to 3D and more complex phenomena like red spot
in Jupiter, giant ocean waves [20], turbulence where spurious oscillations are more serious when
a significant gradient is involved.
References
1. I. Daubechies, Orthonormal basis of compactly supported wavelets, Comm Pure Appl Math 41 (1988),
906 –966.
2. J. Liandrat, V. Perrier, and P. Tchmitchian, Numerical resolution of non-linear partial differential
equations using wavelet approach, Wavelets and Applications, Boston, MA, 1992, pp. 227–238.
3. A. Latto and E. Tenenbaum, Les ondelettes à support compact et la solution de l’èquation de Burgers,
CR Acad Sci Paris série I (1990), 903–909.
4. J. Weiss, Wavelets and the study of two dimensional turbulence, In: Proceedings of the USA-French
workshop on wavelets and turbulence, Princeton University, June 1991.
5. M. Farge, N. Kevlahan, V. Perrier, and E. Goirand, Wavelets and turbulence, Proceedings of the IEEE,
1996, p. 84.
6. J. Frohlich and K. Schneider, An adaptive-vaguelette algorithm for the solution of PDEs, J Comp Phys
130 (1997), 174 –190.
7. R. Glowinski, W. Lawton, M. Ravachol, and E. Tenenbaum, Wavelet solutions of linear and
non-linear elliptic, parabolic and hyperbolic problems in 1D, Computing methods in applied sciences
and engineering, SIAM, Philadelphia, PA, 1990, chapt. 4, pp. 55–120.
8. B. Engquist, S. Osher, and S. Zhong, Fast wavelet based algorithms for linear evolution equations,
SIAM J Sci Comp 15(4) (1994), 755–775.
9. J. D. Lambert, Computational methods for ordinary differential equations, Wiley, London, 1973.
10. P. D. Lax and B. Wendroff, Systems of conservation laws, Comm Pure Appl Math 13 (1960), 217.
11. P. D. Lax and B. Wendroff, On the stability of difference schemes, Comm Pure Appl Math 15 (1962),
363.
12. P. D. Lax and B. Wendroff, Difference schemes for hyperbolic equations with high order of accuracy,
Comm Pure Appl Math 17 (1964), 381.
13. J. Gazdag, Numerical convective schemes based on accurate computation of space derivatives,
J Comput Phys 13 (1973), 100 –113.
10982426, 2006, 2, Downloaded from [Link] by <Shibboleth>-member@[Link], Wiley Online Library on [19/12/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME ACCURATE METHOD FOR PDES 295
14. J. Donea, A Taylor-Galerkin method for convective transport problems, Int J Numer Methods Eng 20
(1984), 101–119.
15. J. Donea, L. Quartapelle, and V. Selmin, An analysis of time discretization in finite element solution
of hyperbolic problems, J Comput Phys 70 (1987), 463.
16. O. M. Nilsen, Wavelets in scientific computing, PhD thesis, Technical University of Denmark,
Lyngby, 1998.
17. C. W. Hirt, Heuristic stability theory for finite-difference equations, J Comput Phys 2 (1968),
339 –355.
18. M. Holmstrom, Solving hyperbolic PDEs using interpolating wavelets, SIAM J Sci Comput 21(2)
(1999), 405– 420.
19. M. Holmstrom and J. Walden, Adaptive wavelets methods for hyperbolic PDEs, J Sci Comput 13(1)
(1998), 19 – 49.
20. B. V. Rathish Kumar and Mani Mehra, Time-accurate solutions of Korteweg-de Vries equation using
wavelet Galerkin method, Appl Math Comput 162(1) (2004), 447– 460.