Abstract

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.

1. Introduction

In this paper, we seek a numerical solution of linear parabolic partial differential equation (PDE) in one space dimension given by 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 for and . PDE (1) is classified as follows.(i)If and , PDE (1) is called diffusion (or heat) equation.(ii)If and , PDE (1) is called diffusion-reaction equation.(iii)If and , PDE (1) is called convection-diffusion equation.(iv)If and , PDE (1) is called convection-diffusion-reaction equation.

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–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–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 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.

2. 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].

2.1. Basic Properties of Daubechies Wavelets

For a positive integer , consider two functions defined by where and are two specific sequences [8] such that for . The functions and are called scaling function and wavelet function, respectively, and is called their order. These functions are compactly supported with and they are available in wavelet toolbox of MATLAB for . They satisfy the following properties: The integer dilates and translates of and are defined as for .

Now, for all , define Then we have in the sense that where is the -orthogonal projection of onto .

2.2. Approximation of Function Spaces Using Wavelet Integrals

We mention here that for an open interval and for an integer , 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 is the space of all square (Lebesegue) integrable functions which is also a Hilbert space with inner product 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 such that and the Sobolev space can be approximated by the restrictions of translates and dilates of to [4]. Hence can be approximated by their integrals.

Now we have, supp. Let be any integer and let be defined by as in (6). Now we define to be the set of restrictions of all functions in to . In fact, we take We shift the support of from to by using the transformation where and .

Now, considering instead of , 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( and a basis for can be constructed as Now, for , we define the spaces(i), where (ii), where

and(iii), where Now by Proposition  4.2 in [12], the spaces , , and are finite dimensional subspaces of , , and , respectively, and the sets are bases for , , and , respectively.

Figure 1 shows the basis functions for , , and for and .

(a)

(b)

(c)

(a)
(b)
(c)

Figure 1 

Basis functions for , , and for .

The derivatives of the above basis functions are the interpolated scaling functions themselves, with a little difference, as follows:(i), for ,(ii), for ,(iii), for .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.

For any (resp., or ) and (resp., or ), we have where , , and is a positive constant depending on and .

Now, by Proposition  4.3 in [12], we can find a finite dimensional subspace of given by a basis for which is given by

Remark 1. The dimensions of the space are , that of and are each and of is .

3. 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 , , , or (according to boundary conditions (26)(i), (26)(ii), (26)(iii), and (26)(iv), resp.), integrating by parts with respect to in , and using the natural boundary conditions, we get(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

4. Spatial Approximation

In [12], it is established that is sufficient for the solution of second-order (in space) problems. Let be any integer and let be the scaling function. Let , , be the approximating subspaces of the semidiscrete space approximations of the variational problems formulated above. Considering the bases of =, , , and constructed in Section 2, the approximate solution of variational problems (27), (28), (29), and (30) can be taken as where  + and , for problem (27),  and , for problem (28),  and , for problem (29),  and , for problem (30).Applying Galerkin method to each of the variational problems (27), (28), (29), and (30) with the approximate solution (31), we get a system of first-order ordinary differential equations in : where is the stiffness matrix given by is the mass matrix given by and is the force vector given by(i), for problem (27),(ii), for problem (28),(iii), for problem (29),(iv), for problem (30), where The solution of the system of ordinary differential equations (32) gives the complete solution of the original problem. The value of the vector at is determined from the system of linear equations obtained by multiplying (25) by 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).

5. 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 of . (In [16], is .)

Using approximation (37) for time and in (32), we obtain Rearranging the terms, we get where The solution at time is obtained in terms of the solution at time by inverting the matrix . The solution at time is , obtained from (36).

5.2. Stability and Accuracy

Error is introduced into the solution 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 , 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 in (39) is diagonal. Otherwise, it is called an implicit scheme.

For different values of , the time discretization scheme (37) is classified as follows.(1)The forward difference (or Euler) scheme: ; explicit and conditionally stable; order of accuracy = .(2)The backward difference scheme: ; implicit and unconditionally stable; order of accuracy = .(3)The Galerkin scheme: ; implicit and unconditionally stable; order of accuracy = .(4)The Crank-Nicholson scheme: ; implicit and unconditionally stable; order of accuracy = .

6. Numerical Experiments

Here the methodology for the solution of IBVP (24)–(26) described above is demonstrated with examples of diffusion, diffusion-reaction, convection-diffusion, and convection-diffusion-reaction equations with all the four types of boundary conditions, using all the four time discretization schemes. We compare the computed solutions with the exact solutions and obtain and norm errors. The stiffness matrices for diffusion and diffusion-reaction equations are symmetric and those for convection-diffusion and convection-diffusion-reaction equations are not symmetric. All the systems of linear equations are solved by using Gaussian elimination method and the integrations are carried out by Simpson’s quadratures. All the computations are performed in MATLAB.

6.1. Test Problem (42)

Here we solve a diffusion equation () with the first mixed boundary conditions (MBC1). Here we take , the space domain , the initial value function , and the boundary value functions . If we solve this problem by the method of separation of variables, we get the exact solution as

For , 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 . Figure 2 shows the exact solution and wavelet solution for (unstable) and (stable) at .

(a)

(b)

(a)
(b)

Figure 2 

Test problem (42) . wavelet solution and exact solution at .

For , 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 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 .

Figure 3 

Test problem (42) (). Convergence of 2( wavelet solutions with decreasing time steps at .

For , wavelet is used for spatial discretization. This scheme is also unconditionally stable. Figure 4 shows the computed solutions at times , , , , and with 0.05.

Figure 4 

Test problem (42) (, ). wavelet solutions at 5, 10, 15, 20, and 25 time steps.

For , we use and wavelets for spatial discretization. Table 1 shows , , and norm errors at for with various values of .

6.2. Test Problem (43)

Here we solve a diffusion-reaction equation () with the second mixed boundary conditions (MBC2). Here we take , and the space domain . The right-hand function , the initial value function , and the boundary value functions and are obtained as per the exact solution:

For , wavelet is used for spatial discretization. This scheme is conditionally stable. The largest eigenvalue of the associated problem is 235.57. Therefore, the maximum time step is . Figure 5 shows the exact solution and wavelet solution for (unstable) and (stable) at .