Abstract

A numerical approach is proposed for solving multidimensional parabolic diffusion and hyperbolic wave equations subject to the appropriate initial and boundary conditions. The considered numerical solutions of the these equations are considered as linear combinations of the shifted Bernoulli polynomials with unknown coefficients. By collocating the main equations together with the initial and boundary conditions at some special points (i.e., CGL collocation points), equations will be transformed into the associated systems of linear algebraic equations which can be solved by robust Krylov subspace iterative methods such as GMRES. Operational matrices of differentiation are implemented for speeding up the operations. In both of the one-dimensional and two-dimensional diffusion and wave equations, the geometrical distributions of the collocation points are depicted for clarity of presentation. Several numerical examples are provided to show the efficiency and spectral (exponential) accuracy of the proposed method.

1. Introduction

Many of the physical models and chemical reactions can be formulated in terms of parabolic and hyperbolic partial differential equations (PDEs) [1, 2]. Since analyzing these models and reactions has considerable importance in applied science and engineering, we should investigate the solutions of the associated governed PDEs by modern tools [3, 4]. Analytical approaches for solving PDEs have received the researchers attention for solving a large number of PDEs, but numerical methods are more favorable with respect to these approaches. One of the basic advantages of the numerical methods, which made them more popular with respect to the analytical approaches, is based on the implementation of the operational matrices of differentiation and high accurate Gauss quadrature rules instead of direct differentiation and integration. By using these operational matrices and quadratures, computations time in numerical algorithms may be decreased significantly. Therefore, in most of the numerical methods, one of the above-mentioned tools may be applied for speeding up the operations in solving PDEs [5].

In a numerical point of view, the methods that are based on the operational matrices of differentiation may be divided into the collocation and Tau matrix methods. In Tau matrix methods, all the known and unknown functions should be approximated in terms of a specific complete basis (such as orthogonal polynomials or trigonometric functions). Since the mentioned basis is complete, one can factorize this basis and simplify the considered PDEs into the associated system of algebraic equations. On the other hand, in every considered PDE one may have some boundary conditions and we should impose these conditions in the procedure of the Tau matrix method. Our numerical implementation experiences show that if the boundary conditions are first imposed to the main equations [5] (via transforming the basic PDEs into the associated integro-PDEs), then we may have more stable numerical solutions. However, in some research works (e.g., Taylor matrix method [6, 7], Euler matrix method [8], and Chebyshev Tau method [9]) boundary conditions are imposed after discretizing the main equations. Such ideas may have similar numerical results with [5] but produce more algebraic equations and this may take more computations time and therefore make these methods deficient. Another disadvantage of Tau matrix method is that these methods apply the operational matrices of product for solving PDEs with variable coefficients [10]. Such idea again reduces the accuracy of the Tau matrix method. Moreover, Tau matrix schemes in general have difficulties in the treatment of nonlinear multidimensional PDEs, since no straightforward algebraic system of equations can be efficiently achieved for such problems [11]. However, the collocation schemes for PDEs are rarely ready to implement and they can overcome the aforementioned challenges.

In recent years, some research works have been focused on the application of collocation methods for solving linear one-dimensional parabolic and hyperbolic (specially Telegraph equations) PDEs such as Chebyshev collocation method [12] and Bessel collocation method [13, 14]. But application of this scheme for solving multidimensional PDEs has had few results. This partially motivates us to propose a new collocation scheme, which is based on the Bernoulli polynomials, for solving linear multidimensional diffusion and wave PDEs. It should be noted that the proposed collocation method can be generalized for solving other linear multidimensional parabolic and hyperbolic PDEs. In this paper we consider the following linear one-dimensional diffusion equationwith the initial conditiontogether with the Dirichlet boundary conditionsAlso, we will consider the following linear one-dimensional wave equationwith the initial conditionstogether with the Dirichlet boundary conditions

It should be noted that all the known functions in above-mentioned PDEs such as , , , , and are some smooth functions. Numerical solutions of the above-mentioned equations are considered in the literature through some classical and modern techniques such as finite difference methods (FDMs) [1518], finite element methods (FEMs) [19], finite volume methods (FVMs) [20], dual reciprocity boundary integral equation method [21], interpolating scaling function method [22], and Haar wavelet scheme [23]. In this paper, we treat the above equations (together with their 2D generalization forms) by a numerical method which is based on Bernoulli polynomials as the basis and the well-known Chebyshev Gauss Lobatto (CGL) collocation points. In this regard, the considered equations are collocated and then transformed into the associated systems of linear algebraic equations which can be solved through some iterative methods such as GMRES.

The structure of the remainder of this paper is as follows. In the next section, we will review some preliminaries regarding the shifted Bernoulli polynomials and shifted operational matrices which play important roles throughout the paper. In Section 3, implementation of the Bernoulli collocation method (BCM) for solving both of the 1D equations (1)–(3) and (4)–(6) will be provided and in Section 4, generalization of BCM for solving 2D diffusion and wave equations will be implemented. In Section 5, several numerical examples are given to illustrate the spectral accuracy of the proposed method. Finally, some conclusions regarding the suggested BCM are conveyed in Section 6.

2. Preliminaries

Bernoulli polynomials form a complete basis in the interval . Since in our considered PDEs we have arbitrary intervals, we should introduce the shifted Bernoulli polynomials, which keep the property of completeness, and the shifted operational matrices of differentiation. The shifted Bernoulli polynomials can be constructed via the relation for any arbitrary positive integer index in the interval . It should be recalled that the standard Bernoulli polynomials, in the interval , have the following properties [24]:

Therefore the shifted Bernoulli polynomials have the following properties in the interval :

Now we assume that . The shifted operational matrix of differentiation associated with the shifted Bernoulli polynomials can be derived from the first relation in (8):where is the shifted operational matrix of differentiation in the interval . Because of the existence of two-variable functions in (1)–(3) and (4)–(6), we must extend the mentioned matrix through the Kronecker multiplication. For this purpose, assume thatwhere for all and for all . Evidently, , where .

In the following lemma, operational matrices of differentiation associated with the two-variable functions will be introduced.

Lemma 1. Suppose that and , where denotes the identity matrix of dimension and . One can conclude that(i)(ii)

Proof. See [24].

Since in 2D diffusion and wave equations, which will be considered in Section 4, we have three-variable functions, again one can develop the subject of operational matrices via the Tensor product. For this goal, we should suppose that where . Similar to the previous lemma, operational matrices of differentiation corresponding to the three-variable functions will be introduced.

Lemma 2. Assume that , , and , where . The following relations hold.  (i)(ii)(iii)

Proof. The proof of this lemma is similar to the previous lemma.

3. Bernoulli Collocation Method for 1D Equations

This section is devoted to implement BCM for solving 1D diffusion equations (1)–(3) and also 1D wave equations (4)–(6). In this regard, we first discretize 1D diffusion equation and then localize 1D wave equation in two separate subsections via a geometrical point of view. In both cases, we should assemble the associated algebraic equations to form a final linear system which can be solved by some efficient iterative solvers such as GMRES method. It should be noted that BCM has received considerable attentions during recent years. The interested readers can refer to [25] and the references therein for application of BCM in solving complex ODEs, nonlinear two-point BVPs, delay Pantograph ODEs, and other types of applied mathematics problems.

3.1. 1D Diffusion Equations

In this subsection, the aim is to reduce (1)–(3) into the associated system of linear algebraic equations in the form of by implementing the BCM. This process has three steps.(1)Collocating the main equation (1).(2)Collocating the boundary conditions (3).(3)Collocating the initial condition (2).

At first we assume thatwhere is defined in the previous section and

It should be noted that our aim is to find the unknown vector . From the previous section we have

Also, we should define some suitable collocation points. In this paper, we choose the CGL points in the following form:

Figure 1 plots such collocation points for in the unit square (i.e., ) for the diffusion equation. At the first step, we collocate the main equation via the interior collocation points.

Step 1 (interior collocation points: main equation (1)). In this part of the paper, we have and (see the purple rhombics in Figure 1). So we should define with rows and columns. Moreover, we should define a column vector with elements. Therefore, one can refer to the first step of Algorithm 1 for this purpose.

Step : Collocating the main equation (1):
for
for
;
;
end
end
Step : Collocating the boundary conditions (3):
for
;
;
;
;
end
Step : Collocating the initial condition (2):
for
;
end
Step : Assembling the coefficient matrices and right hand side column vectors:
;
;
maxit = ;
tol = ;
= gmres (, , tol, maxit)

Step 2 (boundary conditions collocation (3)). In this part of the paper, we have for the left boundary condition and for the right boundary condition. In both cases, we have . Therefore, we should define (left boundary condition: see the red triangles in Figure 1) and (right boundary condition: see the blue down triangles in Figure 1) with rows and columns. Moreover, we should define column vectors (left boundary condition) and (right boundary condition) with elements. Thus, one can refer to the second step of Algorithm 1 for this goal.

Step 3 (initial condition collocation (2)). In this part of the paper, we have and also (see the black squares in Figure 1). In this case we should define matrix with rows and columns. Moreover, we should define column vector with elements. Therefore, one can refer to the third step of Algorithm 1 for this aim.

Step 4 (assembling the matrices and column vectors). Finally, by assembling the coefficient matrices and right hand side column vectors, we will reach to the system of linear algebraic equations in the form of , which can be solved by some appropriate iterative methods such as GMRES algorithm. Algorithm 1 can describe BCM for solving 1D diffusion equations (1)–(3).

3.2. 1D Wave Equations

This subsection is similar to the previous subsection and we have just an extra initial condition . Also, the aim is to reduce (4)–(6) into the associated system of linear equations in the form of by implementing the BCM. Similar to the previous subsection, this process has three steps.(1)Collocating the main equation (4).(2)Collocating the boundary conditions (6).(3)Collocating the initial conditions (5).

Again, we assume a similar approximate solution in the form of (17) and consider the collocation points in the rectangle that was previously introduced in (20). Similar to the previous subsection, one can write

It should be noted that Figure 2 plots such collocation points for in the unit square for the wave equation. The readers can see the difference between this figure and the previous figure in distribution of collocation points. At the first step, we collocate the main equation via the interior collocation points.

Step 1 (interior collocation points: main equation (4)). In this part of the paper, we have and (see the purple rhombics in Figure 2). So we should define with rows and columns. Moreover, we should define a column vector with elements. Therefore, one can refer to the first step of Algorithm 2 for this purpose.

Step : Collocating the main equation (4):
for
for
;
;
end
end
Step : Collocating the boundary conditions (6):
for
;
;
;
;
End
Step : Collocating the initial conditions (5):
for
;
;
;
;
end
Step : Assembling the coefficient matrices and right hand side column vectors:
;
;
maxit = ;
tol = ;
= gmres(, , tol, maxit)

Step 2 (boundary conditions collocation (6)). In this stage, we have for the left boundary condition and for the right boundary condition. In both cases, we have . Therefore, we should define (left boundary condition: see the red triangles in Figure 2) and (right boundary condition: see the blue down triangles in Figure 2) with rows and columns. Moreover, we should define column vectors (left boundary condition) and (right boundary condition) with elements. Therefore, one can refer to the second step of Algorithm 2 for this goal.

Step 3 (initial conditions collocation (5)). For the initial condition , we have and also (see the black squares in Figure 2) and for the initial condition we have and also (see the green circles in Figure 2). In these cases, we should define and matrices with rows and columns. Moreover, we should define and column vectors with elements. Therefore, one can refer to the third step of Algorithm 2 for this aim.
In this case, assembling the coefficient matrices and right hand side column vectors is similar to the previous subsection. Algorithm 2 can describe BCM for solving 1D wave equations (4)–(6).

4. Bernoulli Collocation Method for 2D Equations

In this section, we generalize the BCM for solving linear 2D diffusion and wave equations. Similar to the previous section, we depict the collocation points associated with the initial conditions and boundary conditions in each case of the diffusion and wave equations for clarity of presentation. Therefore, we consider the following linear 2D diffusion equation:with the initial conditiontogether with the Dirichlet boundary conditions

Also, we will consider the following linear 2D wave equation:with the initial conditiontogether with the Dirichlet boundary conditions

In the next subsection, we will extend BCM for solving 2D diffusion equations (22)–(24).

4.1. BCM for Solving 2D Diffusion Equations

We suppose that the approximate solution of (22)–(24) can be written in the following form ofwhere and in which

According to discussions in the second section, one can deduce that

Similar to the 1D cases, one can consider the following CGL collocation points in 2D equations in the rectangle cube :

Figure 3 shows only the collocation points associated with the initial condition (23) and also the boundary conditions (24) for 2D diffusion equation for . It should be noted that depicting the collocation points associated with the main equation (22) may make the readers confused. Moreover, in this figure stands for the variable, stands for the variable, and stands for the variable.

BCM for solving (22)–(24) has also three steps: main equation collocation; initial condition collocation; and boundary conditions collocation. Despite the 1D cases, we first start with the initial conditions.

Step 1 (initial condition collocation (23)). For the initial condition we have and also , (see the black squares in Figure 3). In this case we should define matrix with rows and columns. Moreover, we should define column vector with elements. Therefore, one can refer to the first step of Algorithm 3 for this aim.

Step : Collocating the initial condition (23):
for
for
;
;
end
end
Step : Collocating the boundary conditions (24):
variable discretization:
for
for
;
;
;
;
end
variable discretization:
for
for
;
;
;
;
end
Step : Collocating the main equation (22):
for
for
for
;
;
end
end
end
Step : Assembling the coefficient matrices and right hand side column vectors:
;
;
maxit = ;
tol = ;
= gmres(, , tol, maxit)

Step 2 (boundary conditions collocation (24)). Since we have four boundary conditions, we should discretize them in two different ways. For the variable, we have for the left boundary condition and for the right boundary condition. In both cases, we have and . Therefore, we should define (left boundary condition: see the red circles in Figure 3) and (right boundary condition: see the green circles in Figure 3) with rows and columns. Moreover, we should define column vectors (left boundary condition) and (right boundary condition) with elements. For the variable, we have for the bottom boundary condition and for the top boundary condition. In both cases, we have and . Therefore, we should define (bottom boundary condition: see the blue circles in Figure 3) and (top boundary condition: see the yellow circles in Figure 3) with rows and columns. Moreover, we should define column vectors (bottom boundary condition) and (top boundary condition) with elements. Therefore, one can refer to the second step of Algorithm 3 for this goal.

Step 3 (interior collocation points: main equation (22)). In this part of the paper, we have , , and . So we should define with rows and columns. Moreover, we should define a column vector with elements. Therefore, one can refer to the third step of Algorithm 3 for this purpose.
In this case, assembling the coefficient matrices and right hand side column vectors is similar to the previous subsection. Algorithm 3 can describe BCM for solving 2D diffusion equations (22)–(24).

4.2. BCM for Solving 2D Wave Equations

Similar to the previous subsection, we should consider an approximate solution in the form of (28) and also the collocation points in the rectangle cube that is introduced in (31). Moreover, one can write

Distribution of collocation points in 2D wave equations is similar to 2D diffusion equations and just has a difference (see the brown squares in Figure 4). Such brown squares are related to the initial condition .

For clarity of presentation, we just provide Algorithm 4 for describing BCM for solving (25)–(27). In this algorithm, we first discretize initial conditions (26) and then localize boundary conditions (27) and finally collocate the main equation (25) via the shifted CGL collocation points.

Step : Collocating the initial condition (26):
for
for
;
;
;
;
end
end
Step : Collocating the boundary conditions (27):
variable discretization:
for
for
;
;
;
;
end
variable discretization:
for
for
;
;
;
;
end
Step : Collocating the main equation (25):
for
for
for
;
;
end
end
end
Step : Assembling the coefficient matrices and right hand side column vectors:
;
;
maxit = ;
tol = ;
= gmres(, , tol, maxit)

5. Numerical Examples

In this section, we will test our method for solving both of the considered diffusion and wave equations. In the first example, we consider a constructed 1D diffusion equation and in the second example, we provide a 1D wave which is selected from [6]. Moreover, in the third and fourth examples, we will provide two test problems that are selected from [5] in the case of 2D diffusion equations. The presented idea (i.e., BCM) is easy to implement and can be applied for solving other 2D parabolic and hyperbolic equations such as 2D Telegraph equations [26]. All of the programs associated with the implementation of BCM for solving the considered examples are written in Matlab 2015b in a Laptop PC with 12 GB Ram with Cash of 6. Readers of this paper can communicate for receiving the associated codes for each example via email of the corresponding author. It should be noted that the tolerance for the GMRES solver is set to be and also the iterations associated for applying this solver is set to be “size (A, 1).” GMRES algorithm for solving the linear algebraic systems of the considered problems has a good performance and no preconditioning technique is needed for solving the algebraic systems.

Example 1. As the first example we consider (1)–(3) with the following assumptions: with the exact solution in the unit square . We have implemented the BCM for solving this constructed example by assuming several values of starting from 4 up to 14. It is observed that as gets greater values, the associated errors decreased. For this example and the next example, we have defined the error function , where is the approximate solution that is computed via applying BCM for solving the considered problems. Figures 5 and 6 depict the associated with the obtained numerical solution of the this example for and , respectively. From these figures, one can see the robustness of the proposed method for solving 1D diffusion equations.

Example 2 (see [6]). As the second example we consider (4)–(6) with the following assumptions: with the exact solution in the rectangle . Again, we have implemented the BCM for solving this provided example by assuming several values of starting from 4 up to 14. Since this test problem has been selected from [6], we make some error comparisons with the Taylor matrix method (TMM) in Table 1 at some special points in the considered computational domain. From this table one can see that the errors associated with the BCM are more stable and accurate than TMM [6]. Since in [6] the authors considered the Maclaurin approximation for the computational interval , this assumption may decrease the accuracy of the TMM. Figures 7 and 8 illustrate associated with the obtained numerical solution of this example via BCM for and , respectively. From these figures, one can conclude that the proposed BCM is an efficient numerical approach for solving 1D wave equations. For instance, in Figure 8, the error function values are around approximately, while the associated error function values of the TMM [6] are around (see the fourth figure of [6] for this aim).

Example 3 (see [5]). As the third example we consider (22)–(24) with the following assumptions: with the exact solution in the unit cube . The numerical results of this example are provided with the numerical results of the next example in Table 2. Moreover, numerical solution , which is obtained via the proposed BCM, together with the exact solution is depicted in Figure 9.

Example 4 (see [5]). As the fourth example we consider (22)–(24) with the following assumptions: with the exact solution in the unit cube . For solving this example and the previous example, we applied several values of such as 4, 6, 8, 10, 12, and 14 and achieve the numerical solution by the suggested BCM. For these numerical solutions, we have defined the discrete error function for both of the third and fourth examples and make a comparison with the Bernoulli matrix method (BMM) [5] in Table 2. From this table, one can conclude that the proposed BCM is more accurate with respect to the numerical solutions achieved via BMM [5]. Moreover, numerical solution , which is obtained via the proposed BCM, together with the exact solution is depicted in Figure 10.

6. Conclusions

Applying analytical methods for solving multidimensional parabolic and hyperbolic PDEs usually is not efficient, because symbolic computations (such as direct differentiation and integration) in higher dimensions are time consuming. Therefore, numerical methods which are based on the operational matrices of differentiation (instead of direct differentiation) and high accurate Gauss quadrature rules (instead of direct integration) are more favorable for solving PDEs. In this paper, we apply a numerical approach, which is based on the Bernoulli polynomials and their operational matrices of differentiation and also the CGL collocation points for solving multidimensional diffusion and wave equations. In this regard, the above-mentioned equations will be reduced to the associated systems of linear algebraic equations which can be solved by robust iterative solvers such as GMRES algorithm. Moreover, geometrical distribution of collocation points for the considered PDEs is depicted for clarity of presentation. To show the robustness of the proposed method, some test problems are provided. In all of the considered numerical examples, as gets greater, more accurate solutions will be achieved. Also, more accurate results (even in larger computational domain) in the second example are obtained with respect to the Taylor matrix method (TMM) [6] which confirm high accuracy of the BCM. The presented method is easy to implement in any software such as Matlab and Maple and can be generalized for solving three-dimensional PDEs. In this case, BCM for collocating the space variables should be combined with method of lines (MOL) [26, 27] to reduce three-dimensional PDEs into the associated systems of ODEs that can be integrated by some efficient ODEs solvers.

Competing Interests

The authors declare that they do not have any conflict of interests in their submitted paper.