Bernoulli Collocation Method for Solving Linear Multidimensional Diffusion and Wave Equations with Dirichlet Boundary Conditions

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


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 2 Advances in Mathematical Physics 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 equation with the initial condition together with the Dirichlet boundary conditions Also, we will consider the following linear one-dimensional wave equation with the initial conditions together with the Dirichlet boundary conditions It should be noted that all the known functions in abovementioned PDEs such as (), (),  1 (),  2 (), 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) [15][16][17][18], 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.
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 () = [ 0 ()  1 () ⋅ ⋅ ⋅   ()] 1×(+1) .Similar to the previous lemma, operational matrices of differentiation corresponding to the three-variable functions will be introduced.
Lemma 2. Assume that M =  (+1) ( Proof.The proof of this lemma is similar to the previous lemma.

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 AU =  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.

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 AU =  by implementing the BCM.This process has three steps.
At first we assume that where (, ) 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  = 8 in the unit square (i.e., [  ,   ] × [  ,   ] = [0, 1] × [0, 1]) 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 2 ≤  ≤  and 2 ≤  ≤  + 1 (see the purple rhombics in Figure 1).So we should define  main with ( − 1) rows and ( + 1) 2 columns.Moreover, we should define a column vector  main with ( − 1) elements.Therefore, one can refer to the first step of Algorithm 1 for this purpose.
Step 2 (boundary conditions collocation (3)).In this part of the paper, we have  = 1 for the left boundary condition and  =  + 1 for the right boundary condition.In both cases, we have 2 ≤  ≤  + 1.Therefore, we should define  bcl (left boundary condition: see the red triangles in Figure 1) and  bcr (right boundary condition: see the blue down triangles in Figure 1) with  rows and ( + 1) 2 columns.Moreover, we should define column vectors  bcl (left boundary condition) and  bcr (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  = 1 and also 1 ≤  ≤  + 1 (see the black squares in Figure 1).In this case we should define  in matrix with (+1) rows and (+1) 2 columns.Moreover, we should define  in column vector with (+1) 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 AU = , which can be solved by some appropriate iterative methods such as GMRES algorithm.

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 AU =  by implementing the BCM.Similar to the previous subsection, this process has three steps.
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  = 8 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 2 ≤  ≤  and 3 ≤  ≤  + 1  (see the purple rhombics in Figure 2).So we should define  main with (−1) 2 rows and (+1) 2 columns.Moreover, we should define a column vector  main with ( − 1) 2 elements.Therefore, one can refer to the first step of Algorithm 2 for this purpose.
Step 2 (boundary conditions collocation ( 6)).In this stage, we have  = 1 for the left boundary condition and  =  + 1 Step 1: Collocating the main equation ( 4 for the right boundary condition.In both cases, we have 3 ≤  ≤  + 1.Therefore, we should define  bcl (left boundary condition: see the red triangles in Figure 2) and  bcr (right boundary condition: see the blue down triangles in Figure 2) with (−1) rows and (+1) 2 columns.Moreover, we should define column vectors  bcl (left boundary condition) and  bcr (right boundary condition) with ( − 1) elements.Therefore, one can refer to the second step of Algorithm 2 for this goal.
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).

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 condition Advances in Mathematical Physics 7 Also, we will consider the following linear 2D wave equation: In the next subsection, we will extend BCM for solving 2D diffusion equations ( 22)-(24).

BCM for Solving 2D Diffusion
Equations.We suppose that the approximate solution of ( 22)-( 24 According to discussions in the second section, one can deduce that   (, , ) ≈  , (, , ) =  (, , ) M,   (, , ) ≈  , (, , ) =  (, , ) ( M) Figure 3 shows only the collocation points associated with the initial condition (23) and also the boundary conditions (24) for 2D diffusion equation for  = 6.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 2 (boundary conditions collocation (24)).Since we have four boundary conditions, we should discretize them in two different ways.For the  variable, we have  = 1 for the left boundary condition and  =  + 1 for the right boundary condition.In both cases, we have 1 ≤  ≤  + 1 and 2 ≤  ≤  + 1.Therefore, we should define  bcl (left boundary condition: see the red circles in Figure 3) and  bcr (right boundary condition: see the green circles in Figure 3)  Step 3 (interior collocation points: main equation ( 22)).In this part of the paper, we have 2 ≤  ≤ , 2 ≤  ≤ , and 2 ≤  ≤  + 1.So we should define  main with ( − 1) 2 rows and ( + 1) 3 columns.Moreover, we should define a column vector  main with ( − 1) 2 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).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 (, ,   )/ = (, ).

BCM for
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.

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 10 −(−2) 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.
Step 1: Collocating the initial condition ( 23 Example 1.As the first example we consider (1)-( 3) with the following assumptions: with the exact solution (, ) = exp( + ) 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  = 7, 9 and  = 11, 13, 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:  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 [  ,   ] = [0, ], 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  = 6,8 and  = 10, 12, 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  12 (, ) are around 10 −11 approximately, while the associated error function values  12 (, ) of the TMM [6] are around 10 −7 (see the fourth figure of [6] for this aim).
Example 3 (see [5]).As the third example we consider ( 22)-  with the exact solution (, , ) = exp( The numerical results of this example are provided with the numerical results of the next example in Table 2.Moreover, numerical solution  7 (1, , ), which is obtained via the proposed BCM, together with the exact solution (1, , ) is depicted in Figure 9.
Example 4 (see [5]).As the fourth example we consider ( 22)-( 24) with the following assumptions:   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 (  ,   , 1) −   (  ,   , 1)      (37) 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  9 (1, , ), which is obtained via the proposed BCM, together with the exact solution (1, , ) is depicted in Figure 10.

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,

Figure 1 :
Figure 1: CGL collocation points in the unit square for  = 8 for 1D diffusion equation.

Figure 2 :
Figure 2: CGL collocation points in the unit square for  = 8 for 1D wave equation.

Figure 3 :
Figure 3: CGL collocation points in the unit cube for  = 6 for 2D diffusion equation.

Figure 5 :Figure 6 :
Figure 5: Error histories associated with the obtained numerical solution of the first example for  = 7, 9.

Figure 7 :Figure 8 :
Figure 7: Error histories associated with the obtained numerical solution of the second example for  = 6, 8.

Table 1 :
Numerical results of the BCM and TMM associated with the error   (, ) in the second example.with the exact solution (, ) = sin() cos() in the rectangle [  ,   ] × [  ,   ] = [0, ] × [0, 1].Again, we have implemented the BCM for solving this provided example by assuming several values of  starting from 4 up to 14.

Table 2 :
Comparisons of the BCM and BMM associated with the errors   in the third and fourth example.