A boundary value approach for solving three-dimensional elliptic and hyperbolic partial differential equations.

In this article, the boundary value method is applied to solve three dimensional elliptic and hyperbolic partial differential equations. The partial derivatives with respect to two of the spatial variables (y, z) are discretized using finite difference approximations to obtain a large system of ordinary differential equations (ODEs) in the third spatial variable (x). Using interpolation and collocation techniques, a continuous scheme is developed and used to obtain discrete methods which are applied via the Block unification approach to obtain approximations to the resulting large system of ODEs. Several test problems are investigated to elucidate the solution process.

Background

This paper is devoted to the numerical computation of the three dimensional elliptic and hyperbolic PDEs. The general three dimensional problem (without cross derivatives) is given asdefined in the domian with boundary . Associated with (1) are the initial conditionsor the boundary conditionsand the Dirichlet boundary conditionsThe unknown function u, the variable coefficients a, b, c, d, e, f and the forcing term G are assumed to be continuously differentiable and have the required partial derivatives on .

The numerical approximation of (1) has received great attention from researchers in the past few decades. This is because most modeled physical and engineering processes results in PDEs of the form (1). The standard convection diffusion equation is obtained if and for all (see Ge and Zhang 2002). Many transport processes including fluid flows and heat transfer can be modeled by a convective diffusive equation which describes the convection and diffusion of various physical quantities such as heat and momentum (see Roache 1976). Several numerical schemes have been developed for (1) which include the compact difference scheme (Ge and Zhang 2002; Spotz 1995; Spotz and Carey 1996; Zhang 1998a, b, c, d; Zhang et al. 2000), compact alternating direct implicit scheme (Cui 2010; Liao and Sun 2010), finite difference methods (Jain et al. 1992; Mohanty and Singh 2006), method of lines (Brastos 2007; Dehghan and Shokri 2008) , Jacobi elliptic function method (Bhrawy and Abdelkawy 2013), the use of collocation and radial basis functions (Shakeri and Dehghan 2008), Riccati–Bernoulli sub ODE method (Yang et al. 2015), the use of the expansion methods (He 1998; Roshid and Rahman 2014; Alam et al. 2014; Roshid et al. 2013, 2014a, b, c) among others.

The method of lines approach is commonly used for solving PDEs whereby the PDEs is converted into a system of ODEs by replacing the appropriate derivatives by finite difference approximations. Our objective is to convert the 3D PDEs into a system of ODEs by replacing two of the spatial derivatives by central difference approximations. The resulting system is then solved using the BVM. Specifically, we discretize the y and z variable with mesh spacingsWe then define the vectorsandwhere and .

Also, we replace the partial derivatives , , and occuring in (1) by central difference approximationsProblem (1) leads to the semidiscretized systemwhich can be written in the formsubject to the initial conditionsor the boundary conditionswhere and is a matrix arising from the semidiscretized system (2) which is expressed in the form (3) and whose solution is sought using the BVMs. We note that is a vector of constants. The BVMs are a class of linear multistep methods (LMM) with step number k and whose k additional conditions are not only imposed at the beginning of the integration process but also at the end so that they form a discrete analog of the continuous boundary value problems. Thus, they are used for the numerical approximation of both initial and boundary value problems. They have been used for the solution of first order 1D initial and boundary value problems and their convergence and linear stability properties have been fully discussed in Brugnano and Trigiante (1998). Recently, Biala and Jator (2015) developed BVMs for the direct solution of systems of the general second order ODEs. One main feature of the BVMs is that they can be used in the same way for solving both initial and boundary value problems. Therefore, such methods are the best candidate for solving the semididscretized PDEs in (3). The BVMs simultaneously generates approximate solutions to the exact solution on the entire interval of integration. This approach has the advantage of producing smaller global errors than those produced by the step-by-step methods due to the presence of accumulated errors at each step in the step-by-step method.

The paper is organized as follows: in "Description of the BVM", we derive a continuous scheme which is used to formulate the BVM as well as investigate the properties of the BVM. The computational complexities associated with the method is addressed in "Computational procedures". Several numerical test examples are given in "Numerical examples" to show the accuracy of the method. We give some concluding remarks in "Conclusion".

Description of the BVM

We propose a BVM for (1) in which on the partition the four step is given by the equationswhere and , , and are coefficients to be uniquely determined. We note that denote the numerical approximation to the analytical solution and .

Development of the continuous BVM

In this section, we discuss the construction of the continuous scheme, via the interpolation and collocation approach (Jator and Li 2012), from which (6) is derived.

Theorem 1

Let the continuous representationsatisfy the following conditionsthen the continuous representation (7) is equivalent towhere we define the matrix V asis obtained by replacing thejth column ofV byW whereTdenotes the transpose,are basis functions and W is a vector given by

Proof

We define the polynomial basis functionswhere and are coefficients to be uniquely determined.

Substituting (10) into (7), we havewhich may be written asand expressed in the formwhereImposing conditions (8) on (11), we obtain a system of seven equations, which can be expressed as where is a vector of seven undetermined coefficients. Using the Crammer’s rule, the elements of L can be obtained and are given aswhere is obtained by replacing the jth column of V by W. We rewrite (11) as (9) using the newly found elements of L.

BVM and its block extension

The coefficients given in (6) are specified by evaluating (9) at and evaluating at to obtain

Remark 1

We note that the method (12) is locally obtained on and is applied to simultaneously obtain approximations to the semidiscretized system (3) over the whole 3D space in which case , and . Also, we note that the first three formulas in (12) are of while the derivative formulas are of .

Convergence analysis

We discuss the convergence of the BVMs in the following theorem

Theorem 2

Let be an approximation of the solution vector for the system obtained on a partition from the method (12). If, , where the exact solution is several times differentiable on and if, then, the BVM is convergent of order 6, which implies that .

We begin the proof by compactly writing (12) in matrix form with the introduction of the following matrix notations. Let A be a matrix defined by and are null and identity matrices respectively.

Similarly, let B be a matrix defined bywhere are given as and are null matrices.

We also define the vectorswhere L(h) is the local truncation error vector of the formulas in (12).

The exact form of the system formed by (12) is given byand the approximate form of the system is given bywhere is the approximate solution of . Subtracting (13) from (14) and letting and using the Mean value theorem, we have the error systemwhere J is the Jacobian matrix and its entries , , and are defined asLet be a matrix of dimension 2L so that (15) becomesand for sufficiently small h, is a monotone matrix and thus nonsingular. Thereforeandwhich shows that the methods are convergent and the global error is of order .

Computational procedures

The method (12) can also be expressed in block form aswhere the positive integer is the number of blocks,and , , , and are matrices each of dimension 8 whose entries are given by the coefficients of (12).

Equation (1) is converted to (3) by discretizing the partition and , given bywhereare constant stepsizes of the partition and respectively, , M and N are positive integers and m, n are the grid index in the y and z direction respectively.

The resulting system of ODEs (3) is then solved on the partition . The block unification of (17) lead to a large system of finite difference equations which is then solved to provide all solutions of (3) on .

The following algorithm summarizes the numerical integration of (3) for some set of points x on .

Step 1 Use the block unification of (17) for to obtain on the domian and for , is obtained on the domain , and on the domains , for , , we obtain .

Step 2 The unified block given by the system obtained in Step 1 is a large system with dimensions with , , ,

Step 3 The system obtained in Step 2 is solved using the feature NSolve in Mathematica 10.0 for linear problems and FindRoot (which incorporates the Newton’s method) for nonlinear problems.

Step 4 The solution of (1) is approximated by the solution in Step 3 as ,

Numerical examples

In this section, some examples are investigated to show the reliability and efficiency of the proposed scheme in this paper.

Test 1

We consider the Laplace equation with non zero forcing term G(x, y, z) and with zero boundary values on the entire given in Zhang (1998c) and whose solution is given aswhere . Table 1 shows the errors in the norm with different mesh sizes. Figure 1 shows the plot of the exact, approximate and error function when

Table 1

Errors in the norm for test problem 1

L_4	Δy=Δz=0.25	Δy=Δz=0.125	Δy=Δz=0.0625	Δy=Δz=0.03125
0.1	1.903e−04	4.905e−05	1.242e−05	3.109e−06
0.2	1.413e−03	3.641e−04	9.208e−05	2.305e−05
0.3	4.238e−03	1.091e−03	2.754e−04	6.896e−05
0.4	8.610e−03	2.213e−03	5.576e−04	1.396e−04
0.5	1.401e−02	3.586e−03	9.028e−04	2.260e−04
0.6	1.972e−02	5.010e−03	1.263e−03	3.158e−04
0.7	2.506e−02	6.291e−03	1.592e−03	3.980e−04
0.8	2.953e−02	7.340e−03	1.855e−03	4.636e−04
0.9	3.280e−02	8.142e−03	2.029e−03	5.092e−04
1.0	3.472e−02	8.616e−03	2.144e−03	5.356e−04

Fig. 1

Graphical evidence when and for test problem 1

Test 2

Next, we consider the following Sine-Gordon equation given in Cui (2010)where and whose theoretical solution is . In our computations, we have chosen and the norms are given in Table 2 with different meshsizes. Figure 2 also shows the plot of the exact, approximate and error function when

Table 2

Errors in the norm for test problem 2

L_4	Δy=Δz=0.5	Δy=Δz=0.25	Δy=Δz=0.125	Δy=Δz=0.0625
0.2	1.266e−01	3.343e−02	8.471e−03	2.125e−03
0.4	4.016e−01	1.035e−02	2.609e−02	26.535e−03
0.6	6.692e−01	1.637e−02	4.067e−02	1.015e−02
0.8	8.183e−01	1.811e−01	4.406e−02	1.112e−02
1.0	8.115e−01	1.815e−01	4.406e−02	1.013e−02
1.2	8.344e−01	1.814e−01	4.406e−02	1.015e−02
1.4	7.892e−01	1.788e−01	4.380e−02	1.187e−02
1.6	8.224e−01	1.810e−01	4.380e−02	1.055e−02
1.8	8.429e−01	1.764e−01	4.371e−02	1.034e−02
2.0	8.308e−01	1.815e−01	4.408e−02	1.154e−02

Fig. 2

Graphical evidence when and for test problem 2

Test 3

We also consider the following equation given in Cui (2010)where and whose theoretical solution is with corresponding Dirichlet boundary conditions. The norms are given in Table 3 with different meshsizes. Figure 3 also shows the plot of the exact, approximate and error function when

Table 3

Errors in the norm for test problem 3

L_4	Δy=Δz=0.5	Δy=Δz=0.25	Δy=Δz=0.125	Δy=Δz=0.0625
-0.8	5.340e−04	1.493e−04	3.896e−05	9.775e−06
-0.6	2.042e−03	5.911e−04	1.493e−03	3.763e−05
-0.4	4.878e−03	1.257e−03	3.167e−04	7.991e−05
-0.2	7.896e−03	2.007e−03	5.273e−04	1.320e−05
0.0	1.017e−02	3.005e−03	7.539e−04	1.895e−04
0.2	1.502e−02	3.966e−03	9.964e−04	2.494e−04
0.4	1.739e−02	4.742e−03	1.231e−03	3.097e−04
0.6	2.005e−02	5.389e−03	1.436e−03	3.648e−04
0.8	2.324e−02	5.975e−03	1.514e−03	3.793e−04
1.0	1.995e−02	5.048e−03	1.349e−03	3.449e−04

Fig. 3

Graphical evidence when and for test problem 3

Test 4

We also consider the singularly perturbed convection diffusion equation given in Mohanty and Singh (2006)defined in the domain with boundary and subject to the Dirichlet boundary conditions and whose theoretical solution iswhere . We have solved the problem using to 1.0. The norms are given in Table 4 with different meshsizes. Figure 4 also shows the plot of the exact, approximate and error function when

Table 4

Errors in the norm for test problem 4

ϵ	Δy=Δz=0.25	Δy=Δz=0.125	Δy=Δz=0.0625	Δy=Δz=0.03125
0.1	1.492e−01	2.086e−02	1.801e−03	2.767e−04
0.2	2.830e−02	5.918e−03	1.108e−03	2.626e−04
0.34	1.890e−2	4.807e−03	1.113e−03	2.742e−04
0.4	1.655e−02	4.744e−03	1.146e−03	2.843e−04
0.5	1.680e−02	4.805e−03	1.174e−03	2.916e−04
0.6	1.737e−02	4.867e−03	1.195e−03	2.916e−04
0.7	1.780e−02	4.919e−03	1.211e−03	3.970e−04
0.8	1.813e−02	4.962e−03	1.224e−03	3.048e−04
0.9	1.839e−02	4.998e−03	1.235e−03	3.075e−04
1.0	1.860e−02	5.028e−03	1.243e−03	3.097e−04

Fig. 4

Graphical evidence when for test problem 4

Test 5

The singularly perturbed elliptic boundary value problem given in Mohanty and Singh (2006)defined in the domain with boundary and subject to the Dirichlet boundary conditions and where the forcing term G is set to satisfy the exact solution given as . We have solved the problem using . The norms are given in Table 5 with different meshsizes and for different values of . Figure 5 also shows the plot of the exact, approximate and error function when

Table 5

Errors in the norm for different values of for test problem 5

α	Δy=Δz=0.25	Δy=Δz=0.125	Δy=Δz=0.0625	Δy=Δz=0.03125
0	1.359e−04	3.664e−05	9.521e−06	2.395e−06
1	1.383e−04	3.725e−05	9.768e−05	3.226e−06
2	1.380e−04	3.951e−05	9.858e−05	1.892e−06
3	1.370e−04	4.582e−05	1.649e−05	8.940e−05
4	1.361e−04	5.382e−05	2.623e−05	1.892e−05
5	1.359e−04	6.466e−05	3.936e−05	3.226e−05

Fig. 5

Graphical evidence when for for test problem 5

Conclusion

In this paper, we have developed a highly accurate 3D problem solver. This has been achieved by the discretization of two of the spatial variables and the construction of a continuous BVM via the interpolation and collocation approach for solving the resulting semidiscretized system. The results given in "Numerical examples" show that the approach is highly efficient and highly accurate.