Wavelet Method for Numerical Solution of Parabolic Equations

We derive a highly accurate numerical method for the solution of parabolic partial differential equations in one space dimension using semidiscrete approximations. The space direction is discretized by wavelet-Galerkin method using some special types of basis functions obtained by integrating Daubechies functions which are compactly supported and differentiable. The time variable is discretized by using various classical finite difference schemes. Theoretical and numerical results are obtained for problems of diffusion, diffusion-reaction, convection-diffusion, and convection-diffusion-reaction with Dirichlet, mixed, and Neumann boundary conditions. The computed solutions are highly favourable as compared to the exact solutions.


Introduction
In this paper, we seek a numerical solution of linear parabolic partial differential equation (PDE) in one space dimension given by − ( ) + + = ( , ) , < < , > 0, with initial condition and with boundary conditions of Dirichlet, mixed, or Neumann type, where , , and are functions of and , in general. Here we assume that > 0 for < < and > 0. PDE (1) is classified as follows. These equations have applications in a number of fields, for example, heat transfer, fluid mechanics, ground water pollutants, financial mathematics, and so on.
Several methods exist for the solution of parabolic problems, for example, [1][2][3][4][5][6][7]. But still there is a need for modification of the solution methodology in case of (i) Neumann and mixed boundary conditions, (ii) time-dependent boundary conditions, and (iii) time-dependent source term . The algorithm in this paper is suitable for such situations.
In the present paper, we apply semidiscrete approximation for the solution of initial-boundary value problem (IBVP) governed by the PDE (1). Here we use wavelet-Galerkin (variational) method to discretize the space direction and the time variable is discretized by using various classical finite difference schemes. Wavelets in consideration here are Daubechies compactly supported wavelets [8] which are differentiable.
Wavelet applications to the solution of PDE problems are relatively new. Some recent applications are [5,6,[9][10][11][12] among many more. To discretize a PDE problem by wavelet-Galerkin method, the Galerkin bases are constructed from orthonormal bases of compactly supported wavelets. This can be done in a number of ways. If wavelets are used directly in such construction such as one in [4], the resulting basis functions do not satisfy various types of boundary conditions. The basis functions obtained by integrating wavelet functions in a suitable manner can be 2 Journal of Computational Engineering made to satisfy various types of boundary conditions; see Deka and Choudhury [12], which is a generalization of [13] for general domains. Also, this approach makes one free from computing connection coefficients [14] for the solution of higher dimensional problems and is convenient to use low order wavelets resulting in low computational efforts. The present paper is a generalization of [12] to the solution of parabolic problems. Basis functions satisfying various types of boundary conditions can also be obtained without integrating wavelet functions; see Choudhury and Deka [10]. For the demonstration of the algorithm presented in the paper, some numerical results are computed which show that the computed solutions compare highly favourably with the exact solutions. All the numerical tests are performed in MATLAB. The paper is organized as follows.
Section 2 describes wavelet integral functions and their approximation properties elaborately including a short introduction to wavelet functions. In Section 3, we formulate the problems to be solved. Sections 4 and 5 discuss the spatial and time approximations of the problems using wavelet-Galerkin method and classical finite difference schemes, respectively. In Section 6, we perform some numerical tests to demonstrate the method presented. Section 7 contains the necessary conclusion.

Wavelet Integrals and Their Approximation Properties
In this section, we discuss the approximations of various function (Sobolev) spaces in finite dimension, useful for wavelet-Galerkin solutions of PDE problems, using integrals of Daubechies scaling functions (the wavelet integrals) [12], after a short account of Daubechies scaling and wavelet functions. For details about Daubechies functions, we refer [8,15].

Basic Properties of Daubechies Wavelets.
For a positive integer , consider two functions , ∈ 2 (R) defined by where { } ∈Z and { } ∈Z are two specific sequences [8] such that = = 0 for ∉ {0, 1, . . . , }, = 2 − 1. The functions and are called scaling function and wavelet function, respectively, and is called their order. These functions are compactly supported with supp( ) = supp( ) = [0, ] and they are available in wavelet toolbox of MATLAB for 1 ⩽ ⩽ 45. They satisfy the following properties: The integer dilates and translates of and are defined as for , ∈ Z. Now, for all ∈ Z, define Then we have in the sense that where is the 2 -orthogonal projection of 2 (R) onto .

Approximation of Function Spaces Using Wavelet Integrals.
We mention here that for an open interval ( , ) and for an integer ⩾ 1, the space is called the Sobolev space of order , which is a Hilbert space with inner product and associated norm ‖ ⋅ ‖ , where the derivatives are in weak (or distributional) sense. It is to be noted here that 0 ( , ) = 2 ( , ) is the space of all square (Lebesegue) integrable functions which is also a Hilbert space with inner product ⟨ , V⟩ = ∫ V and associated norm ‖ ⋅ ‖. Now, we define the following subspaces of ( , ) which will be used for the solution of PDE problems in this paper: Let be any positive integer and let and be the scaling function and wavelet function, respectively. Then there exists an integer , 0 ⩽ < such that , ∈ (R) and the Sobolev space ( , ) can be approximated by the restrictions of translates and dilates of to ( , ) [4]. Hence +1 ( , ) can be approximated by their integrals. Now we have, supp( ) = supp( ) = [0, ]. Let ⩾ 0 be any integer and let be defined by as in (6). Now we define (0, ) to be the set of restrictions of all functions in to (0, ). In fact, we take We shift the support of from [0, ] to [ , ] by using the transformation where 0 ⩽ ⩽ and ⩽ ⩽ . Now, considering ∈ [ , ] instead of ∈ [0, ], we define the space ( , ) as where Since ( , ) is bounded, the space ( , ) is finite dimensional and is a closed subspace of ( , ). By Proposition 4.1 in [12], dim( ( , )) = 2 + − 1 and a basis for ( , ) can be constructed as Now, for , ∈ ( , ), we define the spaces and Now by Proposition 4.2 in [12], the spaces 0 ( , ), * ( , ), and + ( , ) are finite dimensional subspaces of are bases for 0 ( , ), * ( , ), and + ( , ), respectively. Figure 1 shows the basis functions for 0 (0, 1), * (0, 1), and + (0, 1) for = 2 and = 2.
The derivatives of the above basis functions are the interpolated scaling functions themselves, with a little difference, as follows: By Theorem 4.2 in [12] and a standard duality argument [13], the error estimates of the above approximations in norm can be obtained as follows.

Problem Formulation
Here, we consider the variational formulation of the following parabolic IBVP: with the initial condition: and with the boundary conditions: Here, the boundary conditions where the (space) derivative of is present are natural boundary conditions and the others are essential boundary conditions. Multiplying (24) by a function V ∈ 1 0 ( , ), 1 * ( , ), (i) the weak form of problem (24)-(26)(i) as (ii) the weak form of problem (24)-(26)(ii) as (iii) the weak form of problem (24)-(26)(iii) as and (iv) the weak form of problem (24)-(26)(iv) as Journal of Computational Engineering 5

Spatial Approximation
In [12], it is established that ⩾ 2 is sufficient for the solution of second-order (in space) problems. Let ⩾ 2 be any integer and let be the scaling function. Let , ⩾ 0, be the approximating subspaces of the semidiscrete space approximations of the variational problems formulated above. Considering the bases {Φ , } of = 0 ( , ), * ( , ), + ( , ), and ( , ) constructed in Section 2, the approximate solution of variational problems (27), (28), (29), and (30) can be taken as where = [ ] is the stiffness matrix given by = [ ] is the mass matrix given by and = [ ] is the force vector given by The solution of the system of ordinary differential equations (32) gives the complete solution of the original problem. The value of the vector at = 0 is determined from the system of linear equations obtained by multiplying (25) by V and integrating in ( , ) with respect to , where the vector = [ ] is given by Thus, the IBVP (24)-(26) becomes a single initial value problem (IVP) given by (32) and (36).

Time Discretization
5.1. The Scheme. The most commonly used method for the solution of IVP given by (32) and (36) is a -family of approximation [16]: where refers to the value of at time = = Δ , and there is a similar meaning oḟ. (In [16], is .) Using approximation (37) for time and +1 in (32), we obtain Rearranging the terms, we get The solution +1 at time +1 is obtained in terms of the solution at time by inverting the matrix̂+ 1 . The solution 0 at time = 0 = 0 is (0), obtained from (36).

Stability and Accuracy.
Error is introduced into the solution +1 at each time step when approximation (37) is used. If the error grows unboundedly, the solution scheme is said to be unstable. Otherwise it is said to be stable. Accuracy of a numerical scheme is a measure of the closeness between the approximate and the exact solutions, whereas stability is a measure of the boundedness of the approximate solution with time. The size of the time step Δ can influence both the accuracy and the stability. A smaller time step results in an improved accuracy. A numerical scheme is said to be conditionally stable, if it is stable only when certain restrictions on the time step are satisfied. The above -family of time approximation schemes, for which < 1/2, is stable only if where is the largest eigenvalue of the problem associated with the original PDE problem. A numerical scheme is called an explicit scheme, if the matrix̂+ 1 in (39) is diagonal. Otherwise, it is called an implicit scheme.
For different values of , the time discretization scheme (37) is classified as follows.

Numerical Experiments
Here For = 0, 2( = 0) wavelet is used for spatial discretization. As this scheme is conditionally stable, we have to find an upper limit of the time step Δ using the stability condition (41). The largest eigenvalue of the associated problem is 526.47. Therefore, the maximum time step is 2/526.47 ≈ 0.0038. Figure 2 shows the exact solution and 2( = 0) wavelet solution for Δ = 0.004 (unstable) and Δ = 0.003 (stable) at = 1. For = 1, 2( = 1) wavelet is used for spatial discretization. As this scheme is unconditionally stable, there is no restriction on Δ . However, to obtain a sufficiently accurate solution, Δ must be taken small enough. Figure 3 shows the exact solution and the computed solutions at time = 1 with Δ = 0.5, 0.25, 0.125, and 0.0625. Here it can be seen that the computed solutions converge to the exact solution as Δ → 0. For = 2/3, 2( = 2) wavelet is used for spatial discretization. This scheme is also unconditionally stable. Figure 4 shows the computed solutions at times 5Δ , 10Δ , 15Δ , 20Δ , and 25Δ with Δ =0.05. For = 1/2, we use 2 and 3 wavelets for spatial discretization. Table 1 shows ∞ , 2 , and 1 norm errors at = 1 for = 0, 1, 2, 3 with various values of Δ . ( , ) = − sin .

Conclusion
In this paper, we have derived a method for numerical solution of linear parabolic partial differential equations in one space dimension. In this algorithm, we apply wavelet-Galerkin method in the space direction and some classical finite difference schemes in the time direction. The basis functions for wavelet-Galerkin method have been constructed by integrating Daubechies scaling functions as in [12]. Indeed, wavelet functions can also be used instead of scaling functions [13]. The most important merits of this algorithm are that it is useful for problems with mixed and Neumann boundary conditions besides Dirichlet one. Also the boundary conditions and the source term can be timedependent. The method gives high favourable accuracy when compared to the exact solutions for problems of this types. The matter here left to see is what will happen if wavelets are used to discretize both the space and the time directions. The extension of the algorithm in multidimension is the matter of further research.