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

Authors: Bashar Zogheib; Emran Tohidi; Stanford Shateyi
Journal: Advances in Mathematical Physics
 (2017)

---

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


---

## Body

## 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 equation(1)∂ψ∂τ=∂2ψ∂η2+ωη,τ,η,τ∈aη,bη×aτ,bτ,with the initial condition(2)ψη,aτ=fη,η∈aη,bη,together with the Dirichlet boundary conditions(3)ψaη,τ=g1τ,ψbη,τ=g2τ,τ∈aτ,bτ.Also, we will consider the following linear one-dimensional wave equation(4)∂2ψ∂τ2=∂2ψ∂η2+ωη,τ,η,τ∈aη,bη×aτ,bτ,with the initial conditions(5)ψη,aτ=fη,∂ψ∂τη,aτ=kη,η∈aη,bη,together with the Dirichlet boundary conditions(6)ψaη,τ=g1τ,ψbη,τ=g2τ,τ∈aτ,bτ.It should be noted that all the known functions in above-mentioned PDEs such asfη, kη, g1τ, g2τ, 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–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 Section3, 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 interval0,1. 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 βiη=Biη-aη/bη-aη for any arbitrary positive integer index i in the interval aη,bη. It should be recalled that the standard Bernoulli polynomials, in the interval 0,1, have the following properties [24]: (7)Bi′x=iBi-1x,∀i≥1,∫01Bixdx=0,∀i≥1,B0x=1.Therefore the shifted Bernoulli polynomials have the following properties in the intervalaη,bη:(8)βi′η=ibη-aηβi-1η,∀i≥1,∫aηbηβiηdη=0,∀i≥1,β0η=1.Now we assume thatβη=β0ηβ1η⋯βNη1×N+1. The shifted operational matrix of differentiation associated with the shifted Bernoulli polynomials can be derived from the first relation in (8):(9)β0ηβ1η⋮βN-1ηβNη′=1bη-aη000⋯0100⋯0020⋯0⋮⋮⋱⋱⋮00⋯N0︷Mηβ0ηβ1η⋮βN-1ηβNη,where Mη is the shifted operational matrix of differentiation in the interval aη,bη. 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 that(10)βη,τ=β0η,τβ1η,τ⋯βNη,τ1×N+12,where βiη,τ=β0η,τβ1η,τ⋯βNη,τ for all i=0,1,…,N and βm,nη,τ=βmηβnτ for all m,n=0,1,…,N. Evidently, βη,τ=βη⊗βτ, where βτ=β0τβ1τ⋯βNτ1×N+1.In the following lemma, operational matrices of differentiation associated with the two-variable functions will be introduced.Lemma 1.
Suppose thatM^=IN+1⊗MτT and M~=MηT⊗IN+1, where IN+1 denotes the identity matrix of dimension N+1 and Mτ=bη-aη/bτ-aτMη. One can conclude that(i)
⁡(11)∂βη,τ∂τ=βτη,τ=βη,τM^,∂2βη,τ∂τ2=βττη,τ=βη,τM^2,(ii)
⁡(12)∂βη,τ∂η=βηη,τ=βη,τM~,∂2βη,τ∂η2=βηηη,τ=βη,τM~2.Proof.
See [24].Since in 2D diffusion and wave equations, which will be considered in Section4, 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 (13)βη,ξ,τ=βη⊗βξ⊗βτ1×N+13,where βξ=β0ξβ1ξ⋯βNξ1×N+1. Similar to the previous lemma, operational matrices of differentiation corresponding to the three-variable functions will be introduced.Lemma 2.
Assume thatM^^=IN+12⊗MτT, M~~=MηT⊗IN+12, and M˙˙=IN+1⊗MξT⊗IN+1, where Mξ=bη-aη/bξ-aξMη. The following relations hold. (i)
⁡(14)∂βη,ξ,τ∂τ=βτη,ξ,τ=βη,ξ,τM^^,∂2βη,ξ,τ∂τ2=βττη,ξ,τ=βη,ξ,τM^^2.(ii)
⁡(15)∂βη,ξ,τ∂η=βηη,ξ,τ=βη,ξ,τM~~,∂2βη,ξ,τ∂η2=βηηη,ξ,τ=βη,ξ,τM~~2.(iii)
⁡(16)∂βη,ξ,τ∂ξ=βξη,ξ,τ=βη,ξ,τM˙˙,∂2βη,ξ,τ∂ξ2=βξξη,ξ,τ=βη,ξ,τM˙˙2.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 AU=b 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 AU=b 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 that(17)ψη,τ≈ψNη,τ=∑m=0N∑n=0Numnβmηβnτ=βη,τU,where βη,τ is defined in the previous section and (18)U=u00u01⋯u0Nu10u11⋯u1N⋯uN0uN1⋯uNNN+12×1T.It should be noted that our aim is to find the unknown vectorU. From the previous section we have(19)ψτη,τ≈ψN,τη,τ=βη,τM^U,ψηηη,τ≈ψN,ηηη,τ=βη,τM~2U.Also, we should define some suitable collocation points. In this paper, we choose the CGL points in the following form:(20)ηi=bη-aη21-cosi-1πN+aη,i=1,2,…,N+1,τj=bτ-aτ21-cosj-1πN+aτ,j=1,2,…,N+1.Figure1 plots such collocation points for N=8 in the unit square (i.e., aη,bη×aτ,bτ=0,1×0,1) for the diffusion equation. At the first step, we collocate the main equation via the interior collocation points.Figure 1
CGL collocation points in the unit square forN=8 for 1D diffusion equation.Step 1 (interior collocation points: main equation (1)).
In this part of the paper, we have2≤i≤N and 2≤j≤N+1 (see the purple rhombics in Figure 1). So we should define Amain with NN-1 rows and N+12 columns. Moreover, we should define a column vector bmain with NN-1 elements. Therefore, one can refer to the first step of Algorithm 1 for this purpose.Algorithm 1:Implementation of BCM for solving (1)–(3).
Step1: Collocating the main equation (1): fori=2:N forj=2:N+1 A m a i n i - 2 ∗ N + j - 1 , : = β η i , τ j ∗ M ^ - M ~ 2; b m a i n i - 2 ∗ N + j - 1 , 1 = ω η i , τ j; end end Step2: Collocating the boundary conditions (3): forj=2:N+1 A b c l j - 1 , : = β a η , τ j; A b c r j - 1 , : = β b η , τ j; b b c l j - 1,1 = g 1 τ j; b b c r j - 1,1 = g 2 τ j; end Step3: Collocating the initial condition (2): fori=1:N+1 A i n i , : = β η i , a τ; b i n i , 1 = f η i end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A i n; b = b m a i n ; b b c l ; b b c r ; b i n; maxit =sizeA,1; tol =10-N-2; U = gmres (A,b, [], tol, maxit)Step 2 (boundary conditions collocation (3)).
In this part of the paper, we havei=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 2≤j≤N+1. Therefore, we should define Abcl (left boundary condition: see the red triangles in Figure 1) and Abcr (right boundary condition: see the blue down triangles in Figure 1) with N rows and N+12 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with N 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 havej=1 and also 1≤i≤N+1 (see the black squares in Figure 1). In this case we should define Ain matrix with N+1 rows and N+12 columns. Moreover, we should define bin column vector with N+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 ofAU=b, 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∂ψ/∂τη,aτ=kη. Also, the aim is to reduce (4)–(6) into the associated system of linear equations in the form of AU=b 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 aη,bη×aτ,bτ that was previously introduced in (20). Similar to the previous subsection, one can write (21)ψττη,τ≈ψN,ττη,τ=βη,τM^2U.It should be noted that Figure2 plots such collocation points for N=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.Figure 2
CGL collocation points in the unit square forN=8 for 1D wave equation.Step 1 (interior collocation points: main equation (4)).
In this part of the paper, we have2≤i≤N and 3≤j≤N+1 (see the purple rhombics in Figure 2). So we should define Amain with N-12 rows and N+12 columns. Moreover, we should define a column vector bmain with N-12 elements. Therefore, one can refer to the first step of Algorithm 2 for this purpose.Algorithm 2:Implementation of BCM for solving (4)–(6).
Step1: Collocating the main equation (4): fori=2:N forj=3:N+1 A m a i n i - 2 ∗ N + j - 2 , : = β η i , τ j ∗ M ^ 2 - M ~ 2; b m a i n i - 2 ∗ N + j - 2 , 1 = ω η i , τ j; end end Step2: Collocating the boundary conditions (6): forj=3:N+1 A b c l j - 2 , : = β a η , τ j; A b c r j - 2 , : = β b η , τ j; b b c l j - 2,1 = g 1 τ j; b b c r j - 2,1 = g 2 τ j; End Step3: Collocating the initial conditions (5): fori=1:N+1 A i n i , : = β η i , a τ; A i n 1 i , : = β η i , a τ M ^; b i n i , 1 = f η i; b i n 1 i , 1 = k η i; end Step4: Assembling the coefficient matrices and right hand side column vectors: A = [ A m a i n ; A b c l ; A b c r ; A i n ; A i n 1 ]; b = b m a i n ; b b c l ; b b c r ; b i n ; b i n 1; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], tol, maxit)Step 2 (boundary conditions collocation (6)).
In this stage, we havei=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 3≤j≤N+1. Therefore, we should define Abcl (left boundary condition: see the red triangles in Figure 2) and Abcr (right boundary condition: see the blue down triangles in Figure 2) with N-1 rows and N+12 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with N-1 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ψη,aτ=fη, we have j=1 and also 1≤i≤N+1 (see the black squares in Figure 2) and for the initial condition ∂ψ/∂τη,aτ=kη we have j=2 and also 1≤i≤N+1 (see the green circles in Figure 2). In these cases, we should define Ain and Ain1 matrices with N+1 rows and N+12 columns. Moreover, we should define bin and bin1 column vectors with N+1 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. Algorithm2 can describe BCM for solving 1D wave equations (4)–(6).


## 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 AU=b 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 that(17)ψη,τ≈ψNη,τ=∑m=0N∑n=0Numnβmηβnτ=βη,τU,where βη,τ is defined in the previous section and (18)U=u00u01⋯u0Nu10u11⋯u1N⋯uN0uN1⋯uNNN+12×1T.It should be noted that our aim is to find the unknown vectorU. From the previous section we have(19)ψτη,τ≈ψN,τη,τ=βη,τM^U,ψηηη,τ≈ψN,ηηη,τ=βη,τM~2U.Also, we should define some suitable collocation points. In this paper, we choose the CGL points in the following form:(20)ηi=bη-aη21-cosi-1πN+aη,i=1,2,…,N+1,τj=bτ-aτ21-cosj-1πN+aτ,j=1,2,…,N+1.Figure1 plots such collocation points for N=8 in the unit square (i.e., aη,bη×aτ,bτ=0,1×0,1) for the diffusion equation. At the first step, we collocate the main equation via the interior collocation points.Figure 1
CGL collocation points in the unit square forN=8 for 1D diffusion equation.Step 1 (interior collocation points: main equation (1)).
In this part of the paper, we have2≤i≤N and 2≤j≤N+1 (see the purple rhombics in Figure 1). So we should define Amain with NN-1 rows and N+12 columns. Moreover, we should define a column vector bmain with NN-1 elements. Therefore, one can refer to the first step of Algorithm 1 for this purpose.Algorithm 1:Implementation of BCM for solving (1)–(3).
Step1: Collocating the main equation (1): fori=2:N forj=2:N+1 A m a i n i - 2 ∗ N + j - 1 , : = β η i , τ j ∗ M ^ - M ~ 2; b m a i n i - 2 ∗ N + j - 1 , 1 = ω η i , τ j; end end Step2: Collocating the boundary conditions (3): forj=2:N+1 A b c l j - 1 , : = β a η , τ j; A b c r j - 1 , : = β b η , τ j; b b c l j - 1,1 = g 1 τ j; b b c r j - 1,1 = g 2 τ j; end Step3: Collocating the initial condition (2): fori=1:N+1 A i n i , : = β η i , a τ; b i n i , 1 = f η i end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A i n; b = b m a i n ; b b c l ; b b c r ; b i n; maxit =sizeA,1; tol =10-N-2; U = gmres (A,b, [], tol, maxit)Step 2 (boundary conditions collocation (3)).
In this part of the paper, we havei=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 2≤j≤N+1. Therefore, we should define Abcl (left boundary condition: see the red triangles in Figure 1) and Abcr (right boundary condition: see the blue down triangles in Figure 1) with N rows and N+12 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with N 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 havej=1 and also 1≤i≤N+1 (see the black squares in Figure 1). In this case we should define Ain matrix with N+1 rows and N+12 columns. Moreover, we should define bin column vector with N+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 ofAU=b, 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∂ψ/∂τη,aτ=kη. Also, the aim is to reduce (4)–(6) into the associated system of linear equations in the form of AU=b 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 aη,bη×aτ,bτ that was previously introduced in (20). Similar to the previous subsection, one can write (21)ψττη,τ≈ψN,ττη,τ=βη,τM^2U.It should be noted that Figure2 plots such collocation points for N=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.Figure 2
CGL collocation points in the unit square forN=8 for 1D wave equation.Step 1 (interior collocation points: main equation (4)).
In this part of the paper, we have2≤i≤N and 3≤j≤N+1 (see the purple rhombics in Figure 2). So we should define Amain with N-12 rows and N+12 columns. Moreover, we should define a column vector bmain with N-12 elements. Therefore, one can refer to the first step of Algorithm 2 for this purpose.Algorithm 2:Implementation of BCM for solving (4)–(6).
Step1: Collocating the main equation (4): fori=2:N forj=3:N+1 A m a i n i - 2 ∗ N + j - 2 , : = β η i , τ j ∗ M ^ 2 - M ~ 2; b m a i n i - 2 ∗ N + j - 2 , 1 = ω η i , τ j; end end Step2: Collocating the boundary conditions (6): forj=3:N+1 A b c l j - 2 , : = β a η , τ j; A b c r j - 2 , : = β b η , τ j; b b c l j - 2,1 = g 1 τ j; b b c r j - 2,1 = g 2 τ j; End Step3: Collocating the initial conditions (5): fori=1:N+1 A i n i , : = β η i , a τ; A i n 1 i , : = β η i , a τ M ^; b i n i , 1 = f η i; b i n 1 i , 1 = k η i; end Step4: Assembling the coefficient matrices and right hand side column vectors: A = [ A m a i n ; A b c l ; A b c r ; A i n ; A i n 1 ]; b = b m a i n ; b b c l ; b b c r ; b i n ; b i n 1; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], tol, maxit)Step 2 (boundary conditions collocation (6)).
In this stage, we havei=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 3≤j≤N+1. Therefore, we should define Abcl (left boundary condition: see the red triangles in Figure 2) and Abcr (right boundary condition: see the blue down triangles in Figure 2) with N-1 rows and N+12 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with N-1 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ψη,aτ=fη, we have j=1 and also 1≤i≤N+1 (see the black squares in Figure 2) and for the initial condition ∂ψ/∂τη,aτ=kη we have j=2 and also 1≤i≤N+1 (see the green circles in Figure 2). In these cases, we should define Ain and Ain1 matrices with N+1 rows and N+12 columns. Moreover, we should define bin and bin1 column vectors with N+1 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. Algorithm2 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:(22)∂ψ∂τ=∂2ψ∂η2+∂2ψ∂ξ2+ωη,ξ,τ,η,ξ,τ∈aη,bη×aξ,bξ×aτ,bτ,with the initial condition(23)ψη,ξ,aτ=fη,ξ,η,ξ∈aη,bη×aξ,bξ,together with the Dirichlet boundary conditions(24)ψaη,ξ,τ=g1ξ,τ,ψbη,ξ,τ=g2ξ,τ,ξ,τ∈aξ,bξ×aτ,bτ,ψη,aξ,τ=h1η,τ,ψη,bξ,τ=h2η,τ,η,τ∈aη,bη×aτ,bτ.Also, we will consider the following linear 2D wave equation:(25)∂2ψ∂τ2=∂2ψ∂η2+∂2ψ∂ξ2+ωη,ξ,τ,η,ξ,τ∈aη,bη×aξ,bξ×aτ,bτ,with the initial condition(26)ψη,ξ,aτ=fη,ξ,∂ψη,ξ,aτ∂τ=kη,ξ,η,ξ∈aη,bη×aξ,bξ,together with the Dirichlet boundary conditions(27)ψaη,ξ,τ=g1ξ,τ,ψbη,ξ,τ=g2ξ,τ,ξ,τ∈aξ,bξ×aτ,bτ,ψη,aξ,τ=h1η,τ,ψη,bξ,τ=h2η,τ,η,τ∈aη,bη×aτ,bτ.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 of(28)ψη,ξ,τ≈ψNη,ξ,τ=∑m=0N∑n=0N∑p=0Numnpβmηβnξβpτ=βη,ξ,τU,where βη,ξ,τ=βη⊗βξ⊗βτ and U=u0..Tu1..T⋯uN..TT in which (29)ui..=ui00ui01⋯ui0Nui10ui11⋯ui1N⋯uiN0uiN1⋯uiNN.According to discussions in the second section, one can deduce that(30)ψτη,ξ,τ≈ψN,τη,ξ,τ=βη,ξ,τM^^U,ψηηη,ξ,τ≈ψN,ηηη,ξ,τ=βη,ξ,τM~~2U,ψξξη,ξ,τ≈ψN,ξξη,ξ,τ=βη,ξ,τM˙˙2U.Similar to the 1D cases, one can consider the following CGL collocation points in 2D equations in the rectangle cubeaη,bη×aξ,bξ×aτ,bτ:(31)ηi=bη-aη21-cos⁡i-1πN+aη,i=1,2,…,N,ξj=bξ-aξ21-cos⁡j-1πN+aξ,j=1,2,…,N,τk=bτ-aτ21-cos⁡k-1πN+aτ,k=1,2,…,N.Figure3 shows only the collocation points associated with the initial condition (23) and also the boundary conditions (24) for 2D diffusion equation for N=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 x stands for the η variable, y stands for the ξ variable, and z stands for the τ variable.Figure 3
CGL collocation points in the unit cube forN=6 for 2D diffusion equation.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ψη,ξ,aτ=fη,ξ we have k=1 and also 1≤i≤N+1, 1≤j≤N+1 (see the black squares in Figure 3). In this case we should define Ain matrix with N+12 rows and N+13 columns. Moreover, we should define bin column vector with N+12 elements. Therefore, one can refer to the first step of Algorithm 3 for this aim.Algorithm 3:Implementation of BCM for solving (22)–(24).
Step1: Collocating the initial condition (23): fori=1:N+1 forj=1:N+1 A i n i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ; b i n i - 1 ∗ N + 1 + j , 1 = f η i , ξ j; end end Step2: Collocating the boundary conditions (24): → x variable discretization: forj=1:N+1 fork=2:N+1 A b c l j - 1 ∗ N + k - 1 , : = β a η , ξ j , τ k; A b c r j - 1 ∗ N + k - 1 , : = β b η , ξ j , τ k; b b c l j - 1 ∗ N + k - 1 , 1 = g 1 ξ j , τ k; b b c r j - 1 ∗ N + k - 1 , 1 = g 2 ξ j , τ k; end → y variable discretization: fori=2:N fork=2:N+1 A b c b i - 2 ∗ N + k - 1 , : = β η i , a ξ , τ k; A b c t i - 2 ∗ N + k - 1 , : = β η i , b ξ , τ k; b b c l i - 2 ∗ N + k - 1 , 1 = h 1 η i , τ k; b b c r i - 2 ∗ N + k - 1 , 1 = h 2 η i , τ k; end Step3: Collocating the main equation (22): fori=2:N forj=2:N fork=2:N+1 A m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 1 , : = β η i , ξ j , τ k ∗ M ^ ^ - M ~ ~ 2 - M ˙ ˙ 2; b m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 1 , 1 = ω η i , ξ j , τ k; end end end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A b c b ; A b c t ; A i n; b = b m a i n ; b b c l ; b b c r ; b b c b ; b b c t ; b i n; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], tol, maxit)Step 2 (boundary conditions collocation (24)).
Since we have four boundary conditions, we should discretize them in two different ways. For thex variable, we have i=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 1≤j≤N+1 and 2≤k≤N+1. Therefore, we should define Abcl (left boundary condition: see the red circles in Figure 3) and Abcr (right boundary condition: see the green circles in Figure 3) with NN+1 rows and N+13 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with NN+1 elements. For the y variable, we have j=1 for the bottom boundary condition and j=N+1 for the top boundary condition. In both cases, we have 2≤i≤N and 2≤k≤N+1. Therefore, we should define Abcb (bottom boundary condition: see the blue circles in Figure 3) and Abct (top boundary condition: see the yellow circles in Figure 3) with NN-1 rows and N+13 columns. Moreover, we should define column vectors bbcb (bottom boundary condition) and bbct (top boundary condition) with NN-1 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 have2≤i≤N, 2≤j≤N, and 2≤k≤N+1. So we should define Amain with NN-12 rows and N+13 columns. Moreover, we should define a column vector bmain with NN-12 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. Algorithm3 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 aη,bη×aξ,bξ×aτ,bτ that is introduced in (31). Moreover, one can write (32)ψττη,ξ,τ≈ψN,ττη,ξ,τ=βη,ξ,τM^^2U.Distribution of collocation points in 2D wave equations is similar to 2D diffusion equations and just has a difference (see the brown squares in Figure4). Such brown squares are related to the initial condition ∂ψη,ξ,aτ/∂τ=kη,ξ.Figure 4
CGL collocation points in the unit cube forN=6 for 2D wave equation.For clarity of presentation, we just provide Algorithm4 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.Algorithm 4:Implementation of BCM for solving (25)–(27).
Step1: Collocating the initial condition (26): fori=1:N+1 forj=1:N+1 A i n i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ; A i n 1 i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ ∗ M ^ ^; b i n i - 1 ∗ N + 1 + j , 1 = f η i , ξ j; b i n 1 i - 1 ∗ N + 1 + j , 1 = k η i , ξ j; end end Step2: Collocating the boundary conditions (27): → x variable discretization: forj=1:N+1 fork=3:N+1 A b c l j - 1 ∗ N + k - 2 , : = β a η , ξ j , τ k; A b c r j - 1 ∗ N + k - 2 , : = β b η , ξ j , τ k; b b c l j - 1 ∗ N + k - 2 , 1 = g 1 ξ j , τ k; b b c r j - 1 ∗ N + k - 2 , 1 = g 2 ξ j , τ k; end → y variable discretization: fori=2:N fork=3:N+1 A b c b i - 2 ∗ N + k - 2 , : = β η i , a ξ , τ k; A b c t i - 2 ∗ N + k - 2 , : = β η i , b ξ , τ k; b b c l i - 2 ∗ N + k - 2 , 1 = h 1 η i , τ k; b b c r i - 2 ∗ N + k - 2 , 1 = h 2 η i , τ k; end Step3: Collocating the main equation (25): fori=2:N forj=2:N fork=3:N+1 A m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 2 , : = β η i , ξ j , τ k ∗ M ^ ^ 2 - M ~ ~ 2 - M ˙ ˙ 2; b m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 2 , 1 = ω η i , ξ j , τ k; end end end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A b c b ; A b c t ; A i n ; A i n 1; b = b m a i n ; b b c l ; b b c r ; b b c b ; b b c t ; b i n ; b i n 1; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], tol, maxit)


## 4.1. BCM for Solving 2D Diffusion Equations

We suppose that the approximate solution of (22)–(24) can be written in the following form of(28)ψη,ξ,τ≈ψNη,ξ,τ=∑m=0N∑n=0N∑p=0Numnpβmηβnξβpτ=βη,ξ,τU,where βη,ξ,τ=βη⊗βξ⊗βτ and U=u0..Tu1..T⋯uN..TT in which (29)ui..=ui00ui01⋯ui0Nui10ui11⋯ui1N⋯uiN0uiN1⋯uiNN.According to discussions in the second section, one can deduce that(30)ψτη,ξ,τ≈ψN,τη,ξ,τ=βη,ξ,τM^^U,ψηηη,ξ,τ≈ψN,ηηη,ξ,τ=βη,ξ,τM~~2U,ψξξη,ξ,τ≈ψN,ξξη,ξ,τ=βη,ξ,τM˙˙2U.Similar to the 1D cases, one can consider the following CGL collocation points in 2D equations in the rectangle cubeaη,bη×aξ,bξ×aτ,bτ:(31)ηi=bη-aη21-cos⁡i-1πN+aη,i=1,2,…,N,ξj=bξ-aξ21-cos⁡j-1πN+aξ,j=1,2,…,N,τk=bτ-aτ21-cos⁡k-1πN+aτ,k=1,2,…,N.Figure3 shows only the collocation points associated with the initial condition (23) and also the boundary conditions (24) for 2D diffusion equation for N=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 x stands for the η variable, y stands for the ξ variable, and z stands for the τ variable.Figure 3
CGL collocation points in the unit cube forN=6 for 2D diffusion equation.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ψη,ξ,aτ=fη,ξ we have k=1 and also 1≤i≤N+1, 1≤j≤N+1 (see the black squares in Figure 3). In this case we should define Ain matrix with N+12 rows and N+13 columns. Moreover, we should define bin column vector with N+12 elements. Therefore, one can refer to the first step of Algorithm 3 for this aim.Algorithm 3:Implementation of BCM for solving (22)–(24).
Step1: Collocating the initial condition (23): fori=1:N+1 forj=1:N+1 A i n i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ; b i n i - 1 ∗ N + 1 + j , 1 = f η i , ξ j; end end Step2: Collocating the boundary conditions (24): → x variable discretization: forj=1:N+1 fork=2:N+1 A b c l j - 1 ∗ N + k - 1 , : = β a η , ξ j , τ k; A b c r j - 1 ∗ N + k - 1 , : = β b η , ξ j , τ k; b b c l j - 1 ∗ N + k - 1 , 1 = g 1 ξ j , τ k; b b c r j - 1 ∗ N + k - 1 , 1 = g 2 ξ j , τ k; end → y variable discretization: fori=2:N fork=2:N+1 A b c b i - 2 ∗ N + k - 1 , : = β η i , a ξ , τ k; A b c t i - 2 ∗ N + k - 1 , : = β η i , b ξ , τ k; b b c l i - 2 ∗ N + k - 1 , 1 = h 1 η i , τ k; b b c r i - 2 ∗ N + k - 1 , 1 = h 2 η i , τ k; end Step3: Collocating the main equation (22): fori=2:N forj=2:N fork=2:N+1 A m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 1 , : = β η i , ξ j , τ k ∗ M ^ ^ - M ~ ~ 2 - M ˙ ˙ 2; b m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 1 , 1 = ω η i , ξ j , τ k; end end end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A b c b ; A b c t ; A i n; b = b m a i n ; b b c l ; b b c r ; b b c b ; b b c t ; b i n; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], tol, maxit)Step 2 (boundary conditions collocation (24)).
Since we have four boundary conditions, we should discretize them in two different ways. For thex variable, we have i=1 for the left boundary condition and i=N+1 for the right boundary condition. In both cases, we have 1≤j≤N+1 and 2≤k≤N+1. Therefore, we should define Abcl (left boundary condition: see the red circles in Figure 3) and Abcr (right boundary condition: see the green circles in Figure 3) with NN+1 rows and N+13 columns. Moreover, we should define column vectors bbcl (left boundary condition) and bbcr (right boundary condition) with NN+1 elements. For the y variable, we have j=1 for the bottom boundary condition and j=N+1 for the top boundary condition. In both cases, we have 2≤i≤N and 2≤k≤N+1. Therefore, we should define Abcb (bottom boundary condition: see the blue circles in Figure 3) and Abct (top boundary condition: see the yellow circles in Figure 3) with NN-1 rows and N+13 columns. Moreover, we should define column vectors bbcb (bottom boundary condition) and bbct (top boundary condition) with NN-1 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 have2≤i≤N, 2≤j≤N, and 2≤k≤N+1. So we should define Amain with NN-12 rows and N+13 columns. Moreover, we should define a column vector bmain with NN-12 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. Algorithm3 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 aη,bη×aξ,bξ×aτ,bτ that is introduced in (31). Moreover, one can write (32)ψττη,ξ,τ≈ψN,ττη,ξ,τ=βη,ξ,τM^^2U.Distribution of collocation points in 2D wave equations is similar to 2D diffusion equations and just has a difference (see the brown squares in Figure4). Such brown squares are related to the initial condition ∂ψη,ξ,aτ/∂τ=kη,ξ.Figure 4
CGL collocation points in the unit cube forN=6 for 2D wave equation.For clarity of presentation, we just provide Algorithm4 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.Algorithm 4:Implementation of BCM for solving (25)–(27).
Step1: Collocating the initial condition (26): fori=1:N+1 forj=1:N+1 A i n i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ; A i n 1 i - 1 ∗ N + 1 + j , : = β η i , ξ j , a τ ∗ M ^ ^; b i n i - 1 ∗ N + 1 + j , 1 = f η i , ξ j; b i n 1 i - 1 ∗ N + 1 + j , 1 = k η i , ξ j; end end Step2: Collocating the boundary conditions (27): → x variable discretization: forj=1:N+1 fork=3:N+1 A b c l j - 1 ∗ N + k - 2 , : = β a η , ξ j , τ k; A b c r j - 1 ∗ N + k - 2 , : = β b η , ξ j , τ k; b b c l j - 1 ∗ N + k - 2 , 1 = g 1 ξ j , τ k; b b c r j - 1 ∗ N + k - 2 , 1 = g 2 ξ j , τ k; end → y variable discretization: fori=2:N fork=3:N+1 A b c b i - 2 ∗ N + k - 2 , : = β η i , a ξ , τ k; A b c t i - 2 ∗ N + k - 2 , : = β η i , b ξ , τ k; b b c l i - 2 ∗ N + k - 2 , 1 = h 1 η i , τ k; b b c r i - 2 ∗ N + k - 2 , 1 = h 2 η i , τ k; end Step3: Collocating the main equation (25): fori=2:N forj=2:N fork=3:N+1 A m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 2 , : = β η i , ξ j , τ k ∗ M ^ ^ 2 - M ~ ~ 2 - M ˙ ˙ 2; b m a i n i - 2 ∗ N ∗ N - 1 + j - 2 ∗ N + k - 2 , 1 = ω η i , ξ j , τ k; end end end Step4: Assembling the coefficient matrices and right hand side column vectors: A = A m a i n ; A b c l ; A b c r ; A b c b ; A b c t ; A i n ; A i n 1; b = b m a i n ; b b c l ; b b c r ; b b c b ; b b c t ; b i n ; b i n 1; maxit =sizeA,1; tol =10-N-2; U = gmres(A,b, [], 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 10-N-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.Example 1.
As the first example we consider (1)–(3) with the following assumptions: (33)ωη,τ=0;fη=expη;g1τ=expτ;g2τ=expτ+1,with the exact solution ψη,τ=expη+τ in the unit square aη,bη×aτ,bτ=0,1×0,1. We have implemented the BCM for solving this constructed example by assuming several values of N starting from 4 up to 14. It is observed that as N gets greater values, the associated errors decreased. For this example and the next example, we have defined the error function eNη,τ=ψη,τ-ψNη,τ, where ψNη,τ is the approximate solution that is computed via applying BCM for solving the considered problems. Figures 5 and 6 depict the eNη,τ associated with the obtained numerical solution of the this example for N=7,9 and N=11,13, respectively. From these figures, one can see the robustness of the proposed method for solving 1D diffusion equations.Figure 5
Error histories associated with the obtained numerical solution of the first example forN=7,9.Figure 6
Error histories associated with the obtained numerical solution of the first example forN=11,13.Example 2 (see [6]).
As the second example we consider (4)–(6) with the following assumptions: (34)ωη,τ=0;fη=sin⁡η;kη=0;g1τ=0;g2τ=0,with the exact solution ψη,τ=sin⁡ηcos⁡τ in the rectangle aη,bη×aτ,bτ=0,π×0,1. Again, we have implemented the BCM for solving this provided example by assuming several values of N 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 aη,bη=0,π, this assumption may decrease the accuracy of the TMM. Figures 7 and 8 illustrate eNη,τ associated with the obtained numerical solution of this example via BCM for N=6,8 and N=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 e12η,τ are around 10-11 approximately, while the associated error function values e12η,τ of the TMM [6] are around 10-7 (see the fourth figure of [6] for this aim).Table 1
Numerical results of the BCM and TMM associated with the erroreNη,τ in the second example.
( η , τ ) BCM TMM N = 8 N = 10 N = 12 N = 8 N = 10 N = 12 ( 0.6,0.6 ) 1.62 × 1 0 - 6 6.90 × 1 0 - 11 2.51 × 1 0 - 11 8.26 × 1 0 - 7 7.08 × 1 0 - 9 4.63 × 1 0 - 11 ( 0.7,0.7 ) 1.37 × 1 0 - 6 1.56 × 1 0 - 9 2.29 × 1 0 - 11 3.18 × 1 0 - 6 3.69 × 1 0 - 8 3.27 × 1 0 - 10 ( 0.8,0.8 ) 1.12 × 1 0 - 6 2.84 × 1 0 - 9 2.52 × 1 0 - 11 1.01 × 1 0 - 5 1.52 × 1 0 - 7 1.75 × 1 0 - 9 ( 0.9,0.9 ) 1.22 × 1 0 - 6 2.52 × 1 0 - 9 2.22 × 1 0 - 11 2.76 × 1 0 - 5 5.24 × 1 0 - 7 7.61 × 1 0 - 9 ( 1,1 ) 1.40 × 1 0 - 6 1.59 × 1 0 - 9 1.81 × 1 0 - 11 6.70 × 1 0 - 5 1.55 × 1 0 - 6 2.76 × 1 0 - 9Figure 7
Error histories associated with the obtained numerical solution of the second example forN=6,8.Figure 8
Error histories associated with the obtained numerical solution of the second example forN=10,12.Example 3 (see [5]).
As the third example we consider (22)–(24) with the following assumptions: (35)ωη,ξ,τ=0;fη,ξ=expη+ξ;g1ξ,τ=expξ+2t;g2ξ,τ=expξ+2t+1;h1η,τ=expη+2t;h2η,τ=expη+2t+1,with the exact solution ψη,ξ,τ=expη+ξ+τ in the unit cube aη,bη×aξ,bξ×aτ,bτ=0,1×0,1×0,1. The numerical results of this example are provided with the numerical results of the next example in Table 2. Moreover, numerical solution ψ71,ξ,τ, which is obtained via the proposed BCM, together with the exact solution ψ1,ξ,τ is depicted in Figure 9.Table 2
Comparisons of the BCM and BMM associated with the errorsEN in the third and fourth example.
N Third example Fourth example BCM BMM BCM BMM 4 3.2 × 1 0 - 2 3.5 × 1 0 - 2 7.3 × 1 0 - 3 2.8 × 1 0 - 2 6 1.6 × 1 0 - 3 2.2 × 1 0 - 3 2.9 × 1 0 - 3 6.1 × 1 0 - 3 8 8.2 × 1 0 - 5 1.7 × 1 0 - 4 2.7 × 1 0 - 4 1.4 × 1 0 - 3 10 3.4 × 1 0 - 7 1.4 × 1 0 - 5 4.2 × 1 0 - 5 3.7 × 1 0 - 4 12 1.8 × 1 0 - 9 1.1 × 1 0 - 6 7.9 × 1 0 - 6 9.1 × 1 0 - 5 14 5.4 × 1 0 - 11 9.2 × 1 0 - 8 1.2 × 1 0 - 6 2.3 × 1 0 - 5Figure 9
Numerical solutionψ71,ξ,τ and exact solution ψ1,ξ,τ of the third example.Example 4 (see [5]).
As the fourth example we consider (22)–(24) with the following assumptions: (36)ωη,ξ,τ=0;fη,ξ=cos⁡πηsin⁡πξ;g1ξ,τ=sin⁡πξexp-2π2τ;g2ξ,τ=-sin⁡πξexp-2π2τ;h1η,τ=0;h2η,τ=0,with the exact solution ψη,ξ,τ=cos⁡πηsin⁡πξexp-2π2τ in the unit cube aη,bη×aξ,bξ×aτ,bτ=0,1×0,1×0,1. For solving this example and the previous example, we applied several values of N such as 4, 6, 8, 10, 12, and 14 and achieve the numerical solution ψNη,ξ,τ by the suggested BCM. For these numerical solutions, we have defined the discrete error function (37)EN=max1≤i≤N+1,1≤j≤N+1ψηi,ξj,1-ψNηi,ξj,1for 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 ψ91,ξ,τ, which is obtained via the proposed BCM, together with the exact solution ψ1,ξ,τ is depicted in Figure 10.Figure 10
Numerical solutionψ9(1,ξ,τ) and exact solution ψ(1,ξ,τ) of the fourth example.


## 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, asN 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 x,y,z 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.



---

*Source: 5691452-2017-02-22.xml*