1. Introduction
The simplest example of the application of the boundary layer theory is related to the celebrated Blasius [
1] problem. This problem describes the flow around a very thin flat plate.
The first goal of this paper is to numerically solve, with great accuracy, the Magneto-Hydro Dynamics (MHD) boundary layer equation that governs the flow of an incompressible fluid past a flat plate by an implicit non-standard finite difference scheme that is defined on a quasi-uniform mesh that allows imposing the given boundary conditions exactly.
Numerical methods for problems, like the one considered in this paper, can be classified according to the numerical treatment of the boundary condition imposed at infinity. The oldest and simplest treatment is to replace infinity with a suitable finite value, the so-called truncated boundary. However, being the simplest approach, this has revealed, within the decades, some drawbacks that suggest not to apply it, especially if we have to face a given problem without any clue on its solution behaviour. Several other treatments have been proposed in the literature to overcome the shortcomings of the truncated boundary approach. In this research area, some treatments are worthy of consideration: the formulation of so-called asymptotic boundary conditions by de Hoog and Weiss [
2], Lentini and Keller [
3], and Markowich [
4,
5]; the reformulation of the given problem in a bounded domain, as first studied by de Hoog and Weiss and more recently developed by Kitzhofer et al. [
6]; the free boundary formulation that was proposed by Fazio [
7], where the unknown free boundary can be identified with a truncated boundary; the treatment on the original domain via pseudo-spectral collocation methods, see the book by Boyd [
8] or the review by Shen and Wang [
9] for more details on this topic; and, finally, a non-standard finite difference scheme on a quasi-uniform grid that is defined on the original domain by Fazio and Jannelli [
10]. The proposed non-standard finite difference scheme has been successively modified by Fazio and Jannelli [
11] and used in [
12,
13]. Moreover, this method has been applied to the numerical solution of the perpetual American put option problem of financial markets, see Fazio [
14].
The key advantage of the presented approach is that the given BVP is solved on a semi-infinite interval by introducing a stencil that is constructed in such a way that the boundary conditions at infinity are exactly assigned. In this way, the whole infinite domain is taken into account in the mapping where the grid-points are located at the mid-point of each sub-interval, and then the difficulty that is caused by numerical treatment of the last infinite sub-interval is avoided.
Finally, via a mesh refinement, by using Richardson’s extrapolation, we improve the accuracy of the obtained numerical results, showing that they are in excellent agreement with the solutions found in the literature. This study concludes by comparing the current numerical results with those that are given by the integral approximation method (ITM) and the non-integral technique (NIT) used by Singh and Chandarki [
15].
3. The Finite Difference Scheme
Without a loss of generality, we consider the class of BVPs
where
is a
dimensional vector with
for
as components,
, and
. Here, and in the following, we use Lambert’s notation for the vector components [
17] (pp. 1–5).
In order to solve the problem (
2) on the original domain, we first discuss quasi-uniform grids maps from a reference finite domain and then introduce, on the original domain, a non-standard finite difference scheme that allows for us to impose the given boundary conditions exactly. Let us consider the smooth strict monotone quasi-uniform maps
, the so-called grid generating functions, see Boyd [
8] (pp. 325–326) or Canuto et al. [
18] (p. 96),
and
where
,
, and
is a control parameter. Accordingly, a family of uniform grids
defined on interval
generates one parameter family of quasi-uniform grids
on the interval
. The two maps (
3) and (
4) are referred to as logarithmic and algebraic maps, respectively. As far as the authors’ knowledge is concerned, van de Vooren and Dijkstra [
19] were the first to use these kinds of maps. We notice that more than half of the intervals are in the domain with a length approximately equal to
c and
for (
3), while
for (
4). For both maps, the equivalent mesh in
x is nonuniform, with the most rapid variation occurring with
. The logarithmic map (
3) gives slightly better resolution near
than the algebraic map (
4), while the algebraic map gives much better resolution than the logarithmic map as
. In fact, it is easily verified that
for all
, but
.
The problem under consideration can be discretized by introducing a uniform grid
of
nodes in
with
and
with
, so that
is a quasi-uniform grid in
. The last interval in (
3) and (
4), namely
, is infinite, but the point
is finite, because the non integer nodes are defined by
with
and
. These maps allow for us to describe the infinite domain by a finite number of intervals. The last node of such grid is placed on infinity, so right boundary conditions are correctly taken into account.
We approximate the values of the scalar variable
and its derivative at mid-points of the grid
, for
, while using non-standard difference discretizations
We emphasize that the key advantage of our non-standard finite difference formulation is to overcome the difficulty of the numerical treatment of the boundary conditions at infinity. In fact, the formulae (
5) use the value
, but not
and, then, the boundary conditions at infinity are taken into account in a natural way.
For the class of BVPs (
2), a non-standard finite difference scheme on a quasi-uniform grid can be defined using the approximations given by (
5) above, and it can be written, as follows
where
for
. The finite difference formulation (
6) has an order of accuracy
. It is evident that (
6) is a nonlinear system of
equations in the
unknowns
. For the solution of (
6), we can apply the classical Newton’s method along with the simple termination criterion
where
, for
and
, is the difference between two successive iterate components and
is a fixed tolerance.
4. Numerical Results and Comparison
In this Section, we present the numerical results that were obtained by solving the mathematical model (
1) using the non-standard finite difference scheme (
6) on the quasi-uniform grid that was defined by the logarithmic map (
3) with control parameter
. The control parameter
c defines the grid points distribution in the original physical infinite domain and it is chosen in such a way that
is not too large, being
. The chosen value
leads to obtaining a finer grid points distribution close to
and then a better resolution where the solution has a transient behavior, with an increasing spatial resolution going toward infinity.
Now, let us rewrite the model (
1) as a first-order system, as follows
with
or, in an equivalent form,
where
is a three-dimensional vector with components
for
, and
and
, with
.
For all tests, we consider the control parameter
, a fixed tolerance
, and, as a first guess for Newton’s iteration, the following initial data
In
Table 1, we report the numerical results that were obtained for the missing initial data
for increasing points number
N and
. Here, “iter ”stands for the number of the Newton’s iterations. It is well known that Newton’s method could have a computational cost that is quite high, depending on the initial iterate and the analytic form of the function whose root is being sought. A suitable choice of the initial iterate can contribute to reducing the computational cost. In this context, the chosen initial guess leads to obtain a low iterations number: 6 or 7 iterations, depending on the values of the grid point number.
In
Figure 1, we show the numerical solution that was obtained for
and
. The recovered value of the second order derivative of the solution at the origin is
, obtained in 6 iterations.
Table 2 and
Table 3 list the obtained numerical results for different values of parameter
, with
and
. For the sake of brevity, we have chosen to only report the values of the wall shear stress, which is the second derivative value at the origin. Within the same table, we can compare our results with those that were reported by Singh and Chandarki [
15]. The last column, named
Err Rel, reports the values of the relative errors, as computed by
where, DTM, NIT, and FD represent the numerical solutions obtained, for different values of
, by the integral approximation method (ITM) and the non-integral technique (NIT) used by Singh and Chandarki [
15], and by the proposed numerical method, respectively. In both cases, the value of the relative errors decrease as the values of
increase. Note that the relative errors are all of same the order and, moreover, for
, we obtain the smaller relative error. This is due to the fact that
affects the non-linearity of the model under study. The problem with
corresponding to the Blasius problem and, in this case, the computed missing initial condition
can be compared with the value
that wascomputed by Fazio [
7] by a free boundary formulation of the Blasius problem.
We apply Richardson’s extrapolation, using several refinements of the computational domain, in order to improve the accuracy of the computed solution. On the computational domain of the problem, we build a quasi-uniform grid with a mesh-points number equal to
and proceed with subsequent grid refinements by constructing meshes with grid-point numbers
for
, where
with refinement factor
. On each grid, the numerical solution
,
is computed using the non-standard finite difference method. We adopt a continuation strategy in order to reduce the calculations; in fact, we use the final solution
that was obtained on the grid
g as an initial guess for calculating the solution
on the grid
, where the new grid values are approximated by linear interpolations. We define the level of the Richardson’s extrapolation by the index
k and the two numerical solutions related to the grids
g and
at the extrapolated level
k by
and
. We use the following formula to calculate a more accurate approximation
In
Table 4, we report the extrapolated values with
grid points for
. The last extrapolated value is
and it can be considered as our benchmark value for
. We can conclude that the reported extrapolated value is correct up to 9 decimal places.
5. Concluding Remarks
In this paper, the mathematical model, which describes the MHD boundary layer flow of an incompressible fluid past a flat plate, is solved by the non-standard implicit finite difference method on a quasi-uniform grid for the different magnetic parameters . We solved the given BVP on a semi-infinite interval by introducing a stencil that is built in such a way that the boundary conditions at infinity are exactly assigned. Moreover, we improved the accuracy of the numerical results by Richardson’s extrapolation, which allowed for us to find a more accurate solution.
The values of the second-order derivative of the solution at the origin, for different values of parameter
, are reported in the
Table 2 and
Table 3. The results are compared with those by Singh and Chandarki [
15] in order to verify the accuracy of the proposed method. Moreover, in the case of the Blasius problem, we compare the missing initial condition with the one that was computed by Fazio [
7] by a free boundary formulation of the Blasius problem. The computed values are found to be really accurate.