A Numerical Algorithm for Solving Higher-Order Nonlinear BVPs with an Application on Fluid Flow over a Shrinking Permeable Infinite Long Cylinder

Wepresent an efficient iterative power seriesmethod for nonlinear boundary-value problems that treats the typical divergence problem and increases arbitrarily the radius of convergence. This method is based on expanding the solution around an iterative initial point. We employ this method to study the unsteady, viscous, and incompressible laminar flow and heat transfer over a shrinking permeable cylinder. More precisely, we solve the unsteady nonlinear Navier–Stokes and energy equations after reducing them to a system of nonlinear boundary-value problems of ordinary differential equations. The present method successfully captures dual solutions for both the flow and heat transfer fields and a unique solution at a specific critical unsteadiness parameter. Comparisons with previous numerical methods and an exact solution verify the validity, accuracy, and efficiency of the present method.


Introduction
Numerous phenomena in engineering and applied science fields are governed by nonlinear boundary-value problems (BVPs).Therefore, BVPs have received a huge attention from mathematicians, physicists, and engineers for the sake of finding and analyzing their solutions.Generally speaking, finding the analytical solutions for nonlinear BVPs is far from trivial and often is impossible.Therefore, many numerical techniques have been developed to solve such type of problems.These methods include Adomian's decomposition method, homotopy perturbation method, variational iteration method, optimal homotopy asymptotic method, operational matrices techniques based on various orthogonal polynomials and wavelets, finite difference method, and spectral methods; the reader is referred to [1][2][3][4][5][6] and references therein.
The fluid dynamics and heat transfer of a viscous incompressible fluid flowing past stretching surfaces, such as a sheet or tube, have attracted considerable interests of many researchers because of their importance in many industrial applications such as the quality of certain products.One of the most interesting conditions for stretching surfaces problems is the velocity at the surface, where it mainly figures the characteristics of the fluid based upon two essential factors, fluid viscosity and suction parameter.A remarkable interest of several researchers concentrated on tracking the existence of dual solutions for the flow within a certain range of unsteadiness and suction parameters [7][8][9][10][11][12][13][14][15][16][17].Although, the literature reveals numerous research papers discussing the flow over a stretching sheet and moving plate [18][19][20][21][22][23][24][25][26], there are only few studies focusing on the problem of flowing past a stretching cylinder or tube; see [8][9][10] and references therein.
It is established that the differential equations describing the fluid flow BVPs are highly nonlinear and demand extremely accurate numerical schemes.In this work, we show that an iterative procedure, based on successive power series expansions, provides one such high accuracy numerical scheme.Often, using power series solutions turns to be useless because the resulting solution diverges at a finite radius of convergence.The divergence is intrinsic to the nature of the solution in the sense that it persists to exist even with an infinite power series expansion.The method presented here 2 Complexity solves this problem showing that the radius of convergence can be delayed arbitrarily to any large value.This value could, in principle, approach infinity achieving exact solutions.
Recognizing that a powerful numerical scheme based on this method is already established [34][35][36][44][45][46][47][48], we nonetheless present a thorough investigation of the error associated with this method with the aim of showing how we can systemically reduce errors to infinitesimal values while having the Central Processing Unit (CPU) time within a reasonable range.We will show robustness and efficiency of the method via a highly demanding fluid flow boundaryvalue problem.Therefore, solving the problem of finite radius of convergence will open the door wide for applying the power series method to much larger class of differential equations, particularly the nonlinear ones.
Briefly, the present technique is based on iterative power series expansions of the solution.The domain of the independent variable, say , is divided into a number  of segments each of width Δ, where Δ is smaller than the radius of convergence.A power series solution is obtained by expanding the solution around the left end of the first segment using the initial conditions given with the problem.Similarly, a power series solution is obtained by expanding around the start of the second segment but now using the first series to calculate the initial conditions.This is repeated  times till a solution at  =  × Δ is obtained.In the limit  → ∞ and Δ → 0 the series solution becomes an exact solution.This scheme is effectively equivalent to an iterative procedure of repeated iterative calculation of the recursion relations of the power series in the first segment.Another aim of this paper is to apply this method to study the unsteady flow and heat transfer characteristics of fluid flow over a shrinking permeable infinite long cylinder.We will show that the iterative numerical scheme resulting from this method is exceeding the efficiency of typical numerical methods used.In addition, we managed to find an exact solution which enabled us to calculate accurately the error for a finite value of number of iterations .It should be noted that the present work is a part of the Master thesis [49].
The rest of the paper is organized as follows.In Section 2, a mathematical representation of our method is illustrated using a general form of nonlinear ordinary differential equation, while henceforth we call it iterative power series method.Sections 3 and 4 display the implementation of the iterative power series method on the heat and mass transfer model.Section 5 focuses on analyzing the validity of this present technique and demonstrating its efficiency by drawing comparisons with the achieved exact analytical solution and other numerical methods.In Section 6, we analyze and discuss the properties of the solutions obtained.We end with a summary of our main conclusions in Section 7.

General Scheme of the Iterative Power Series (IPS) Method
In this section, we give a brief description of the present method.Consider a general ordinary differential equation of the form with  initial conditions where  () is the th derivative of (),   are real constants, and () is a known function.The factor ! is introduced, without loss of generality, for the constants   to correspond to the coefficients of the power series expansion below.At first, we divide the interval [ 0 , ] into a number of  identical segments each of width Δ = ( −  0 )/.Then we expand () in a power series around the beginning of each interval, namely, where   () is the power series expansion around the start of the th segment and    are the coefficients of the power series.Recursion relations between the coefficients are obtained upon substituting the power series solution, (3), into the differential equation, (1), which can be expressed in terms of the first  coefficients corresponding to the initial conditions where {   } denotes the set of coefficients   0 ,   1 , . . .,   −1 .The essential idea of the IPS method is to calculate the coefficients { +1  } of the (+1)th power series from the th series according to (2), which upon using (3) reads and is simplified to and then imposes the condition  max > .Here, ( +  ) is the binomial function.The last equation is the basis for the IPS algorithm.Starting from the initial conditions { 0  } for the power series of the zeroth interval, an iterative application of (7) leads to the coefficients of the th interval, namely, {   } which give the solution at the desired point,  =  0 + Δ, Both analytical and numerical schemes may be deduced from this algorithm.For the numerical scheme, the value of Δ used is inserted as a number Δ = ( ∞ −  0 )/.On the other hand, leaving  as a variable results in an analytical solution in terms of a power series in  which is equivalent to a functional transformation on the zeroth order series; that is, the coefficients of the th series are functional transformation of the ( − 1)th series.In such a case the last power series for the th interval corresponds to  such that functional transformations and all power series expansions of the zeroth up to ( − 1)th intervals will be included in the th expansion.
The coefficient   0 of each th expansion represents the value of the solution at  =  0 + Δ, which gives a discrete representation of ().Therefore, in the limit  → ∞, the discrete representation turns to a continuous one and thus we conjecture that the exact solution is obtained in the limits of

Heat and Mass Transfer Model
In this section we will employ the present numerical technique on heat and mass transfer over a shrinking permeable cylinder described in Zaimi et al. [8] and Elnajjar et al. [10].For completeness, we redescribe precisely the physical model.The flow is considered an unsteady, laminar, viscous, and incompressible fluid with uniform velocity  and uniform temperature  ∞ over a permeable shrinking circular cylinder.The cylinder is assumed to be infinitely long and the flow has constant properties.The diameter of the cylinder is assumed to be time dependent with the radius () =  0 √1 − , where  0 is a positive constant,  is the constant of expansion/contraction strength, and  is the time.
Clearly, the cylinder's radius is shrinking with time if  is positive and stretching with time if  is negative.Notice that since the flow is axisymmetric, the flow field should be a function of the radial coordinate, , and the longitudinal.The governing equations for the unsteady and incompressible fluid without body force are the continuity, momentum, and energy equations.These equations in cylindrical coordinate system, (, ), are given by where  and  are the polar coordinates in the radial and axial directions, respectively,   and   are the fluid velocity components in the radial and axial directions, respectively, and  is the fluid temperature.The function  represents the fluid pressure and the parameters ], , and  denote the fluid viscosity, the fluid density, and the fluid thermal diffusivity, respectively.Notice that we assumed that there is no azimuthal velocity component.The assumed boundary conditions associated with (10) for the velocity components and the temperature are given by where   is the constant surface temperature and  0 is a positive constant.The similarity transformations which convert (10) into nonlinear ordinary differential equations are given by [8,10] where   () = / and  is the similarity variable given by In addition, it should be noted that  represents the dimensionless stream function and  represents the normalized temperature.Applying the above similarity transformations, (10) and the boundary conditions (11) reduce to subject to where  =  2 0 /4] is the unsteadiness parameter representing the strength of contraction/expansion,  = − 0 /2] is the suction parameter, and Pr = ]/ is the Prandtl number.
Our main target is to solve ( 14) and ( 15), subject to the boundary conditions (16), using the present technique in the ranges 0 ≤  ≤ 7 and −4 ≤  ≤ 0 at Pr = 0.7.In this study, we will analyze the normalized skin friction coefficient,   (1), and the normalized heat transfer rate, −  (1).

IPS Method for Heat and Mass Transfer Model
The following is a detailed implementation of the IPS method used to solve ( 14) and ( 15) subject to the boundary conditions (16).It is worth mentioning that ( 14) includes only the variable , while (15) includes both  and .Therefore, it is more convenient to find the variable  from ( 14) and then solve (15).Furthermore, for the sake of simplicity, we render ( 14) to an initial-value problem; that is, with where  must be chosen using the shooting method [10] so that the solution satisfies the boundary condition   (∞) = 0. Notice that by setting different initial values for  in the shooting method, dual solutions are obtained.
As mentioned in Section 2, we start by expanding () in power series around the initial point and the derivatives,  0  ,  0  , and  0  can be calculated simply by differentiating this series.Substituted in (14), the coefficients,  0  , for  ≥ 3, can be found recursively in terms of the initial conditions  0 0 ,  0 1 , and  0 2 through the recursion relations.The first two recursion relations are given by Recalculating  0 (),  0  (), and  0  () at  = 1 + Δ gives Now,  1 0 ,  1 1 , and  1 2 play the role of the initial conditions for the next series expansion, where we expand the solution and its derivatives in power series around Resubstituting these power series expansions in the differential equation, we get the new recursion relations 2 ).The next iterative step is to calculate  1 () and its first two derivatives at  = 1 + 2Δ which will give the initial conditions for the new power series.Repeating this process  times, the general forms of the first two recursion relations of ( 14) are found to be where In summary, the IPS procedure can be reduced to the following algorithm: where  3 =  3 ( 0 ,  1 ,  2 ) and  4 =  4 ( 0 ,  1 ,  2 ) are the recursion relations obtained from the differential equation.
We have removed the superscripts that indicate the index of the iteration for convenience.The scheme is thus described simply as follows: one starts with (26) to calculate (Δ), followed by updating the initial conditions according to (27) and then using the updated values back in ( 26) and so on.The procedure has to be repeated  times with Δ = ( − 1)/.Similarly, for the energy equation subject to where  must be chosen using the shooting method [10] so that the solution satisfies the boundary condition (∞) = 0. We expand () in power series around the same initial point and similarly for the derivatives,  0  and  0  .Substituted in (15), the coefficients,  0  , for  ≥ 2, are found recursively in terms of the initial conditions  0 0 and  0 1 through the recursion relations.Employing the IPS method, the first two recursion relations are found to be where

Validation
In this section we aim at demonstrating the performance and efficiency of the present numerical scheme.Firstly, in order to obtain accurate numerical results, we have to pay attention to the selection of the numerical algorithm parameters, ,  max , and  ∞ .To achieve this target we focus on studying the following: with where  is updated using the shooting method to satisfy the condition   (∞) = 0. Notice that several recent studies such as [10][11][12][13][14][15][16] reported the existence of a critical value of  (named   ) at which the problem has no solution (for  >   ), only one solution (at  =   ), and dual solutions (for  <   ).It is worth mentioning that we succeeded in finding the explicit analytical form of the first solution for (33) subject to (34) under a condition  = −1/; that is This exact solution will play a crucial role in proving the advantages of our numerical scheme.The approximate solutions of the problems ( 33) and (34) at three different iterations  = 1, 2, and 3, together with the exact solution obtained by (35) when  = −1 and  = 1, are displayed in Figure 1.It is clearly seen that increasing the number of iterations in the IPS method delays the divergence point.
To achieve an "optimal choice" of  ∞ , we solve the problem with  ∞ = 7, 8, . . ., 17. Table 1 shows the values of  up to 50 digits corresponding to the values of  ∞ .It is clearly  seen that the value of  stabilizes at around  ∞ = 15; hence we choose  ∞ = 15 as the optimal value for the rest of the calculations in the entire paper.It should be noted that most of the used numerical schemes for such type of problems, [8][9][10]17], had chosen  ∞ = 7 to represent the infinity which, definitely, gave lower order of accuracy.This conclusion can be easily tested via the exact solution (35) which gives   (7) = −0.00247875and   (15) = −8.31529×10−7 .However, we will only show the interval up to  ∞ = 8 for the rest of the coming figures.
The upper bound of the error in the IPS method can be estimated as follows.At each iterative step an error of order Δ  max +1 results from terminating the power series at  max .This error will be magnified  times due to the iterative procedure.As a result, the upper bound of the error of the IPS method is Table 2 presents the CPU time in seconds, which is machine-dependent, versus the upper bound of the error,  IPS , of the IPS method for the first solution at  = 2 and  = −1, where  max varies from 3 to 8.
The exact solution, (35), provides a unique possibility of calculating the error of the IPS method and comparing it with that of other numerical methods.Figure 2 presents a comparison between the error of the IPS method and the explicit Runge-Kutta method of order four (ERK4) for the problem at  = 1 and  = −1.The advantages of the present technique over the other one are notable.
A further comparison is done in Figure 3 on the normalized skin friction coefficient,   (1), as a function of  with Zaimi et al. [8] and Elnajjar et al. [10] for the case when  = 0. Excellent agreements are obtained.It should be mentioned herein that Elnajjar et al. [10] used a combination of the implicit Runge-Kutta method and the shooting method while Zaimi et al. [8] implemented the shooting method described in the book by Jaluria and Torrance [11].
It is worth mentioning here that our technique successfully showed its definite capability to exceed the machine precision which is 10 −14 .A clear sign on this is the systematic reduction in the error in Table 2 when compared with the exact solution.Moreover, all computations are performed with at least 14 decimal digits of precision, knowing that all computations are operated using Mathematica software 10.4 and carried out on a Lenovo PC with the following The error is defined as the difference between the numerical solution and the exact solution (35).specifications: model: Z470, processor: Core(TM) i5-2430M CPU @ 2.40 GHz, system type: 64-bit, and installed memory (RAM): 4.00 GB.

Results and Discussion
In this section, we discuss the effects of both suction parameter, , and unsteadiness parameter, , on the velocity profile,   (), the normalized skin friction coefficient,   (1), the temperature profile, (), and the heat transfer rate, −  ().
The numerical simulations are conducted at a fixed Prandtl number, Pr = 0.7, while the ranges considered for the other parameters are 0 ≤  ≤ 7 and −4 ≤  ≤ −1.
Figure 4 shows the first and second solutions of the velocity profiles for  = 0, 1, 3, 5, 7 with a fixed unsteadiness parameter,  = −2.It is clearly seen that the first solution for the fluid velocity inside the boundary layer region increases as  increases, while the second solution shows an opposite trend.In addition, the two solutions of the velocity profile  become steeper with higher magnitudes as  increases.These observations emphasize the effect of increasing the suction parameter of the cylinder's wall which is to decrease the boundary layer thickness.Consequently, increasing the suction parameter causes an increment in the normalized skin friction coefficient for the first solution and decrement  in the normalized skin friction coefficient for the second solution, as clearly shown in Figure 5.These findings are consistent with the results reported by Elnajjar et al. [10] and Zaimi et al. [8].
Figure 6 presents the temperature profiles of the fluid flow, (), at  = −2 and  = 0, 1, 3, 5, 7. It is obviously noticeable that both solutions for temperature profiles admit similar behaviour, where they become wider and more relaxed as the suction parameter decreases.These kinds of behaviour inspire us to conclude that the developed thermal boundary layer and the corresponding rate of heat transfer are decreasing as  increases.However, the second solution depicts more relaxed behaviour compared with the first solution.This slight difference between the first and second temperature profiles indicates that the second solution reflects higher thermal boundary layer than the first solution and, thus, a larger rate of the heat transfer as confirmed by Figure 7.The variation of both the normalized skin friction coefficient,   (1), and the heat transfer rate, −  (1), as functions of , are shown, respectively, in Figures 5 and 7 for  = 0, 0.5, 1, 1.5, 2. The results demonstrate the existence of a critical value   in the -domain at which the problem has no solution for  >   , only one solution at  =   , and dual solutions for  <   .Figure 5 shows that |  (1)| increases as  increases which is due to the increase in the surface shear stress coefficient.Moreover, we observe that |  (1)| is decreasing with .However, Figure 7 clearly shows that increasing  will definitely increase the heat transfer rate while increasing  causes a decrease in the heat transfer rate.
Figure 8 displays the first and second solutions of the velocity profiles for  = −1, −2, −3, −4 with a fixed value of the suction parameter,  = 1.Generally speaking, the behaviour of   () is very similar to the case of the variable suction parameter; that is, increasing the unsteadiness parameter produces steeper behaviour in the velocity profiles for the first solution while the second solution shows an opposite trend.In agreement with the case of the variable suction parameter in Figure 4, increasing the unsteadiness parameter will then cause a reduction in the thickness of the boundary layer.
Figure 9 presents the temperature profiles of the fluid flow, (), at  = 1 and  = −1, −2, −3, −4.Clearly, the increase in the unsteadiness parameter or the suction parameter leads to the same trend.
Finally, we end up our discussion with Figure 10 which presents an overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the - domain  for Pr = 0.7 and 0 ≤  ≤ 2. The straight line in this figure represents the occurrence of unique solution of the problem.

Conclusions
In this work, we presented a numerical technique for solving nonlinear BVPs based on iterative power series solutions.
We have demonstrated its efficiency and accuracy through validation against the numerical ERK4.We have shown that our method excels over the ERK4 by orders of magnitude in accuracy.Moreover, the accuracy in our method is systematically controlled such that the error can be reduced to any arbitrary small value.We successfully studied the unsteady viscous flow over a contracting cylinder using the present technique.The velocity and temperature profiles of the ordinary version of Navier-Stokes equations for different suction and unsteadiness parameters are calculated.We have obtained the exact solution of the first solution of the fluid flow under a specific condition,  = −1/, and have employed it to emphasize the efficiency of the present numerical technique.
The convergence analysis of the IPS method is considered for a future work.However, we strongly believe that the method is convergent.This is conjectured by the systematic Complexity reduction in error upon increasing the number of iterations or the number of terms in the power series.We believe that this technique will serve researchers in different fields working on nonlinear systems.In particular, the technique will be very useful for systems described by nonintegrable nonlinear differential equations.

Figure 2 :
Figure 2: Error of the IPS method (a) and ERK4 (b) for the case of  = 1 and  = −1.The error is defined as the difference between the numerical solution and the exact solution(35).

Figure 5 :
Figure 5: Normalized skin friction coefficient,   (1), as a function of  for different values of : the first solution (solid) and second solution (dotted).

2 Figure 7 :
Figure 7: Heat transfer rate, −  (1), as a function of  for different values of : the first solution (solid) and second solution (dotted).

Table 2 :
The upper bound of the error for the first solution at  = 2 and  = −1 versus the CPU time at different values of  max .