Existence and Uniqueness of BVPs Defined on Semi-Infinite Intervals: Insight from the Iterative Transformation Method

This work is concerned with the existence and uniqueness of boundary value problems defined on semi-infinite intervals. These kinds of problems seldom admit exactly known solutions and, therefore, the theoretical information on their well-posedness is essential before attempting to derive an approximate solution by analytical or numerical means. Our utmost contribution in this context is the definition of a numerical test for investigating the existence and uniqueness of solutions of boundary problems defined on semi-infinite intervals. The main result is given by a theorem relating the existence and uniqueness question to the number of real zeros of a function implicitly defined within the formulation of the iterative transformation method. As a consequence, we can investigate the existence and uniqueness of solutions by studying the behaviour of that function. Within such a context the numerical test is illustrated by two examples where we find meaningful numerical results.


Introduction
We consider the numerical solution of the boundary value problems (BVPs) defined on semi-infinite intervals belonging to the class of problems where (⋅, ⋅, ⋅, ⋅) satisfies appropriate smoothness conditions, 0 , 0 and ∞ are given constants. BVPs belonging to (1.1) arise mainly in boundary layer theory, see Schlichting and Gersten [43], but also in fields of current interest like biology, chemistry, engineering, mechanics or physics.
Before considering the application of any numerical method a question naturally arises: Is the BVP (1.1) well-posed? Indeed, few results are known about the existence and uniqueness question. On this subject we can quote a paper by Weyl [49] about the celebrated Blasius problem [6] and several papers concerning the more general Falkner-Skan model [11]. Hence, any new idea for investigating that question is worth to be considered.
In this paper, we prove a theorem relating the existence and uniqueness question to the number of real zeros of a "transformation function" defined within the formulation of the iterative transformation method (ITM). Then, as a numerical test for the existence and uniqueness of solutions of (1.1) we propose the application of the ITM to look for the number of real zeros of that function. In this context, the first application of that test was considered in [14].
The ITM discussed herein allows us to integrate numerically the class of BVPs (1.1). The method in point allows us to transform a BVP to a sequence of initial value problems. Such a transformation has also theoretical interest because it may represent an intermediate step to prove existence and uniqueness theorems (see [36, pp. 7-13]).
A first application of the ITM to a simple problem describing a biological model was given, on an intuitive basis, in [28]. Moreover, other applications arose in connection with: the length determination for a tubular chemical reactor [14], the shock front propagation in rate-type materials [16], the classical Stefan's problem [17], and the spreading of a viscous fluid under the influence of gravity [17]. A further application of the ITM to the numerical solution of BVPs on infinite intervals governed by a third-order ordinary differential equation was described in [20]. A recent review concerning the applications of the ITM can be found, by the interested reader, in Fazio [26].
It was proved in [19] that the ITM is an extension of the non-ITM proposed in [28]. The ITM can also be applied to two-point BVPs and in this case, it may be seen as an extension of the non-ITM discussed in [12]. However, the applicability of non-ITMs is based on the invariance properties of the governing differential equation and of the given boundary conditions. As a consequence, non-ITMs are considered as ad hoc methods (see [33], [38, pp. 35-36], [39, p. 137], [40, p. 218]). For the interested reader, a recent review reporting on the applications of the non-ITM can be found in Fazio [27].
To the best of our knowledge the ITM is the first transformation method applicable to a general class of problems.
The remainder of the paper is organized as follows. In section 2 we define the ITM. In section 3 we provide a simple constructive proof of the main result that relates the existence and uniqueness of solutions of BVPs defined on infi-nite intervals to the number of real zeros of the transformation function. There, we point out how the present approach reduces the question of existence and uniqueness to find the number of real zeros of the transformation function. A similar theorem has been proved by this author for a class of free BVPs governed by second-order differential equations, see [21]. Some guidelines on the definition of the transformation function and a discussion on its sensitivity analysis are developed in section 4. In section 5 we apply the iterative method to two BVPs describing, respectively, the Sakiadis problem and the Falkner-Skan model.
Our main concern is to verify numerically whether the considered problem has a unique solution. In this context, the test indicates the existence and uniqueness of the solution for the first model and nonuniqueness for the second one. The last section is devoted to the final remarks and conclusions. On the basis of the theorem proved in section 3, the proposed numerical test represents a possible way to investigate the existence and uniqueness of solution of BVPs. Further evidence for this conclusion is given by the numerical experiments proposed in section 5 as well as by the study reported in [14].
In the following sections, we apply the ITM to the Sakiadis problem [41,42] where and are appropriate similarity variables and is a parameter. This set is called the Falkner-Skan model, after the names of two English mathematicians who first studied it [11]. Let us notice that for = 0 the BVP (1.3) reduces to the celebrated Blasius problem [6]. As pointed out by Na [39, pp.  The existence and uniqueness question for the problem (1.3) is really a complex matter. Assuming that > 0 and under the restriction 0 < < 1, known as normal flow condition, Hartree [35] and Stewartson [45] proved that the problem (1.3) has a unique solution, whose first derivative tends to one exponentially.
Coppel [8] and Craven and Peletier [9] pointed out that the above restriction on the first derivative can be omitted when 0 ≤ ≤ 1. For each value of the parameter , there exists a physical solution with positive monotone decreasing second derivative, in the domain [0, ∞), that approaches zero as the independent variables goes to infinity, as proved by Weyl [49]. In the case > 1, the Falkner-Skan model loses the uniqueness and a hierarchy of solutions with reverse flow exists.
In fact, for > 1 Craven and Peletier [10] computed solutions for which < 0 for some value of . In each of these solutions, the velocity approaches its limit exponentially in . Here and in the following the term normal flow indicates that the flow velocity has a unique direction, and instead, reverse flow means that the velocity is both positive and negative in the integration interval.
The considered problem has also multiple solutions for min < < 0, as reported by Veldman and van de Vooren [48], with the minimum value of given Hartree [35] was the first to solve numerically the Falkner-Skan model. Then, Cebeci and Keller [7] applied shooting and parallel shooting methods requiring asymptotic boundary condition to be imposed at a changing unknown boundary in the computation process. As a result, they reported convergence difficulties, which can be avoided by moving towards more complex methods. Moreover, to guarantee reasonable accuracy, they were forced to use a small enough step-size and extensive computations for the solution of the IVPs. Na [39, pp. 280-286] described the application of invariant imbedding to the Falkner-Skan model. For this problem, a modified shooting method [1] and finite-difference methods [2,3] were presented by Asaithambi. Sher and Yakhot [44] defined a new approach to solve this problem by shooting from infinity, using some simple analysis of the asymptotic behaviour of the solution at infinity. Kuo [37] used a differential transformation method, to compute a series solution of the Falkner-Skan equation. Asaithambi [4] proposed a faster shooting method by using a recursive evaluation of Taylor coefficients. Zhang and Chen [50] investigated a modification of the shooting method, where the computation of the Jacobian matrix was obtained by solving two IVPs. A Galerkin-Laguerre spectral method was defined and applied to the Falkner-Skan model by Auteri and Quartapelle [5].

The iterative transformation method
By requiring the invariance of (1.1) with respect to a given transformation group we characterize a subclass of problems (see [28] for the stretching and the spiral group and [13] for the translation group). As a consequence non-ITMs are applicable only to special classes of problems. To overcome this drawback we let the function (⋅, ⋅, ⋅, ⋅) in (1.1) to depend also on a numerical parameter ℎ and we use ℎ also to modify the initial conditions to get invariant initial conditions, and in this way, we introduce the problem This allows us to consider the extended stretching group * = ; where and are constants different from zero and is the exponential of the group parameter, and therefore it is not zero. Our intention is to require the invariance of the governing differential equation in (2.1) with respect to (2.2). To this end, we have to require that the following relation , , , Of course, a given problem belonging to (1.1) can always be embedded into From a numerical viewpoint, we have to consider the initial value problem (IVP) in the starred variables defined by where and have to be considered as fixed values. If for every value of ℎ * the problem (2.5) is well posed, then * (0) is uniquely defined. Consequently, the application of invariance considerations allows us to obtain, if ∞ ≠ 0, .
(2.6) If ∞ = 0, then we can use ℎ to get a non-invariant boundary condition at infinity.
In fact, we can use the new side condition Therefore, if ∞ = 0 then we apply the new boundary condition (2.7) and get by . (2.8) In both cases, for a positive value of we get the value of ℎ given by A solution of (1.1) is determined when the value of ℎ = 1 is obtained from (2.9). If and are fixed, then is a function of ℎ * only, i.e., = (ℎ * ) where, in general, the functional form of (⋅) is not known. Thus, we are interested in the zeros of the "transformation function" This function is defined implicitly by the solution of the IVP (2.5).
Along the lines of the analysis sketched above the ITM can be defined as follows: the original BVP is embedded into (2.1) by fixing non-vanishing values of and ; starting with suitable values of ℎ * 0 and ℎ * 1 a root-finder method, secant, regula-falsi, quasi-Newton or bisection, is used to define a sequence ℎ * , for = 2, 3, … , . At each iteration, a related sequence given by Γ(ℎ * ), for = 0, 1, 2, … , is defined according to (2.10); suitable termination criteria have to be used to verify whether Γ(ℎ * ) → 0 as → ∞; a numerical approximation of ( ) can be obtained by rescaling the numerical solution * ( * ) corresponding to ℎ = 1.

Existence and uniqueness
The main result of this paper can be stated as follows: for a given BVP defined on a semi-infinite interval the existence and uniqueness question is reduced to find the number of real zeros of the transformation function. This result is proved below.
Theorem 1 Let us assume that and are fixed and that for every value of ℎ * the IVP (2.5) is well-posed. Then, the BVP (1.1) has a unique solution if and only if the transformation function has a unique real zero; nonexistence (nonuniqueness) of the solution of (1.1) is equivalent to nonexistence of real zeros (existence of more than one real zero) of Γ(⋅).
Proof by invariance considerations. We begin by proving that there exists a one-to-one and onto correspondence between the set of solutions of (1.1) and the set of real zeros of the transformation function. The thesis is an evident consequence of this result.
The mentioned correspondence can be defined as follows. First, for every values of and different from zero, given a solution ( ) of (2.1) we can associate to it the real zero of Γ(⋅) defined by ℎ * = ∞ / * * (∞). The related value of = ℎ * 1/ allows us to verify by substitution in (2.10) that we have defined a real zero of Γ(⋅). Moreover, because of the invariance with respect to (2.2) the related solution of (2.5) is given by * ( * ) = ℎ * 1/ (ℎ * / ). It is easily seen that both a right and a left inverse of our correspondence are fixed by means of the relations defined above. Therefore, the correspondence is one-to-one and onto. Finally, it follows that if one of the two sets is empty then the other one is empty too.  Moreover, by the proof of the theorem above we know in advance that every real zero of the transformation function belongs to IR + and therefore we can restrict the domain of Γ(⋅) to IR + whereupon its range will be reduced to (−1, ∞). We remark that the values of and can always be chosen to consider only positive values of ℎ * .
On the other hand, the sensitivity of Γ(⋅) (and consequently of Ξ(⋅, ⋅, ⋅)) with respect to ℎ * is of paramount interest. In fact, the numerical determination of the roots of Γ(⋅) = 0 is a well-conditioned problem if and only if | Γ/ ℎ * | at the root is bigger than the machine rounding unit (a simple proof of this statement is possible under suitable regularity conditions on Γ(⋅)). Therefore, it is relevant to verify the sensitivity of Γ(⋅) at each zero. Of course, since Γ(⋅) is not given explicitly we cannot evaluate its sensitivity directly. Hence, monitoring of the sensitivity will require an increase in the computational cost. A simple procedure is to apply as root-finder the secant method because in this case, we compute at each iteration a finite difference approximation for the first derivative of Γ(⋅). Some caution has to be used because finite difference approximations of derivatives are prone to rounding errors that, for really small increments, could dominate the approxi-

Numerical tests and results
In this section, we test the ITM with respect to the existence and uniqueness question and obtain a numerical solution for the given BVP.

Sakiadis problem
In order to apply the ITM to (1.2) we introduce, see Fazio [25], the extended In particular, we are interested to compute * * ( * ∞ ), where * ∞ is a suitable truncated boundary, an approximation of the asymptotic value * * (∞), which is used in the definition (2.8) of .
It is evident that our numerical method is based on the behaviour of the transformation function. Our interest is to study the behaviour of this function with respect to its independent variable. We notice that because of the two terms ℎ * 1/2 , which have been introduced in the modified boundary conditions in (5.1), we are allowed to consider only positive values of ℎ * .
From our numerical study concerning the dependence of Γ with respect to the missing initial condition For a negative missing initial condition the numerical results are shown on figure 2. It is evident from figure 2 that the transformation function has only one zero and, by the theorem proved above, this means that the considered problem has one and only one solution. Moreover, we remark that the tangent to the Γ function at its unique zero and the ℎ * axis define a large angle. From a numerical viewpoint, this means that the quest for the ℎ * corresponding to ℎ = 1 is a wellconditioned problem.
As far as the numerical results reported in this section are concerned, the ITM was applied by setting the truncated boundary * ∞ = 10. Moreover, these results were obtained by an adaptive fourth-order Runge-Kutta IVP solver. The adaptive solver uses a relative and absolute error tolerance, for each component This solution was computed by rescaling, note that by rescaling we get * ∞ < ∞ , and this means that we have reduced the computational cost, where the chosen truncated boundary was * ∞ = 10 in our case.  where is a parameter. In the following we set = 4; for the choice = 8 see [18]. In [18], a free boundary formulation of the Falkner-Skan model was considered and numerical results were computed for the Homann flow ( = 1/2) as well as for the Hiemenz flow ( = 1).
From a numerical point of view the request to evaluate (∞) cannot be fulfilled. Several strategies have been proposed in order to provide an approximation of this value. The simplest and widely used one is to introduce, instead of infinity, a suitable truncated boundary. For the sake of simplicity we apply, following Töpfer [47], some preliminary computational tests to find a suitable value for the truncated boundary.
As far as the reverse flow solutions are concerned, in   In figure 6 we plot the behaviour of missed initial condition versus . The solution found by the data in table 3 is plotted in this figure, but not the one found in table 2 because this is very close to the Blasius solution. A good initial choice of the initial iterates of ℎ * , for a given value of , is obtained by employing values close to the one used in a successful attempt made for a close value of .
It is interesting to note that, for values of < min the ITM continued to iterate endlessly, whatever a couple of starting values for ℎ * are selected.
Our extended algorithm has shown a kind of robustness because it is able to get convergence even when, for a chosen value of ℎ * , the IVP solver stops before arriving at the selected truncated boundary getting a wrong value of Γ(ℎ * ) = −1.
On the other hand, the secant method gives an overflow error when this happens for two successive iterates of ℎ * .
The value of min , corresponding to a separation point at = 0, can be found by the ITM by considering as a continuation parameter. In figure 7 we plot the unique solution for the limiting value min , where min is given by equation (1.4).
Let us notice that this is a normal flow solution. The results reported so far have been found by a variable order adaptive multi-step IVP solver that was coupled with the simple secant method. The adap-tive solver uses a relative and absolute error tolerance, for each component of the numerical solution, both equal to 10D − 6. As well known, the secant method is convergent provided that we select two initial iterates sufficiently close to the root, and its convergence is super-linear with an order of convergence equal to (1 + √ 5)/2. As far as a termination criterion for the secant method is concerned, we enforced the conditions |Γ(ℎ * )| ≤ Tol and |ℎ * − ℎ * −1 | ≤ TolR|ℎ * | + TolA , (5.10) with Tol = TolR = TolA = 1D − 06.
For other problems, in boundary layer theory, admitting more that one solution or none, depending on the value of a parameter involved, see [23] or [24].

Final remarks and conclusions
In this paper, we define an ITM to provide numerical evidence for the existence and uniqueness of a solution of BVPs defined on infinite intervals on the basis of the theorem proved in section 3. The leading point is to establish whether the introduced transformation function has only one real zero or not. From a numerical viewpoint, several computations may be necessary to understand and to characterize the behaviour of Γ(⋅). For the proposed numerical test the transformation function was calculated at test-points while at some points a root-finding method was used. In the process we tried to bracket the roots of Γ(⋅) = 0. In this way, it was possible to obtain numerical results that are in very good agreement with the values available in the literature. By setting different values of and we get a different transformation function and consequently we must bracket its zeros once again.
Of course, the asymptotic boundary condition is not easily used numerically and for the sake of simplicity, we have chosen to replace it with the same condition at a suitable truncated boundary. A recent successful way to deal with such an issue is to reformulate the considered problem as a free BVP [15,18,20]; for a survey on this topic see [22]. Recently, Zhang and Chen [50] have used a free boundary formulation to compute the normal flow solutions of the Falkner-Skan model in the full range min < ≤ 40. They applied a modified Newton's method to compute both the initial velocity and the free boundary. A recent approach, proposed by Fazio and Jannelli [29,30] see also [31,32], used to enforce the asymptotic boundary condition exactly, is to apply a quasi-uniform grid, that has the last mesh point at infinity, and suitable non-standard finite difference