On a bivariate spectral relaxation method for unsteady magneto-hydrodynamic flow in porous media

The paper presents a significant improvement to the implementation of the spectral relaxation method (SRM) for solving nonlinear partial differential equations that arise in the modelling of fluid flow problems. Previously the SRM utilized the spectral method to discretize derivatives in space and finite differences to discretize in time. In this work we seek to improve the performance of the SRM by applying the spectral method to discretize derivatives in both space and time variables. The new approach combines the relaxation scheme of the SRM, bivariate Lagrange interpolation as well as the Chebyshev spectral collocation method. The technique is tested on a system of four nonlinear partial differential equations that model unsteady three-dimensional magneto-hydrodynamic flow and mass transfer in a porous medium. Computed solutions are compared with previously published results obtained using the SRM, the spectral quasilinearization method and the Keller-box method. There is clear evidence that the new approach produces results that as good as, if not better than published results determined using the other methods. The main advantage of the new approach is that it offers better accuracy on coarser grids which significantly improves the computational speed of the method. The technique also leads to faster convergence to the required solution.

(SQLM). The HAM has been used extensively by researchers working on such problems Abbas et al. (2008), Ahmad et al. (2008), Ali and Mehmood (2008), Mehmood et al. (2008), Alizadeh-Pahlavan and Sadeghy (2009), Fan et al. (2010), Xu et al. (2007), You et al. (2010). It is an analytic method for approximating solutions of differential equations developed by Liao (2012). The homotopy analysis method is an analytic method where accuracy and convergence are achieved by increasing the number of terms of the solution series. In some cases, such as when a large embedded physical parameter multiplies the nonlinear terms, far too many terms may be required to give accurate results. Retaining too many terms in the solution series is cumbersome, even with the use of symbolic computing software. The use of the HAM further depends on other arbitrarily introduced parameters such as the convergence controlling parameter which must be carefully selected through a separate procedure.
A popular numerical method used by many researchers to solve unsteady boundary layer flow problems is the Keller-box method Ali et al. (2010a, b), Lok et al. (2010), Nazar et al. (2004a, b). The Keller-box method is a finite difference based implicit numerical scheme which was developed by Cebeci and Bradshaw (1984). Recently, Motsa et al. (2014a, b) used spectral based relaxation and quasilinearization schemes to solve unsteady boundary layer problems. These schemes are accurate, easy to implement and are computationally efficient. As observed in Motsa et al. (2014a), the limitation of the spectral quasilinearization method is that the coupled high-order system of differential equations may often lead to very large systems of algebraic equations that may require significant computing resources. In addition, the actual process of developing the solution algorithm is time-consuming in comparison to SRM. This is because with SQLM, the process begins with the quasi-linearization step whereas with SRM the iteration scheme is obtained directly by requiring some terms to be evaluated at the current iteration and others at the previous iteration. The SRM works much like the familiar Gauss-Seidel iteration by decoupling a system of non-linear PDEs into a system of linear PDEs which are then solved in succession. Consequently, the SRM is easy to implement and computationally efficient.
Both the original SRM and SQLM used in Motsa et al. (2014a) use finite differences to discretize derivatives in time. This is a disadvantage because finite difference schemes are known to converge slower than spectral methods. The use of finite differences effectively nullifies the benefits of fast convergence when spectral collocation is used to discretize in space. Furthermore, finite differences require fine grids with very small step sizes to guarantee accuracy, hence there is a huge computation time overhead each time the grid is refined. This paper provides a different approach to the implementation of the spectral relaxation method introduced in Motsa et al. (2014a). The innovation is that the spectral collocation method is used to discretize derivatives in both space and time. As a result, there are uniform convergence benefits in both directions. The scheme uses fewer grid points in space and time and thus it converges very fast. We refer to the improved SRM as the bivariate interpolation spectral relaxation method (BI-SRM). To test the viability of this innovation as a solution method, we have solved the coupled system of third and second order partial differential equations that describe a boundary-layer system. A careful comparison of the new results is made with the earlier SRM, SQLM and the Keller-box results reported in Motsa et al. (2014a). We particularly compare the computational times for the different methods to reach the same level of accuracy.

Model equations
We consider the unsteady and three-dimensional flow of a viscous fluid over a stretching surface investigated by Hayat et al. (2010). The fluid is electrically conducting in the presence of a constant applied magnetic field B 0 . The induced magnetic field is neglected under the assumption of a small magnetic Reynolds number. The flow is governed by the following four dimensionless partial differential equations with the following boundary conditions In the above equations prime denotes the derivative with respect to η, c the stretching parameter is a positive constant. M is the local Hartman number, the local porosity parameter, Sc the Schmidt number, Pr the Prandtl number and γ the chemical reaction parameter. The initial unsteady solution can be found exactly by setting ξ = 0 in the above equations and solving the resulting equations. The closed form analytical solutions are given by . Magagula et al. SpringerPlus (2016) 5:455 Bivariate interpolated spectral relaxation method (BI-SRM) In this section we introduce the Bivariate Interpolated Spectral Relaxation Method (BI-SRM) for solving the system of nonlinear partial differential equations (1)-(3). Applying the relaxation scheme Motsa et al. (2014a) to the system of nonlinear partial differential equations gives the following linear partial differential equations; subject to where the variable coefficients are given by and r and r + 1 denote previous and current iterations respectively. The system of linear partial differential equations (11)-(14) is discretised using the Chebyshev spectral collocation both in space (η) and time (ξ) directions. The Chebyshev collocation method is valid in the domain [−1, 1] in space and time. Therefore, the physical region, ξ ∈ [0, 1] is converted to the region t ∈ [−1, 1] using a linear transformation and similarly, η ∈ [0, L ∞ ] is converted to the region x ∈ [−1, 1]. The system of linear partial differential equations (11)-(14) is decoupled. Therefore, each equation can be solved independently of the other equations in the system. We assume that the solution to Eq. (11) can be approximated by a bivariate Lagrange interpolation polynomial of the form which interpolates f (η, ξ) at selected points in both the η and ξ directions defined by The Chebyshev-Gauss-Lobatto grid points (17) ensures that there is a simple conversion of the continuous derivatives, in both space and time, to discrete derivatives at the grid points. The characteristic Lagrange cardinal polynomial L p (x) is defined as where Similarly, we define the function L j (t). Equation (16) is then substituted into Eq. (11). An important step in the implementation of the solution procedure is the evaluation of the derivatives of L p (x) and L j (t) with respect to x and t respectively. The derivative of f (η, ξ) with respect to ξ at the Chebyshev-Gauss-Lobatto points (x s , t i ), is computed as are the entries of the standard first derivative Chebyshev differen- Trefethen (2000) for i, j = 0, 1, . . . , N t . Similarly, we compute the derivative of f (η, ξ) with respect to η at the Chebyshev-Gauss-Lobatto points (x s , t i ), as follows dx , are the entries of the standard first derivative Chebyshev differentiation matrix of size (N x + 1) × (N x + 1). Therefore, an nth order derivative with respect to η is given by The vector F i is defined as where the superscript T denotes matrix transpose. Collocating using Eqs. (24) and (21) on (11), we get The boundary equations are given by The initial unsteady solution given by equation (7) corresponds to t = t N t = −1. Therefore, we evaluate Eq. (26) for i = 0, 1, . . . , N t − 1. Equation (26) can be expressed as where and F N t is the known initial unsteady solution given by equation (7). Imposing boundary conditions for i = 0, 1, . . . , N t − 1, Eq. (28) can be expressed as the following where Imposing initial boundary conditions and applying the Chebyshev bivariate collocation as described above on Eqs. (12), (13) and (14) we get where and the vectors G N t , N t and N t are the known initial unsteady solutions given by Eqs. (8), (9) and (10) respectively. Imposing boundary conditions for i = 0, 1, . . . , N t − 1, equations (32), (33) and (34) can be expressed as the following N t (N x + 1) × N t (N x + 1) matrix system where and I is the standard (N x + 1) × (N x + 1) identity matrix. We obtain the numerical solutions for g(η, ξ), θ(η, ξ) and φ(η, ξ) by solving matrix equations (38), (39) and (40) iteratively for r = 1, 2, . . . M, where M is the number of iterations to be used. Equations (8), (9) and (10) are used as initial guesses.

Results and discussion
In this section we present the numerical solutions of the three dimensional unsteady three dimensional magneto-hydrodynamic flow and mass transfer in a porous media obtained using the BI-SQLM algorithm. In our computations the η domain was truncated to η ∞ = 20. This value gave accurate results for all the quantities of physical interest. To get accurate solutions, N x = 60 collocation points were used to discretize the space variable η and only N t = 10 collocation points were enough for the time variable ξ.
As earlier mentioned, this problem has been solved before by Motsa et al. (2014a) using the spectral relaxation method (SRM), spectral quasilinearization method (SQLM) and the Keller-box method. The results from their paper combined with the present results of the BI-SRM are shown in Table 1. It can be observed from the table that the Keller-box method takes a significant amount of computational time than the SRM and SQLM. This is because the Keller-box is entirely based on finite difference schemes while the SRM and SQLM only uses finite differences in the time variable. In the space variable both the SRM and SQLM use spectral method. It is well documented from literature   that spectral methods converge very fast when the solution is smooth. This brought about the idea of using spectral methods in both space and time to increase efficiency. The BI-SRM discretizes both the space and time domains using spectral methods. From the results shown in the table it is evident that the BI-SRM is by far superior than the other methods in terms of computational time taken to reach the same level of accuracy.
In Table 1, we also show the number of grid points required by each of the methods to discretize in time. All the finite difference based discretizations required 2000 grid points compared to the spectral discretization of the BI-SRM which required only 10 grid points to reach the same level of accuracy. The grid independence test for the algorithm is shown in Table 2. The skin friction values, Nusselt number and Sherwood numbers are the variables used in carrying out the grid independence test in Table 2.
The residual error graphs of Eqs.
(1)-(3) are presented in Figs. 1, 2, 3 and 4 respectively. In Figs. 1 and 2, we observe that the residual error is reduced with an increase in the iterations of the scheme. The rate of reduction of the residual error appears to be linear. The residual error is minimum at ξ = 0 and is increased sharply near 0 until a certain level is reached after which it is almost constant. The residual error appears to be nearly uniform in 0 < ξ ≤ 1 or increases only slightly. It is also observed that the Table 2 Values of f ′′ (0, ξ), g ′′ (0, ξ), θ ′ (0, ξ) and φ ′ (0, ξ) when = 0.5, M = 2, c = 0.5, Sc = γ = 1, Pr = 1.5 order of magnitude of the residual error can be seen to be small in the 0 ≤ ξ ≤ 1 interval. Lastly, after only two iterations the residual error appears to be less than 0.01 in the entire range of ξ . The small residual error using only a few iteration points to the accuracy of the method. This error can be decreased at a linear rate with an increase in the number of iterations. The decrease in the error with additional iterations suggests that the iteration scheme converges. It should be noted that when ξ = 0, governing equations reduce to a linear system that can be solved directly using the spectral collocation method with discretization only in η without the use of relaxation and iterations. This explains why the best accuracy is observed at ξ = 0. The near uniformity of the residual error in 0 < ξ ≤ 1 can be attributed to the use of Lagrange polynomial basis functions whose error is known to be uniformly distributed in the interpolating region. We can therefore conclude that the method gives accurate results, the rate of convergence of the method is linear and that the method requires only a few iterations to give very accurate results.
We observe that the residual error for the first iteration appears to be very small in Figs. 3 and 4. The residual error is the same for all iterations greater than one. The residual error increases with an increase in ξ. We also observe that the residual error is smaller than the one for f and g even when fewer iterations are used. The observation that the residual error for the momentum and energy is small even for the 1st iteration is perhaps the most interesting finding of the study. This means that when using the proposed approach, the best possible results that can be achieved by the method can be obtained after using just one iteration. Further increase in the number of iterations doesn't improve the accuracy of the solution. After one iteration of the momentum equations for f and g the energy and mass transfer equations reduce to linear homogeneous equations whose solution appears to be marginally influenced by variations in f r and g r for r > 1. Since with just one iteration we obtain extremely accurate results for θ and φ, the implication is that in solving for energy and momentum equations for such a problem, it is not necessary to iterate. It is enough to just use the initial approximation. Is is worth noting that the energy and mass transfer equations are homogeneous equations in θ and φ respectively. It is possible that the findings obtained in this study are only applicable in such equations. This has to be investigated further.
The convergence graphs of Eqs.
(1)-(3) are presented in Figs. 5, 6 respectively. In Figs. 5 and 6, the residual error decreases linearly with an increase in the number of iterations. The residual error is smallest when ξ is near zero and largest when ξ is large. This is seen from the convergence level which is near 10 −20 for ξ = 0.25 and about 10 −15 for ξ = 0.75. The slope of the residual error graphs is the same for all values of ξ. Full convergence is achieved after 13 iterations for both ξ = 0.75 and ξ = 1. For ξ = 0.25 full convergence is achieved after 16 iterations but at a much smaller magnitude of residual error. The decrease in the residual error with increase in iterations suggests that the iteration scheme converges. Small residual error near zero suggests that best accuracy (after full convergence) is observed near zero. The method converges (fewer iterations needed to attain full convergence) at or near ξ = 1. However, the convergence efficiency doesn't translate to better accuracy because, as can be seen for the case of ξ values near zero, the convergence level is 10 −16 . The same slope for all the graphs means that the convergence rates of the method is the same for all values of ξ. The method is convergent and very accurate in whole time domain ξ ∈ [0, 1] which translates to τ ∈ [0, ∞). The method converges with nearly the same convergence rate for all values of ξ. The method gives the best accuracy near ξ = 0 and less accurate, comparatively, at or near ξ = 1. We note that even at ξ = 1, the method gives very accurate results with a residual error norm of about 10 −15 . This is one of the highlights of this investigation.
The accuracy of the solutions for energy and mass transfer equations are not dependent on successive relaxation or iterations of the momentum equations since the convergence of the solutions doesn't improve at all with an increase in the number of iterations. Hence, the results for θ(η) are not dependent on successive approximations for f (η) and g(η). It was observed that it is enough to use one iteration of f (η) to give accurate result for θ(η). In Table 1, we present the results we obtained using the algorithm. The skin friction values, Nusselt number and Sherwood number for various values of ξ are displayed in Table 1. The results obtained using the method match those obtained using other methods. Computing time for the BI-SRM is much smaller than the other methods it is compared with. The method gives valid results when used with few collocation points. In particular, the BI-SRM requires only ten grid points to achieve the same valid results. The table validates the results obtained in this study. BI-SRM is computationally fast in generating valid results when compared with the SRM, SQLM and Keller-Box. We can infer that the BI-SRM is better than finite differences coupled SRM in terms of computational speed and accuracy and the better accuracy could be the result of applying spectral collocation with uniform accuracy level in both η and ξ directions.
In Tables 3 and 4, the residual errors and convergence rates of g and g when = 0.5, M = 2, c = 0.5, Sc = γ = 1, Pr = 1 are displayed. We observe that the convergence rate is linear and the residual error is smaller near ξ = 0. Table 3 The residual errors and convergence rates of f when = 0.5, M = 2, c = 0.5, Sc = γ = 1, Pr = 1.5 Iter.