The main aim of this paper is to provide accurate spectral collocation solutions to systems of boundary value problems of the following form:

Generally, both equations are genuinely nonlinear and the unknowns $u$ and $v$ as well as both parameters $\lambda$ and $\mu$ have mechanical significance.

The parameters $a$ and $b$ can not be simultaneously zero and so a second order equation in system is enabled.

There are some examples of important technological problems in fluid and heat and/or mass transfer which reduces to systems of type (1). For instance, in [20], Magyari and Keller study the heat and mass transfer problems in the boundary layers on continuous stretching surfaces with a given temperature or heat-flux distribution, moving in an otherwise quiescent fluid medium. The effect of non-linear Rosseland thermal radiation in boundary layer is studied by Pantocratoras and Fang in [25] and the list of papers considering similar physical phenomenon can continue. We restrict ourselves to quote only these papers and observe that typically, in the system (1) the first equation is the Falkner-Skan one which is coupled with a second one modelling a sort of heat transfer.

Exact solutions for self-similar boundary-layer flows induced by permeable stretching walls have been provided by Magyari and Keller in [21] and some similar results obtained by HAM have been reported by Liao in [19]. In both papers the authors solve only Falkner-Skan type problems.

From the numerical point of view some variants of shooting method along with domain truncation have been usually used in order to solve such systems of boundary value problems. The method is based on the use of certain empiric asymptotic initial conditions in order to avoid the boundary condition at infinity. High order equations are transformed into first order systems by the compound matrix method. Some iterative and non iterative transformations, which reduce the problems to solutions of a sequence of Cauchy’s related ones, have also been considered (see for instance [7]).

In his well known monograph [3], Boyd observes that the unboundedness of the computational domain extract two penalties. First, the rate of convergence is sub geometric rather than geometric and second, in addition to the truncation $N$ one always has to tackle some kind of scaling (free) parameter. Moreover, the same author notes in [5] that generally speaking singularities degrade the usual exponential rate of convergence of spectral methods to a finite rate of convergence. That is to say that the error falls proportional with $1/N^{k}$ for some natural $k$, the so called "algebraic order of convergence", where $k$ depends on the type of singularity. Fortunately, for particular classes of problems we can tune some collocation methods in order to recover the exponential rate of convergence.

In some of our recent papers [11, 12] we have successfully used the LC based on the roots of Laguerre polynomials, i.e., the Laguerre Gauss nodes, in order to avoid the well known drawback of shooting or domain truncation. The method proved to be fairly efficient in imposing the boundary conditions in origin and at infinity, and easy implementable. The boundary conditions at infinity are directly satisfied by Laguerre functions and the conditions at zero are introduced by a removing technique of independent boundary conditions. The free (scaling) parameter involved in the these algorithm is adjusted in order to improve the convergence of the Newton type processes which solve the nonlinear algebraic systems generated by collocation discretization. One important drawback of this method is the lack of a polynomial transformation between the physical and frequency (phase) spaces.

Fortunately for mapped ChC and LGRC methods we are in the possession of such transforms. Thus, a good heuristic (empirical), but otherwise general strategy, to optimize the scaling parameter is to simply solve the problem and analyze the spectral coefficients of the expansion in the phase space, plotted as absolute values on a log-linear plot, for several different values of the scaling parameter. This strategy has been introduced by Boyd and his coworkers in [4] (see also [16]) and will be exploited below. It enables a clear separation of truncation and roundoff errors as well as a rough estimation of the global accuracy.

The paper is organized as follows. In Section 2 we briefly review the involved spectral collocation methods and in the Section 2.1 we consider a second order linear two-point boundary value problem in order to make clear the idea of scaling parameter optimization. A third order test problem is analyzed in the Section 2.2. In Section 3 we solve the system (1) composed by Falkner-Skan equation and various heat equations. In Section 4 we treat a problem simulating a viscous incompressible flow over a rotating disk and in Section 5 the Blasius viscous flow of a kind of non-Newtonian fluid is analyzed. We have introduced this problem in the present study because it is really a challenging problem for a numerical technique. Some concluding remarks are displayed in Section 5.

For the LC method, including the differentiation matrices, we refer to the paper of Weideman and Reddy [30]. The method is based on the roots of the Laguerre polynomials. In order to solve the two-point boundary value problem (1) by LC we represent each unknown $u(x)$ by the interpolate $p_{N-1}(x)$ defined by

where $\phi_{j}(x)$ are the Lagrangian interpolating polynomials in barycentric form

The interpolating nodes $x_{j},j=2,\ldots,N$ are the first $N-1$ roots of Laguerre polynomial $L_{N-1}(x)$ of degree $N-1$ indexed in increased order of magnitude, i.e. the Laguerre Gauss system. We add to this system the node $x_{1}=0$ only in order to facilitate the incorporation of boundary conditions. We recall that $x_{N}=O(N)$ as $N\rightarrow\infty$.

The interval $[0,\infty)$ can be mapped to itself by change of variables $x=b\widetilde{x}$, where $b$ is any positive real number (the scaling factor). The LC method therefore contains a free-parameter. It means that the Laguerre interpolating process is exact for functions of the form

where $p(x)$ is any polynomial of degree $N-1$ or less.

The usual Laguerre interpolation was systematically analyzed in some well known textbooks (see for instance [2,14]). However, the most important result remains that provided in [2] which reads

where $I_{N}$ is the usual interpolation operator and the weights $\omega$ are defined as $\omega_{0}:=e^{-x}$ and $\omega_{\tau}:=e^{-(1-\tau)x},0<\tau<1$.

The LGRC is thoroughly analyzed in the monograph of Shen et al. [28] as well as in and [1]. This method uses the so called Laguerre-Gauss-Radau nodes which are the roots of the derivatives of some generalized Laguerre polynomials plus the origin. One of the most important advantage of the latter method is the possibility to formulate a rigorous Fourier-Laguerre analysis. The direct and inverse Laguerre polynomial transforms will be fairly useful below. Up to our knowledge the performances of these methods which use the same Laguerre functions but different sets of nodes, have not been compared. Moreover, in [30] the authors state that "to the best of our knowledge it has not been investigated whether there is a significant advantage to be gained by considering this more general basis set", i.e., the set of generalized Laguerre polynomials.

A partial answer to this dilemma is provided in this paper.
However, in the monograph [28], Ch. 7, Unbounded Domains, a fairly analogue estimation of (5) concerning LGRC method is provided. It additionally contains a term $(\ln(N))^{1/2}$ times a weighted norm of the highest derivative of $u$. The function $u$ is supposed continuous on the positive half line and decaying at infinity along with its derivatives.

The fundamentals of the ChC method are available in the well known monograph of Gottlieb and Orszag [14]. J. P. Boyd in his well known monograph [3], Ch. 17 Methods for Unbounded Domains, analyses methods based on the algebraically mapped Chebyshev polynomials. Our experience with this method is available for instance in the textbook [8] and papers [6, 9, 10, 26].

For the ChC mapped method we use the algebraic map

and its inverse. In order to write down any third order differential equation in independent variable $t$ as an equation in new variable $x$, we need the following derivatives:

The corresponding Chebyshev differentiation matrices for this method come from the same paper [30]. It is clear from the beginning that the ChC mapped method adds supplementary algebraic manipulations and possibly additional roundoff errors.

The real positive number $c$ is the scaling (free) parameter of the mapped ChC method. This parameter along with the scaling factor $b$ of the LC methods can be adjusted with respect to physical and/or mathematical considerations.

As a rule, in the first situation, this parameter can be chosen in order to equal $1/L$ where $L$ stands for a typical length scale of a specific problem. Consequently, in boundary layer type problems $L$ can be the thickness of the layer. This fact is illustrated in Section 2.2, Fig. 3.

The scaling parameter can also be adjusted on a mathematical basis. Thus it can be tuned in order to improve the convergence speed of an iterative process of solving nonlinear algebraic systems or in order to improve the smoothness and the decaying rate of the coefficients of a spectral expansion.

The well known Falkner-Skan problem reads

Physically relevant solutions exist only for $-0.19884<\beta\leq 2$. The linear equation of the heat transfer reads

where $\operatorname{Pr}$ stands for the Prandtl number. The nonlinear decoupled system of boundary value problems (12)-(13) has been solved by Rosales-Vera and Valencia in their paper [27] using a Fourier series procedure. Our numerical results for Falkner-Skan problem are reported in Fig. 5 and are in perfect accordance with those reported in [27].

Pantokratoras and Fang considered in their paper [25] a much more involved heat equation than (13), namely

where $N_{R}$ is the radiation number and $\theta_{r}$ is the temperature parameter. They attache the same boundary condition as in (13) and solve the system (12)-(14) with $\beta=0$ by a shooting method based on a Runge-Kutta scheme.

Instead, we solve the same problem with non vanishing $\beta$ by LGRC. Thus, in the left panel of Fig. 6 we report a complete solution, i.e., $f,df/d\eta$, computed $d^{2}f/d\eta^{2}$ and spline smoothed $d^{2}f/d\eta^{2}$ of the Falkner-Skan problem for $\beta=-0.18$. In the right panel we report the solution of the so called Rosseland thermal radiation (14) for various values of the temperature parameter. In the same line with the results from [25] we observe that as $\theta_{r}$ increases the temperature profiles become broader and S -shaped, that is, an inflexion point appears.

The boundary conditions in the origin are introduced by the removing technique of independent boundary conditions introduced by Hoepffner in [15] and exploited in our previous papers [11-13].

In the monographs [23,24] the following problem is suggested

It simulates the flow of a viscous incompressible fluid over a rotating disk. The solutions $f,g$ and $h$ are respectively proportional with the velocity components $u,v$ and $w$. If $\gamma:=1-\varepsilon$, where $0<\varepsilon\ll 1$, one could show that to order $\varepsilon$ the solutions $f(\zeta),g(\zeta)$ and $h(\zeta)$ are given by

We solve the problem (15) by LGRC and confirm the above asymptotic results (16). Our numerical outcomes are displayed in Fig. 7.

In his paper [17], S-J Liao considers that the following problem is a challenging nonlinear one for numerical techniques. It reads

The two-dimensional Blasius viscous flow of power-law fluid on a semi-infinite flat plane can be described by this equation. Our numerical results are displayed in Fig. 8. We have to observe that for $n=1$ the above problem reduces to the classical

Blasius problem. The leftmost curve in Fig. 8 confirms this fact. Up to our knowledge there are no numerical solutions reported for this problem. However, our numerical solution for $n=2$ looks fairly similar with the solutions reported in [17].

The Table 1 containing the output of fsolve (trust-region dogleg algorithm) when we deal with problem (17) shows that this problem is indeed the most challenging with respect to the set of problems considered so far. We have used a continuation with respect to the parameter $n$ and at the value $n=1.50$ the computational effort to solve the nonlinear algebraic system is the most important (see. For $n=1$ it is comparable with that reported in Section 2.2. The exit flag has been equal to 1 and the residual has situated under a value of order $O\left(10^{-20}\right)$ for each and every value of $n$.

It is important to notice that our tentative to solve this problem using LC failed. It has produced non smoothed curves irrespective of the chosen scaling factor and the order of approximation. Moreover, the nonlinearities in the problem (17) in conjunction with the derivatives from (7) has made the ChC method rather inapplicable.

We have compared in this paper three spectral collocation methods in solving third order boundary value problems on the half line. The most reliable and at the same time the most accurate proved to be LGRC.

This method has an important advantage over LC. A Fourier Laguerre type analysis being available a fairly useful Laguerre polynomial transformation becomes applicable. With this transformation the numerical solution can be transferred from the physical space in the phase space where the coefficients of the expansion can be quantified and analyzed. Thus, a new way to find an optimal scaling parameter is established. The method has provided solutions fairly comparable, at least qualitatively, with those obtained by shooting or HAM.

Moreover, it seems that we have solved numerically for the first time a challenging problem from non-Newtonian fluid mechanics.

It is also fairly important to underline that we have solved each nonlinear problem without the empiricism involved in domain truncation coupled with shooting methods.

Acknowledgments The author is indebted to Dr. Wenjie Liu (NTU Singapore) for some helpful messages concerning the LGRC differentiation matrices, to both anonymous referees for useful comments and remarks and to the editor for his efficiency.

1. 1.

   Benacchio, T., Bonaventura, L.: Absorbing boundary conditions: a spectral collocation approach. Int. J. Numer. Meth. Fluids (2013). doi:10.1002/fld. 3768

   Report issue for preceding element
2. 2.

   Bernardy, C., Maday, Y. In: Ciarlet, P., Lions, J.L. (eds.): Spectral methods, vol. 5 (Part 2). NorthHolland (1997)

   Report issue for preceding element
3. 3.

   Boyd, J.P.: Chebyshev and fourier spectral methods. Dover Publications, New-York (2000)

   Report issue for preceding element
4. 4.

   Boyd, J.P., Rangan, C., Bucksbaum, P.H.: Pseudospectral methods on a semi-infinite interval with applications to the hydrogen atom: a comparison of the mapped Fourier-sine method with Laguerre series and rational Chebyshev expansions. J. Comput. Phys. 188, 56-74 (2003)

   Report issue for preceding element
5. 5.

   Boyd, J.P.: Chebyshev spectral methods and the lane-emden problem. Numer. Math. Theor. Meth. Appl. 4, 142-157 (2011)

   Report issue for preceding element
6. 6.

   Dragomirescu, I.F., Gheorghiu, C.I.: Analytical and numerical solutions to an electrohydrodynamic stability problem. Appl. Math. Comput. 216, 3718-3727 (2010). doi:10.1016/j.amc.2010.05.028

   Report issue for preceding element
7. 7.

   Fazio, R.: A novel approach to the numerical solution of boundary value problems on infinite intervals. SIAM J. Numer. Anal. 33, 1473-1483 (1996)

   Report issue for preceding element
8. 8.

   Gheorghiu, C.I.: Spectral methods for differential problems. Casa Cartii de Stiinta Publishing House, Cluj-Napoca, Romania (2007)

   Report issue for preceding element
9. 9.

   Gheorghiu, C.I., Dragomirescu, I.F.: Spectral methods in linear stability. Application to thermal convection with variable gravity field. Appl. Numer. Math 59, 1290-1302 (2009)

   Report issue for preceding element
10. 10.

    Gheorghiu, C.I., Hochstenbach, M.E., Plestenjak, B., Rommes, J.: Spectral collocation solutions to multiparameter Mathieu’s system. Appl. Math. Comput 218, 11990-12000 (2012). doi:10.1016/j.acm.2012.05.068

    Report issue for preceding element
11. 11.

    Gheorghiu, C.I.: Laguerre collocation solutions to boundary layer type problems. Numer. Algor 64, 385-401 (2013). doi:10.1007/s11075-012-9670-y

    Report issue for preceding element
12. 12.

    Gheorghiu, C.I.: Pseudospectral solutions to some singular nonlinear BVPs. Applications in nonlinear mechanics. Numer. Algor. 68, 1-14 (2015). doi:10.1007/s11075-014-9834-z

    Report issue for preceding element
13. 13.

    Gheorghiu, C.I.: Spectral Methods for Non-Standard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond. Springer Cham Heidelberg New York Dordrecht London (2014)

    Report issue for preceding element
14. 14.

    Gottlieb, D., Orszag, S.A.: Numerical analysis of spectral methods: theory and applications. SIAM, Philadelphia (1977)

    Report issue for preceding element
15. 15.

    Hoepffner, J.: Implementation of boundary conditions. http://www.lmm.jussieu.fr/˜hoepffner/ research/realizing.pdf, Accessed 2 Feb 2015 (2010)

    Report issue for preceding element
16. 16.

    Iacono, R., Boyd, J.P.: The Kidder Equation: $u_{xx}+2xu_{x}/\sqrt{1-\alpha u}=0$. Stud. Appl. Math. 135, 6385 (2014)

    Report issue for preceding element
17. 17.

    Liao, S.-J.: A challenging nonlinear problem for numerical techniques. J. Comput. Appl. Math. 181, 467-472 (2005)

    Report issue for preceding element
18. 18.

    Liao, S.-J.: A new branch of solutions of boundary-layer flows over an impermeable stretched plate. Int. J. Heat Mass Transfer 48, 2529-2539 (2005)

    Report issue for preceding element
19. 19.

    Liao, S.-J.: A new branch of solutions of boundary-layer flows over a permeable stretching plate. Int. J. Nonlinear Mech. 42, 819-830 (2007)

    Report issue for preceding element
20. 20.

    Magyari, E., Keller, B.: Heat transfer characteristics of the separation boundary flow induced by a continuous stretching surface. J. Phys. D: Appl. Phys. 32, 2876-2881 (1999)

    Report issue for preceding element
21. 21.

    Magyari, E., Keller, B.: Exact solutions for self-similar boundary-layer flows induced by permeable stretching walls. Eur. J. Mech. B - Fluids 19, 109-122 (2000)

    Report issue for preceding element
22. 22.

    Mandelzweig, V.B., Tabakin, F.: Quasilinearization approach to nonlinear problems in physics with applications to nonlinear ODEs. Comput. Phys. Commun. 141, 268-281 (2001)

    Report issue for preceding element
23. 23.

    Ockendon, H., Ockendon, J.R.: Viscous flow. Cambridge University Press, Cambridge (1995)

    Report issue for preceding element
24. 24.

    O’Neil, M.E., Chorlton, F.: Viscous and compressible fluid dynamics, p. 395. Wiley, New York (1989)

    Report issue for preceding element
25. 25.

    Pantokratoras, A., Fang, T.: Blasius flow with non-linear Rosseland thermal radiation. Meccanica 49, 1539-1545 (2014)

    Report issue for preceding element
26. 26.

    Plestenjak, B., Gheorghiu, C.I., Hochstenbach, M.E.: Spectral collocation for multiparameter eigenvalue problems arising from separable boundary value problems. J. Comput. Phys. 298, 585-601 (2015). doi:10.1016/j.jcp.2015.06.015

    Report issue for preceding element
27. 27.

    Rosales-Vera, M., Valencia, A.: Solutions of Falkner-Skan equation with heat transfer by Fourier series. Int. Commun. Heat Mass 37, 761-765 (2010)

    Report issue for preceding element
28. 28.

    Shen, J., Tang, T., Wang, L.-L.: Spectral methods. Algorithms, analysis and applications. Springer, Berlin (2011)

    Report issue for preceding element
29. 29.

    Wang, L.-L.: Discrete transform of Laguerre function approach. http://www.ntu.edu.sg/home/lilian/ book.htm, Accessed 5 May 2014 (2011)

    Report issue for preceding element
30. 30.

    Weideman, J.A.C., Reddy, S.C.: A MATLAB differentiation matrix suite. ACM Trans. Math. Software 26, 465-519 (2000)

    Report issue for preceding element
31. 31.

    von Winckel, G.: Fast Chebyshev Transform. http://www.mathworks.com/matlabcentral/ fileexchange/4591-fast-chebyshev-transform– 1d-, Accessed 15 May 2015 (2004)

    Report issue for preceding element