A large variety of physical phenomena exhibit thin regions with fairly high variation of parameters which describe their behavior, the so called boundary layers. For instance, viscous flows have such layers next to solid surfaces where the tangential velocity increases rapidly from zero. The use of finite differences or of the more efficient Chebyshev collocation in conjunction with successive stretching is expensive when such phenomena are numerically analyzed. Thus one of the main aim of this study is to introduce a relatively low in cost method which improve considerably the resolution of the boundary layers.

From the mathematical point of view we are mainly concerned with the numerical treatment of some singular nonlinear second order boundary value problems attached to ordinary differential equations defined on the half-line. The linearization of such problems around some constant solutions are also considered.

The first challenge of such problems consist in imposing the boundary conditions at infinity. The most frequent remedy for this difficulty is the domain truncation or a sort of a change of coordinates (”mapping”). Up to our knowledge only for the classical Blasius problem there exists a rigorous error analysis related to the truncated boundary formulation of this problem. The result is due to Rubel LAR ( see also RF ). A general and at the same time rigorous discussion on this topic is available in the well known paper of Lentini and Keller LK . The main difficulty with their asymptotic boundary condition (ABC) method consists in the fact that this method, in the nonlinear case, produces nonlinear boundary conditions. Alternatively, the use of Laguerre functions, which means the Laguerre polynomials times $\exp (-x/2)$, assures the exact fulfillment of boundary conditions at infinity.

A second genuine difficulty is the resolution of narrow regions of sharp variations of unknowns. Tang and Trummer in TT observe that very thin boundary layers must have at least one or two collocation points within boundary layer in order to assure reliable outcomes. But this results in extremely fine discretization if the relative spacing of the grid points remains unchanged. Thus they invoke a coordinate stretching technique in conjunction with Chebyshev collocation in order to solve linear second order singular perturbed BVPs. This is a frequently used manner to solve singular perturbation problems. In Laguerre collocation a scaling parameter, which maps the half-line in itself, will assure a suitable distribution of collocation nodes in regions of interest.

Consequently, in order to overcome these difficulties we use in this paper a Laguerre collocation method based on Laguerre functions. Systematic results in solving linear second order boundary value problems and respectively eigenvalue problems, both formulated on a semi-infinite interval, were first obtained by Weideman and Reddy in their seminal paper WR and then by Boyd et al. in JPB . Boyd et al. carries out a detailed comparison between Fourier-sine method, Laguerre series and rational Chebyshev expansion in solving the hydrogen atom problem (Coulomb-Schrödinger equation). They remark that pseudospectral methods based on Laguerre series and rational Chebyshev expansion are better suited for problems in the interval $(0,\mathrm{\infty })$ than Fourier method because they are less sensitive to the variation of the scaling parameter.

The paper is structured as follows. In Section 2, we shortly review the Laguerre collocation successfully used in our previous paper CIG and write down the general form of the nonlinear algebraic system in which this method cast a specific nonlinear singular boundary value problem formulated on the half-line. In Section 3, we reconsider the Blasius problem with respect to the error introduced by domain truncation. In Sections 4 and respectively 5, we introduce and solve two challenging test problems in physics (fluid mechanics and thermodynamics). The first one is a model equation for low Reynolds number flows due to Cohen, Fokas and Lagerstrom. It is of particular concern here because its solution is not a monotone one on the interval of interest (see Fig 1 and Fig 2 in Sect. 4). This lack of monotonicity is an important difference when compared with the Blasius problem or even the second test problem. This second test problem is the density profile equation from hydrodynamics. As we do not have a software control of the achieved accuracy, we choose to validate the obtained numerical results via a mesh refinement, i.e., by repeatedly increasing the cut off parameter in the Laguerre interpolating representation of the solution. Meaningful results are reported for illustrative purposes in Sections 4 and 5 for both problems. Moreover, in Subsection 5.1 we solve a singular Sturm-Liouville problem obtained by linearization around a constant solution. We find pseudospectrum of this problem as well as the first five eigenvectors and comment on the stability of this solution. Section 6 deals with conclusions supported by evidences of the present work.

We are interested in an accurate solution to the following singular nonlinear boundary value problem

where $m$ is a non negative integer, $a$ stands for a set of real parameters and $f$ is a continuous real map. The real parameters $A$ and $B$ in the boundary conditions do not vanish simultaneously.

Up to our knowledge a general theory for the existence (and uniqueness) of solutions to problems (1) is not available. Thus, in spite of the fact that we do not want to elaborate in detail on this topic, we observe that for some nonlinearity $f$ and boundary conditions the problem could be embedded in the following one

where $f:[0,\mathrm{\infty })\times ℝ\times ℝ\to ℝ$ and $p,\frac{1}{q}:[0,\mathrm{\infty })\to [0,\mathrm{\infty })$ are assumed to be continuous. O’Regan in DOR has proved some existence results establishing first the solution existence on a finite interval and then extending the result by Arzela-Ascoli theorem on the semi-infinite interval.

In order to implement the Laguerre collocation for problem (1) we use the differentiation matrices from WR and implement boundary conditions making use of the procedure of removing technique of independent boundary conditions introduced by Hoepffner in JH . This strategy was successfully used in our previous paper CIG and is detailed there.

Thus we cast (1) into the nonlinear algebraic system

where the vector $𝐱$ contains the nodes $x_2,x_3,\mathrm{\dots },x_N$ which are the roots of the Laguerre polynomial $L_{N-1}\left(x\right)$ of order $N-1$, indexed in increasing order of magnitude. It is well known that $x_N=O\left(N\right)$ as $N\to \mathrm{\infty }.$ These nodes are kept and the first node $x_1=0$ is removed. We remove only one node because only one of the parameters $A$ and $B$ do not vanish in the problems which follow. The interval $[0,\mathrm{\infty })$ can be mapped to itself by change of variables $\eta =b\tilde{\eta }$, where $b$ is any positive real number (the scaling factor). Laguerre method therefore contains a free-parameter which will play a key role in our analysis.

The vector $𝐔$ contains the unknown values of $u$ in the nodes $x_2,x_3,\mathrm{\dots },x_N$ and ${\widetilde{𝐃}}^{\left(i\right)},i=1,2$ are the first and the second order differentiation matrices with boundary conditions enforced. As we use MATLAB codes and symbols, $\mathrm{diag}(\frac{\mathrm{m}}{𝐱})$ stands for the diagonal matrix with $m/x_i,i=2,3,\mathrm{\dots },N.$ Throughout this paper the nonlinear were solved by MATLAB built in function fsolve. The initial guess was chosen in each and every case in order to mimic the correct solution and we have obtained convergent procedures.

The problems of type (1) are particularly difficult because of the singularity $1/x$ in the coefficient of the first order differentiation term. To cope with this difficulty we varied the scaling factor such that a fine enough grid will resolve the solution around singularity.

When the independent variable is denoted by $\eta$ and the dependent one by $f\left(\eta \right)$ this problem reads

In the truncated boundary approach strategy, the boundary condition at infinity is replaced by the same condition at a given finite value $L$. Thus we define $f_L\left(\eta \right)$ as the solution of

The error $e\left(\eta \right)$ related to $f_L\left(\eta \right)$ is defined by

Rubel in LAR found out the following upper bound for this error:

It means that Laguerre collocation method used to solve this problem in our paper CIG has produced errors of order $O\left({10}^{-3}\right)$(see also Fig. 4 in this paper). Roughly $L{\left[f_L\left(L\right)\right]}^{-1}$ approaches unity as $\eta \to L$ and thus the r. h. s. in (7) depends essentially on the drag at large distances away from fixed boundary.

We have to remark that in the actual computed value of $f_L\left(\eta \right)$, and its first two derivatives in $L$, the truncation and round off errors of the collocation are also incorporated. However, as we were not aware of the estimation (7) at the time the paper CIG was edited, this discussion validate once more the Laguerre collocation results.

The following problem was introduced by one of the authors of CFL , namely Lagerstrom (1961)

The differential equation in (8) is the equation for the time independent temperature distribution in an infinite medium. Generally, it is somewhat an unrealistic physical model where $x$ is a radial coordinate in $\left(k+1\right)$-dimensional space and $u$ is the temperature. The first two terms come from Laplace’s operator, the last two represent nonlinear autonomous heat sources unfortunately not validated in any actual physical model. The temperature equals zero on the sphere $x=\epsilon$ and equals $U$ at large distances. The intuitive reasoning indicated that for one-dimensional case $\left(k=0\right)$ the problem is not singular. This was verified by the construction of an exact solution in this case (by quadrature and inversion of a function). If $\alpha >0$ and $\beta >0$ one may make their values unity (even simultaneously) by a scale transformation in $u$ and $x$. However, because $\beta =0$ is an interesting case and $\alpha =0$ occurs in an auxiliary equation, i.e., Stokes equation, we consider arbitrary nonnegative values for $\alpha$ and $\beta$. For $\alpha >0$, (8) has a unique solution. Our numerical approximation for such a solution when $\epsilon =0.1$ and $U=1$ is depicted in Fig. 1.

Singular asymptotic techniques may and have been used only if $k\ge 1$. For $\alpha =0$ the problem may be solved explicitly. The boundary value problem then has a solution which is unique only if $k>1.$ Thus $k=1$ is an important limiting case.

Hsiao in H gave a rigorous discussion of the existence of a solution to (8) for $\beta =0,k=1,$ and $\epsilon \to 0^{+}.$ Such a solution is depicted in Fig. 2 for the same value of $\epsilon ,U,N,$and $b$ as stated above. The case $\beta =0,k=1$ has the advantage of being applicable to ”real” problems, i.e., to some problems in fluid dynamics for which (8) is a model.

Both curves depicted in Fig. 1 and Fig. 2 behaves in perfect accordance with the asymptotic for $x\to \mathrm{\infty }$ predicted by Markowich in PAM . There is still a difference between them for large $x$. It seems that the nonlinearity $\beta {\left(\frac{du}{dx}\right)}^2$ affects to some extent the stability of numerical process but ill conditioning of the differentiation matrices does not appear to be a serious problem as our experiments indicate.

We investigate the monotonously increasing solutions of the problem

where $\rho \left(r\right)$ stands for the density of a fluid, and for the sake of simplicity we can choose the following values of the parameters: $\lambda =1$ without restriction of generality, $N=3$ which corresponds to a physically meaningful case and $\xi$ varying in the range $[0,1]$ such as to reflect different physical situations. The above equation is called the density profile equation and has the origins in the Cahn-Hillard theory which is used in hydrodynamics to study the behavior of non-homogeneous fluids. In their paper KKLW , Kitzhofer, Koch, Lima and Weinmūller find by polynomial collocation the so called ”bubble-type solution”. When such solution exists it has exactly one zero $R$ in $(0,\mathrm{\infty }),$ and $R$ is interpreted as the bubble radius. They also satisfy and The derivative of the solution attains a maximum at some value and tends to $0$ at infinity. Eventually, a bubble-type solution exhibits an interior layer which becomes sharper as $\xi \to 1.$ We also have to observe that the problem (9) always admits the constant solution $\rho \left(r\right)=\xi$ which physically corresponds to the case of a homogeneous fluid, i.e., without bubbles.

In KKLW the authors transform the problem to a finite interval and show that essential singular problem obtained is well-posed. In contrast we solve (9) by Laguerre collocation introduced in the previous sections. In fact, with respect to singularly perturbed problems, in TT , Tang and Trummer observe that ”for good resolution of the numerical solution at least one of the collocation points ought to lie in the boundary layer”. Thus, for $N=96$ and scaling factor $b=10$ we get the repartition of the Laguerre collocation points plotted in Fig. 3. For the singular problem at hand it seems fairly well suited.

The ration between the largest step size and the smallest (the closest to $0$) equals $1.2790e+03.$ With such a density of collocation nodes we obtained first the solution of (9) for $\xi =0.$ It is depicted in Fig. 4.

For $\xi =0.1,\mathrm{\dots },1.0$ the solutions to the same problem are displayed in Fig. 5. The values of unknown function $\rho$ in $0,$ the radius $R$ and elapsed CPU time when the nonlinear algebraic system involved by Laguerre collocation is solved by MATLAB code fsolve are reported in Tab. 1. All computations were carried out using 2010a variant of MATLAB on an HPxw8400 workstation with clock speed of 3.2Ghz.

Thus, unlike the paper KKLW we solve the problem (9) for extreme values of $\xi$, i.e., $0$ and $1.$ The solutions displayed in Fig. 5 reproduce correctly the similar solutions from KKLW . The computed values of $\rho \left(0\right)$ and $R$ from Tab. 1 are fairly closed to the corresponding values displayed in the above quoted paper. It is also worth nothing that the computational effort increase with parameter $\xi$ as is apparent from the last column of Tab. 1. The same difficulty is reported in KKLW . Moreover, for $\xi =0.9$ and $\xi =1.0$ even the continuation could hinder fsolve to find a converging solution. Thus, a particular special attention must be payed for the initial guess in these cases. It is also clear from the last two lines of Tab. 1 that rounding off errors contaminate the solution for large $\xi .$

As a final word we have to remark that our numerical experiments carried out with respect to the problem (9) and based on mapping coupled with Chebyshev collocation provided totally unsatisfactory results.

All in all, the Laguerre collocation has incorporated without any difficulty and avoiding domain truncation the boundary conditions at infinity.

It is largely accepted that in order to obtain a fine resolution for fairly sharp layer problems, at least some of the grid points should lie in the boundary layer no matter how narrow the boundary layer is. In the Laguerre collocation method used above this desiderate is accomplished by making use of the scaling parameter $b$. Its tuning is carried out manually in a reduced number of numerical experiments. However, this could be considered as a drawback of our strategy when compared with the automatic mesh adaptation for sharp layers proposed in KKLW . Fortunately, we have reported results for sharpest layer which were not previously obtained. We also believe that numerical solutions for Cohen, Fokas and Lagerstrom type equation are for the first time reported in this paper.

Anyway, an automatic generation of a suitable scaling parameter, for a specific problem, remains an open issue.