International Journal of Mathematics and Mathematical Sciences
Volume 2019 (2019), Article ID 6423294, 18 pages
https://doi.org/10.1155/2019/6423294

Research Article

On the Multidomain Bivariate Spectral Local Linearisation Method for Solving Systems of Nonsimilar Boundary Layer Partial Differential Equations

Vusi Mpendulo Magagula1,2

¹Faculty of Science and Engineering, University of Eswatini, Private Bag 4, Kwaluseni M201, Matsapha, Eswatini
²DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa

Correspondence should be addressed to Vusi Mpendulo Magagula; gutjwa@gmail.com

Received 10 January 2019; Accepted 15 May 2019; Published 10 June 2019

Academic Editor: Irena Lasiecka

Copyright © 2019 Vusi Mpendulo Magagula. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In this work, a novel approach for solving systems of nonsimilar boundary layer equations over a large time domain is presented. The method is a multidomain bivariate spectral local linearisation method (MD-BSLLM), Legendre-Gauss-Lobatto grid points, a local linearisation technique, and the spectral collocation method to approximate functions defined as bivariate Lagrange interpolation. The method is developed for a general system of nonlinear partial differential equations. The use of the MD-BSLLM is demonstrated by solving a system of nonlinear partial differential equations that describe a class of nonsimilar boundary layer equations. Numerical experiments are conducted to show applicability and accuracy of the method. Grid independence tests establish the accuracy, convergence, and validity of the method. The solution for the limiting case is used to validate the accuracy of the MD-BSLLM. The proposed numerical method performs better than some existing numerical methods for solving a class of nonsimilar boundary layer equations over large time domains since it converges faster and uses few grid points to achieve accurate results. The proposed method uses minimal computation time and its accuracy does not deteriorate with an increase in time.

1. Introduction

Solutions of nonlinear nonsimilar boundary layer partial differential equations have received considerable attention of research scientist in the past few years. Series solution methods have been used by several researchers. These include asymptotic methods [1, 2], extended series solutions [2], and homotopy analysis method (HAM) [3–8]. To achieve accurate results when using series based methods, a large number of terms are required. A complex nonsimilar boundary layer differential equation requires many terms and thus solving the problem becomes tedious. Simplification of a complex model may introduce errors and may lead to inaccurate solutions.

Numerical methods have been used to solve nonsimilar boundary layer equations. Numerical methods like implicit finite difference methods have been used to solve nonsimilar boundary layer equations by Hossain [1, 9–11]. Hussain [2] and Yih [12] used the Keller-Box method. Finite difference methods require many grid points to achieve accurate results. The accuracy of the results deteriorates fast as the time domain increases for parabolic nonlinear partial differential equations. Recently, Magagula et al. [13] showed that pseudospectral methods can achieve accurate results with few grid points compared to Keller-Box method. Results were obtained in a fraction of a second.

Pseudospectral methods have been used to solve boundary layer ordinary differential equations [14–19], and boundary layer partial differential equations [20–22] have been solved using spectral local linearisation method (SLLM) and the bivariate spectral local linearisation method (BSLLM), respectively. Other pseudospectral methods have been used to solve differential equations [23–25]. The bivariate spectral local linearisation method was developed and generalized by Motsa [16] for solving systems of boundary layer ordinary differential equations. This method uses quasilinearisation to linearise nonlinear coupled system of differential equations. The result is a system of linear decoupled differential equations which can be solved independently. The SLLM was extended by Motsa [21] to a system of nonlinear coupled boundary layer partial differential equations, where spectral collocation was used in both space and time. The BSLLM method is suitable for nonlinear partial differential equations with small time domain. Increasing the time domain reduces the accuracy of the method. In some cases, the computed numerical solutions fail to converge even when the time step is decreased. Thus, accurate solutions of systems of differential equations are usually determined over short time intervals. The realization that accurate solutions can be obtained in small time intervals have led to a number of numerical methods that uses the domain-decomposition technique. The domain-decomposition approach together with the spectral collocation technique has been used to solve equations arising from various scientific fields. Funaro [26] and Harald [27] solved elliptic differential equations, Kopriva [28] solved hyperbolic systems, Peyret [29] solved stiff problems in fluid mechanics, and recently, it has been used to solve chaotic system of ordinary differential equations by [30–33] and Lane-Emden type of equations [34]. These methods split large intervals into smaller subintervals so that the equations are solved over a sequence of nonoverlapping subintervals of the domain. Continuity conditions are used to advance the solution across the nonoverlapping subintervals.

In this paper, we present a multidomain numerical approach for finding solutions of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval. The multidomain bivariate spectral local linearisation method (MD-BSLLM) is an extension of the BSLLM. The proposed approach is developed for solving system of coupled nonsimilar boundary layer partial differential equations. The system of equations is solved iteratively using a Jacobi like relaxation approach. The domain is divided into smaller nonoverlapping subintervals on which the Chebyshev spectral collocation method is used to solve the equations. A continuity condition is used to advance the solution across the subintervals. The algorithm is easy to develop and it yields very accurate results using only a few discretization nodes. The accuracy of the method is validated against the series solution for limiting cases. The main aim of the study is to explore the applicability of the multidomain bivariate spectral local linearisation collocation method to systems of nonlinear coupled nonsimilar boundary layer partial differential equations over large time intervals. The results confirm that the method is suitable for solving all types of systems of nonlinear coupled nonsimilar boundary layer partial differential equations over a large time interval.

The article is organized as follows: in Section 2, the multidomain bivariate spectral local linearisation method is described for a general system of partial differential equations. In Section 3, numerical experiments and the series solutions are presented in Section 4. In Section 5, the results are discussed and we conclude in Section 6.

2. Multidomain Bivariate Spectral Local Linearisation Method

In this section, the multidomain bivariate spectral local linearisation method (MD-BSLLM) for approximating solutions of system of nonlinear partial differential equations that model nonsimilar boundary layer equations is introduced and explained in detail. The multidomain bivariate spectral local linearisation method is developed for a general system of nonlinear partial differential equations. Consider a system of nonlinear partial differential equations of the form,

where

The order of differentiation is denoted by ; the required solutions by for and for are the nonlinear operators containing all the spatial derivatives and time derivatives of . The Legendre-Gauss-Lobatto grid points and the corresponding differentiation matrices are defined in the interval . Therefore, the time interval and space region are transformed to using linear transformations and , respectively. Thus . In order to apply the multidomain Legendre-Gauss-Lobatto bivariate spectral local linearisation method on the system of equations under consideration, the interval is first divided into nonoverlapping subintervals for . Each interval is further divided into not necessarily equal divisions. The system of nonlinear partial differential equations under consideration is solved in each subinterval . The solutions in each subinterval are denoted by , for . In the first interval , the solutions , for are obtained subject to the “initial" conditions , for , respectively. For each , at each subinterval , the continuity conditions

are used to implement the MD-BSQLM over the interval . This procedure is repeated to generate a sequence of solutions , for and . It is assumed that, at each subinterval , the solution can be approximated by a bivariate Lagrange interpolation polynomial of the form

for and . The bivariate Lagrange interpolation polynomial interpolates at selected points in both the and directions, for and . The grid points in time are given by the zeros of the polynomial

where is the first derivative of the Legendre polynomial given by

The zeros of can be computed using any root finding method. In this work, we use Newton’s method. Hence given and for , , we have

We note that

and thus

Substituting (15) into (13), we get

Similarly, the grid points in space are given by the zeros of the polynomial

and hence the zeros of can be obtained by using the modified Newton method

The zeros of (11) and (17) are called Legendre-Gauss-Lobatto points [35, 36]. The function is the characteristic Lagrange cardinal polynomial with the Legendre-Gauss-Lobatto points [35, 36]

where

An alternative convenient expression for (19) is given by

The nonlinear operators , for , are first linearised using the quasilinearisation technique as defined by Bellman and Kalaba [37]. The quasilinearisation method is a technique based on the Taylor series expansion of about some previous iteration. It is assumed that the differences between its previous and current solution and all its derivatives are small. Applying the quasilinearisation method yields the following:

where and denote previous and current iterations, respectively, and is a vector of the partial derivatives which is defined as

We define

where the prime denotes differentiation with respect to . The linearised equation (22) can be expressed in a compact form as

Equations (28)-(31) form a system of decoupled linear partial differential equations. They are solved iteratively for Equations (28)-(31) can further be expressed as follows:

where , , and are variable coefficients of , , and , respectively, in the subinterval of the multidomain approach, for and . These variable coefficients correspond to the equation, for . The constant denotes the order of differentiation. Thus, we have

In general, the th right hand side is given by

Equations (32), (33), and (35) are evaluated at the Legendre-Gauss-Lobatto grid points and . The values of the time derivatives are computed at the Legendre-Gauss-Lobatto points , as (for )

where is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size , given by [36]

The values of the space derivatives at the Legendre-Gauss-Lobatto points (for ) are similarly computed as

where is the standard first derivative Legendre-Gauss-Lobatto based differentiation matrix of size as defined in [36]. The entries of the space first derivative matrix with respect to the Legendre-Gauss-Lobatto points are exactly the same as those given by (42). Higher th order derivatives are defined as

where the vector is defined as

and the superscript denotes matrix transpose. Substituting (44) and (45) into (32), (33), and (35) yields

where

The diagonal matrices of the corresponding variable coefficients are given by

Imposing boundary conditions for , (47), (48), and (50) can be expressed as the following matrix system

where

The vector is defined as

for and . The vector corresponds to the initial boundary condition which is always prescribed.

3. Numerical Experiments

In this section, we present systems of equation for various fluid flow models considered in this study. The MD-BSLLM algorithm is applied to each system of equations as a numerical experiment.

3.1. Steady Free Convection Flow Past a Nonisothermal Vertical Porous Cone

In this subsection, the problem of steady two-dimensional laminar free convection flow past a nonisothermal vertical porous cone with variable surface temperature is also considered for numerical experimentation. The governing nonsimilarity system of partial differential equations is expressed in dimensionless form as (see Hossain et. al. [1] for derivation)

where is the Prandtl number, is the dimensionless suction parameter, and and are the dimensionless temperature and stream functions, respectively. The appropriate corresponding boundary conditions are

The skin friction coefficient and the Nusselt number describe the shear-stress and heat flux rate at the surface, respectively, and are defined by [1] as

In this case, the order of differentiation and the number of system of equations . In each subinterval , we have the equation

where

We note that the equations are in terms of the functions and . In order to apply our method, we let and . Therefore, we have

The variable coefficients in (64) are given by

and the variable coefficients for (65) are given by

Equations (62) and (63) can be expressed respectively as the following matrix system:

where

The vectors and are defined as

for .

3.2. Steady Free Convective Flow of a Viscous Incompressible Fluid along a Permeable Vertical Flat Plate

A two-dimensional steady free convective flow of a viscous incompressible fluid along a permeable vertical flat plate in the presence of soluble species is being considered in this subsection. The problem considered here is that of a natural convection boundary layer flow, influenced by the combined buoyancy forces from mass and thermal diffusion from a permeable vertical flat surface with nonuniform surface temperature and nonuniform surface species concentration, but with a uniform rate of suction of fluid through the permeable surface. The transformed boundary layer equations are solved numerically near to and far from the leading edge, using the multidomain bivariate spectral collocation method. The equations were solved by Hussain [2] and are

The appropriate boundary conditions are given by

In this case, the order of differentiation and the number of system of equations . Similarly, in each subinterval , we have

where

The system of equations is in terms of the functions, , , and and in order to apply our method, we set , , and . This leads to

The variable coefficients for (82) are given by

for (83) are given by

and for (84) are given by

Equations (79), (80), and (81) can be expressed, respectively, as the following matrix system:

where

The vectors , and are defined as

for .

4. Series Solution for the Limiting Case

In this section, the systems of equations presented above are solved using the series solution approach. These series solutions for the limiting case when is very large are used to validate the results obtained using the multidomain bivariate spectral local linearisation method. Equations of the form (58)-(59) and (75)-(87) can be decomposed to a sequence of coupled linear ordinary differential equations when . The solutions of the resulting systems of equations can be determined using elementary methods for solving differential equations.

Consider the series solution of (58)-(59). For large values of , the dominant terms in (58) are and . Balancing the orders of magnitude of these terms gives . Since , we may infer that is of as . Assuming that has the same order as for large , then . Thus for large values of , we define the following transformations:

The transformation equation (98) reduces (58)-(59) to

where the primes now denote derivatives with respect to . The corresponding boundary conditions are given by

Using the fact that is large, solutions to (99)-(100) are determined in series form as

Substituting (102) into (99)-(101) and then equating the coefficients of like powers of , we obtain the equations for as

subject to the following boundary conditions:

The analytical solutions of the zeroth order equations (103)-(105) are

For , we sequentially solve the following system of ordinary differential equations:

subject to the following boundary conditions

The series solution of (75)-(87) is now considered. For large values of , the dominant terms in (75) are and . Balancing terms in (75)-(87), for , we get

This leads to the following transformations:

Using the transformations (111)-(112) in (75)-(87), we obtain

The corresponding boundary conditions are given by

For large , we seek solutions to (113)-(116) in series form as

Substituting (117) into (113)-(116) and then equating the coefficients of like powers of , we obtain the equations for as

subject to the following boundary conditions

The sequence of equations obtained when is given as

subject to the following boundary conditions:

5. Results and Discussion

In this section, we present and analyze the results obtained using both the multidomain bivariate spectral local linearisation method (MD-BSLLM) and the series solution. The MD-BSLLM was implemented using MATLAB 2013b. The series solutions were obtained using Mathematica’s NDSolve. The parameters to generate all the solutions are given in captions for the tables and graphs. In this article, the physical quantities of interest in the field of fluid dynamics, the skin friction, Sherwood number, and Nusselt number are presented. These physical quantities were obtained using both methods. We present the values of these quantities for different large values of and . The values of these quantities are in excellent agreement for both methods. This in turn implies that the new approach (MD-BSLLM) presented here can be used to solve other nonsimilar boundary layer partial differential equations. The effect of number of intervals in the solution is also presented in this section. Graphs showing the velocity profiles, temperature profiles, and concentration profiles are also presented for various large values of and . These profiles are in excellent agreement with the results presented by Hussain [2].

We first present the solutions of (58) and (59). The physical quantities, skin friction, and the Nusselt number results obtained by solving the governing equations (58)-(59) are given in Table 1. These results were obtained directly by solving the sequences of ordinary differential equations, for the large limiting case, using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 1. Both solutions are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. We observe that an increase in results in a decrease in both the skin friction and the Nusselt number. It is remarkable that this numerical method was able to produce accurate results for large values of with minimum number of grid points. These results were obtained using , , and converged after five iterations.

Table 1: Comparison of the MD-BSLLM solution and the series solution for and , with , , and .

Multi-Domain BSLLM
Series Solution

5
0.2842719
-3.5066677
0.2842719
-3.5066677

10
0.1428115
-7.0008399
0.1428115
-7.0008399

15
0.0952321
-10.5002490
0.09523208
-10.5002490

20
0.0714271
-14.0001050
0.07142714
-14.0001050

25
0.0571424
-17.5000538
0.05714239
-17.5000538

30
0.0476189
-21.0000311
0.04761886
-21.0000311

35
0.0408162
-24.5000196
0.04081624
-24.5000196

40
0.0357142
-28.0000131
0.03571424
-28.0000131

Figure 1 shows the velocity profiles for different values of . An increase in leads to a decrease in the velocity. The velocity profile is parabolic and hence cannot be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. This velocity profile is in excellent agreement with that obtained by Hussain [2]; thus the MD-BSLLM method can be used as a numerical method for solving large parameter partial differential equations.

Figure 1: Velocity profiles at different values of with , , and .

Figure 2 shows the temperature profile for various values of . The temperature profile behaves like an exponential decay curve, that is, an increase in the pseudosimilarity variable results in a decrease in the temperature. Temperature profiles were also obtained for large values. Decreasing leads to an increase in the temperature profile. Figures 1 and 2 were obtained using and .

Figure 2: Temperature profiles at different values of with , , and .

For the second problem, the skin friction, Sherwood number, and Nusselt number results obtained by solving the governing equations (75)-(87) are given in Table 2. Similarly, these results were obtained directly by solving the sequences of ordinary differential equations, for the large , using Mathematica’s DSolve command. The series solutions and the MD-BSLLM solutions are accurate up to at least eight decimal digits as indicated in Table 2. The MD-BSLLM and series solutions of (75)-(87) are in excellent agreement and thus validating the accuracy of the MD-BSLLM method. Increasing the transpiration parameter results in a decrease of the skin friction, Sherwood number, and Nusselt number. These results were obtained using few grid points, and .

Table 2: Comparison of the multidomain solution and the series solution for , , and : , , , and .

Multi-Domain BSLLM
Series Solution

5
0.3088214
3.5018961
3.0018658
0.3088214
3.5018961
3.0018658

10
0.1547399
7.0002370
6.0002332
0.1547399
7.0002370
6.0002332

15
0.1031717
10.5000702
9.0000691
0.1031717
10.5000702
9.0000691

20
0.0773803
14.0000296
12.0000292
0.0773803
14.0000296
12.0000292

25
0.0619045
17.5000152
15.0000149
0.0619045
17.5000152
15.0000149

30
0.0515872
21.0000088
18.0000086
0.0515872
21.0000088
18.0000086

35
0.0442176
24.5000055
21.0000054
0.0442176
24.5000055
21.0000054

40
0.0386905
28.0000037
24.0000036
0.0386905
28.0000037
24.0000036

Table 3 gives a comparison of the results obtained using the MD-BSLLM and the bivariate spectral local linearisation method (BSLLM) for large values of . The main purpose of the results is to confirm that the accuracy of the BSLLM method deteriorates as increases. We observe that the computed values agree up to an average of three digits. This implies that decomposing the main interval into subintervals and solving the nonsimilar partial differential equations in each subinterval improves the accuracy of the BSLLM method. Table 3 in addition displays the computational time of both methods. For , , and , the BSLLM method takes about a fraction of a second to achieve inaccurate results meanwhile the MD-BSLLM method takes about 18 seconds to achieve accurate results.

Table 3: Comparison of the MD-BSLLM and the BSLLM solution for , , and : , , , , , and .

BSLLM
Multi-Domain BSLLM

5
0.3046603
3.5051984
3.0093278
0.3088214
3.5018961
3.0018658

10
0.1546097
7.0003327
6.0003186
0.1547399
7.0002370
6.0002332

15
0.1030889
10.5001437
9.0001125
0.1031717
10.5000702
9.0000691

20
0.0773803
14.0000294
12.0000289
0.0773803
14.0000296
12.0000292

25
0.0619452
17.4999716
14.9999874
0.0619045
17.5000152
15.0000149

30
0.0516135
20.9999793
17.9999897
0.0515872
21.0000088
18.0000086

CPU Time
0.845
0.845
0.845
17.527
17.527
17.527

Table 4 shows a comparison of the BSLLM method and the MD-BSLLM method with different grid points. For the MD-BSLLM method, we use , , and , while for the BSLLM method, we use and to achieve comparable accurate results. Since the BSLLM method requires more grid points to achieve accurate results over a large domain, then it takes more computational time compared to the MD-BSLLM method over the same interval. In Table 4, the MD-BSLLM method takes about seconds to achieve accurate results; meanwhile the BSLLM method takes about seconds to achieve accurate results. This implies that the MD-BSLLM method takes less computational time compared to the BSLLM method to achieve accurate results.

Table 4: Comparison of the MD-BSLLM and the BSLLM solution for , , and : , , , , , and .

BSLLM
Multi-Domain BSLLM

5
0.3088066
3.5018993
3.0018843
0.3088214
3.5018961
3.0018658

10
0.1547399
7.0002370
6.0002332
0.1547399
7.0002370
6.0002332

15
0.1031717
10.5000702
9.0000691
0.1031717
10.5000702
9.0000691

20
0.0773801
14.0000296
12.0000292
0.0773803
14.0000296
12.0000292

25
0.0619014
17.5000152
15.0000149
0.0619045
17.5000152
15.0000149

30
0.0515660
21.0000087
18.0000086
0.0515872
21.0000088
18.0000086

CPU Time
77.921
77.921
77.921
17.527
17.527
17.527

In Table 5, we observe that, for small values of , the results from both methods (BSLLM and MD-BSLLM) are in good agreement. The BSLLM method takes about seconds while the MD-BSLLM method takes about seconds to achieve accurate results. Therefore, we can conclude that it is not necessary to use the MD-BSLLM method to solve the partial differential equations in small subdomains of . The multidomain approach is suitable for solving equations with large values of . The results in Table 5 were obtained with and for the BSLLM method; meanwhile , , and were used for the MD-BSLLM method for .

Table 5: Comparison of the MD-BSLLM and the BSLLM solution for , , and :, , and .

BSLLM
Multi-Domain BSLLM

0.25
0.6346269
0.6230281
0.6046615
0.6346269
0.6230281
0.6046615

0.50
0.6393410
0.7247205
0.6899358
0.6393410
0.7247205
0.6899358

0.75
0.6379544
0.8353765
0.7818927
0.6379544
0.8353765
0.7818927

1.00
0.6307080
0.9546061
0.8803437
0.6307080
0.9546061
0.8803437

1.25
0.6181647
1.0818907
0.9850250
0.6181647
1.0818907
0.9850250

1.50
0.6011394
1.2166019
1.0955986
0.6011394
1.2166019
1.0955986

1.75
0.5806016
1.3580272
1.2116586
0.5806016
1.3580272
1.2116586

2.00
0.5575697
1.5054016
1.3327435
0.5575697
1.5054016
1.3327435

CPU Time
3.088
3.088
3.088
22.436
22.436
22.436

Table 6 shows the effect of the number of subintervals . Table 6 shows that increasing the number of subintervals yields more accurate results. We note that when , that is, we have one domain, then the MD-BSLLM is behaves exactly as the BSLLM method. For , the results do not match because the MD-BSLLM is exactly the BSLLM method. However, for , only ten subintervals ensure accuracy of the MD-BSLLM method as depicted in Table 6. These results confirm that decomposing the main domain into subintervals improves the accuracy of the BSLLM method.

Table 6: Comparison of the multidomain solution and the series solution for , , and : , , , , , , and .

Multi-Domain BSLLM
Series Solution

1
0.0967112
11.2002848
9.6003085
0.0967241
11.2000579
9.6000569

5
0.0967241
11.2000577
9.6000568
0.0967241
11.2000579
9.6000569

10
0.0967241
11.2000579
9.6000569
0.0967241
11.2000579
9.6000569

15
0.0967241
11.2000579
9.6000569
0.0967241
11.2000579
9.6000569

20
0.0967241
11.2000579
9.6000569
0.0967241
11.2000579
9.6000569

Figures 3, 4, and 5 show the velocity, temperature, and concentration profiles for various values of , respectively. It is observed that an increase in leads to a decrease in the velocity, temperature, and concentration. The velocity profile is parabolic and hence difficult to be obtained by certain numerical methods like the Gauss-Seiedel based pseudospectral relaxation method. The temperature and concentration profiles decay exponentially. The velocity, temperature, and concentration profiles are in excellent agreement with the results by Hossain [1]. This implies that the MD-BSLLM method can be used as a numerical method for solving large parameter nonsimilar, nonlinear partial differential equations.

Figure 3: Velocity profiles at different values of , , , , , , and .

Figure 4: Temperature profiles at different values of , , , , , , and .

Figure 5: Concentration profiles at different values of , , , , , , and .

Figures 6, 7, and 8 show the convergence error graphs for the velocity, temperature, and concentration profiles for various values of , respectively. It is observed that the method converges after approximately three iterations to a very small error. This validates the accuracy and convergence of the method.

Figure 6: Convergence error graph of at different values of , , , , , , and .

Figure 7: Convergence error graph of at different values of , , , , , , and .

Figure 8: Convergence error graph of at different values of , , , , , , and .

6. Conclusion

In this article, a new numerical approach for solving systems of nonsimilar nonlinear partial differential equations over large domains is presented and implemented. The method is generalized for any th order system of partial differential equations. The method was tested on two different numerical experiments. The results were validated by using series solution of the numerical experiments. The results of the numerical method and the series solution were in excellent agreement and hence validating the accuracy of the method over large domains. The idea of decomposing the main domain into subintervals increases the accuracy of the method. The more subintervals we have, the more accurate the method is. However, for small domains, it is recommended to use and implement the BSLLM method other than the MD-BSLLM method since the BSLLM method achieves accurate results within few seconds as opposed to the MD-BSLLM method. This project has added to literature a novel approach for solving nonsimilar nonlinear system of partial differential equations over a large domain. In future, we will develop and solve systems of nonlinear partial differential equations arising from other disciplines.

Data Availability

All the data used to support the findings of this study are freely available from the corresponding author upon request.

Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this article.

Acknowledgments

The author would like to acknowledge funding from the DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS). However, the opinions expressed and conclusions arrived at are those of the author and are not necessarily to be attributed to the Centre of Excellence in Mathematical and Statistical Sciences. Part of this work was presented at the 40th South African Symposium of Numerical and Applied Mathematics, held at Stellenbosch University, from 22nd of March till the 24th of March 2016 and at the Conference of the Southern Africa Mathematical Sciences Association 2016, which was held in Pretoria from the 21st to the 24th of November 2016.

References

M. A. Hossain and S. C. Paul, “Free convection from a vertical permeable circular cone with non-uniform surface temperature,” Acta Mechanica, vol. 151, no. 1-2, pp. 103–114, 2001.

S. Hussain, M. A. Hossain, and M. Wilson, “Natural convection flow from a vertical permeable flat plate with variable surface temperature and species concentration,” Engineering Computations, vol. 17, no. 7, pp. 789–812, 2000.

C. D. Pruett and C. L. Streett, “A spectral collocation method for compressible, non-similar boundary layers,” International Journal for Numerical Methods in Fluids, vol. 13, no. 6, pp. 713–737, 1991.

M. Sajid, T. Hayat, and S. Asghar, “Non-similar analytic solution for MHD flow and heat transfer in a third-order fluid over a stretching sheet,” International Journal of Heat and Mass Transfer, vol. 50, no. 9-10, pp. 1723–1736, 2007.

M. Sajid and T. Hayat, “Non-similar series solution for boundary layer flow of a third-order fluid over a stretching sheet,” Applied Mathematics and Computation, vol. 189, no. 2, pp. 1576–1585, 2007.

S.-J. Liao, “A general approach to get series solution of non-similarity boundary-layer flows,” Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 2144–2159, 2009.

A. Chamkha, R. S. R. Gorla, and K. Ghodeswar, “Non-similar solution for natural convective boundary layer flow over a sphere embedded in a porous medium saturated with a nanofluid,” Transport in Porous Media, vol. 86, no. 1, pp. 13–22, 2011.

S. Abbasbandy and T. Hayat, “On series solution for unsteady boundary layer equations in a special third grade fluid,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 8, pp. 3140–3146, 2011.

M. A. Hossain, M. S. Munir, and H. S. Takhar, “Natural convection flow of a viscous fluid about a truncated cone with temperature dependent viscosity,” Acta Mechanica, vol. 140, no. 3, pp. 171–181, 2000.

M. A. Hossain and H. S. Takhar, “Radiation effect on mixed convection along a vertical plate with uniform surface temperature,” Heat and Mass Transfer, vol. 31, no. 4, pp. 243–248, 1996.

M. A. Hossain, M. A. Alim, and D. A. S. Rees, “The effect of radiation on free convection from a porous vertical plate,” International Journal of Heat and Mass Transfer, vol. 42, no. 1, pp. 181–191, 1999.

K. A. Yih, “MHD forced convection flow adjacent to a non-isothermal wedge,” International Communications in Heat and Mass Transfer, vol. 26, no. 6, pp. 819–827, 1999.

V. M. Magagula, S. S. Motsa, P. Sibanda, and P. G. Dlamini, “On a bivariate spectral relaxation method for unsteady magneto-hydrodynamic flow in porous media,” SpringerPlus, vol. 5, no. 1, pp. 1–15, 2016.

S. S. Motsa, “A New Spectral Local Linearization Method for Nonlinear Boundary Layer Flow Problems,” Journal of Applied Mathematics, vol. 2013, Article ID 423628, 15 pages, 2013.

S. Shateyi and S. S. Motsa, “Successive linearization of heat and mass transfer over an unsteady stretching permeable surface in the presence of thermal radiation and a variable chemical reaction,” in Mass Transfer - Advances in Sustainable Energy and Environment Oriented Numerical Modeling, InTech, 2013.

S. S. Motsa and S. Shateyi, “Numerical analysis of mixed convection magnetohydrodynamic heat and mass transfer past a stretching surface in a micro-polar fluid-saturated porous medium under the influence of ohmic heating,” in Mass Transfer - Advances in Sustainable Energy and Environment Oriented Numerical Modeling, p. 10, InTech, 2013.

S. S. Motsa, Z. G. Makukula, and S. Shateyi, “Spectral Local Linearisation Approach for Natural Convection Boundary Layer Flow,” Mathematical Problems in Engineering, vol. 2013, Article ID 765013, 7 pages, 2013.

S. Shateyi and G. T. Marewo, “On a new numerical analysis of the Hall effect on {MHD} flow and heat transfer over an unsteady stretching permeable surface in the presence of thermal radiation and heat source/sink,” Boundary Value Problems, 2014:170, 17 pages, 2014.

S. Shateyi and G. T. Marewo, “Numerical analysis of unsteady {MHD} flow near a stagnation point of a two-dimensional porous body with heat and mass transfer, thermal radiation, and chemical reaction,” Boundary Value Problems, 2014:218, 18 pages, 2014.

S. S. Motsa, “On the practical use of the spectral homotopy analysis method and local linearisation method for unsteady boundary-layer flows caused by an impulsively stretching plate,” Numerical Algorithms, vol. 66, no. 4, pp. 865–883, 2014.

S. S. Motsa, “On the new bivariate local linearisation method for solving coupled partial differential equations in some applications of unsteady fluid flows with heat and mass transfer,” in Mass Transfer - Advancement in Process Modelling, InTech, 2015.

S. S. Motsa and I. L. Animasaun, “A new numerical investigation of some thermo-physical properties on unsteady mhd non-darcian flow past an impulsively started vertical surface,” Thermal Science, no. 19, pp. 249–258, 2015.

S. M. Sandile and I. L. Animasaun, “Paired quasi-linearization analysis of heat transfer in unsteady mixed convection nanofluid containing both nanoparticles and gyrotactic microorganisms due to impulsive motion,” Journal of Heat Transfer, vol. 138.11, Article ID 114503, 2016.

S. S. Motsa and I. L. Animasaun, “Unsteady boundary layer flow over a vertical surface due to impulsive and buoyancy in the presence of thermal-diffusion and diffusion-thermo using bivariate spectral relaxation method,” Journal of Applied Fluid Mechanics, vol. 9, no. 5, pp. 2605–2619, 2016.

S. S. Motsa and I. L. Animasaun, “Bivariate spectral quasi-linearisation exploration of heat transfer in the boundary layer flow of micropolar fluid with strongly concentrated particles over a surface at absolute zero due to impulsive,” International Journal of Computational Science and Mathematics, vol. 9, no. 5, pp. 455–473, 2018.

D. Funaro, “A multidomain spectral approximation of elliptic equations,” Numerical Methods for Partial Differential Equations, vol. 2, no. 3, pp. 187–205, 1986.

P. P. Harald, E. K. Lawrence, A. S. Mark, and A. T. Saul, “A multidomain spectral method for solving elliptic equations,” Computer Physics Communications, vol. 152, no. 3, pp. 253–273, 2003.

D. A. Kopriva, “A spectral multidomain method for the solution of hyperbolic systems,” Applied Numerical Mathematics, vol. 2, no. 3-5, pp. 221–241, 1986.

R. Peyret, “The Chebyshev multidomain approach to stiff problems in fluid mechanics,” Computer Methods Applied Mechanics and Engineering, vol. 80, no. 1-3, pp. 129–145, 1990.

S. S. Motsa, P. Dlamini, and M. Khumalo, “A new multistage spectral relaxation method for solving chaotic initial value systems,” Nonlinear Dynamics, vol. 72, no. 1-2, pp. 265–283, 2013.

S. S. Motsa, P. G. Dlamini, and M. Khumalo, “Solving hyperchaotic systems using the spectral relaxation method,” Abstract and Applied Analysis, vol. 2012, Article ID 203461, 18 pages, 2012.

P. G. Dlamini, M. Khumalo, and S. S. Motsa, “A note on the multi-stage spectral relaxation method for chaos control and synchronization,” International Journal of Nonlinear Sciences and Numerical Simulation, vol. 15, no. 5, pp. 289–298, 2014.

S. N. Hassan and P. Rebelo, “Multistage spectral relaxation method for solving the hyperchaotic complex systems,” The Scientific World Journal, vol. 2014, 2014.

M. S. Sydney et al., “A multi-domain spectral collocation approach for solving lane-emden type equations,” Numerical Simulation-From Brain Imaging to Turbulent Flows, 2016.

M. Shamsi and M. Dehghan, “Recovering a time-dependent coefficient in a parabolic equation from overspecified boundary data using the pseudospectral legendre method,” Numerical Methods for Partial Differential Equations, vol. 23, no. 1, pp. 196–210, 2007.

C. G. Canuto, M. Y. Hussaini, A. M. Quarteroni, and T. A. Zang, Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics (Scientific Computation), Springer, New York, NY, USA, 2007.

R. E. Bellman and R. E. Kalaba, Quasilinearization and Nonlinear Boundary-Value Problems, 1965.