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

Authors: Laila Y. Al Sakkaf; Qasem M. Al-Mdallal; U. Al Khawaja
Journal: Complexity
 (2018)

---

## Abstract

We present an efficient iterative power series method 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.


---

## Body

## 1. 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–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–17]. Although, the literature reveals numerous research papers discussing the flow over a stretching sheet and moving plate [18–26], there are only few studies focusing on the problem of flowing past a stretching cylinder or tube; see [8–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 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.Among the many different methods of solving nonlinear differential equations [27–33], the power series method is the most straight forward and efficient [34]. It has been used as a powerful numerical scheme for many problems [35–43] including chaotic systems [44–47]. Many numerical algorithms and codes have been developed based on this method [34–36, 44–48]. However, the abovementioned finiteness of radius of convergence is a serious problem that hinders the use of this method to wide class of differential equations, in particular the nonlinear ones. For instance, the nonlinear Schrödinger equation (NLSE) with cubic nonlinearity has sech⁡(x) as a solution. Using the power series method to solve this equation produces the power series of sech⁡(x), which is valid only for x<π/2.Recognizing that a powerful numerical scheme based on this method is already established [34–36, 44–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 boundary-value 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 I 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 I times till a solution at η=I×Δ is obtained. In the limit I→∞ 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 I. 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 Section2, a mathematical representation of our method is illustrated using a general form of nonlinear ordinary differential equation, while henceforth we call ititerative 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.


## 2. 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(1)Ffη,f′η,f′′η,…,fmη,gη=0,η∈η0,η∞,with m initial conditions(2)fiη0=ai×i!,i=0,1,2,…,m-1,where f(i) is the ith derivative of f(η), ai are real constants, and g(η) is a known function. The factor i! is introduced, without loss of generality, for the constants ai to correspond to the coefficients of the power series expansion below. At first, we divide the interval [η0,η] into a number of I identical segments each of width Δ=(η-η0)/I. Then we expand f(η) in a power series around the beginning of each interval, namely,(3)fiη=∑n=0nmaxaniη-η0+iΔn,iΔ≤η≤i+1Δ,0≤i<I,where fi(η) is the power series expansion around the start of the ith segment and ani 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 m coefficients corresponding to the initial conditions(4)ani=aniaki,0≤k<m,n≥m,where {aki} denotes the set of coefficients a0i,a1i,…,am-1i. The essential idea of the IPS method is to calculate the coefficients {aki+1} of the (i+1)th power series from the ith series according to (2),(5)aki+1=1k!dkdηkfiηη=η0+i+1Δ,which upon using (3) reads(6)aki+1=1k!∑n=0nmaxan+ki×dkdΔkΔn+kand is simplified to(7)aki+1=∑n=0nmaxan+ki×n+kkΔnand then imposes the condition nmax>k. Here, n+kk is the binomial function. The last equation is the basis for the IPS algorithm. Starting from the initial conditions ak0 for the power series of the zeroth interval, an iterative application of (7) leads to the coefficients of the Ith interval, namely, {akI} which give the solution at the desired point, η=η0+IΔ,(8)fIη=∑n=0nmaxanI×η-η0+IΔn.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)/I. 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 ith series are functional transformation of the (i-1)th series. In such a case the last power series for the Ith interval corresponds to I such that functional transformations and all power series expansions of the zeroth up to (I-1)th intervals will be included in the Ith expansion.The coefficienta0i of each ith expansion represents the value of the solution at η=η0+iΔ, which gives a discrete representation of f(η). Therefore, in the limit I→∞, the discrete representation turns to a continuous one and thus we conjecture that the exact solution is obtained in the limits of I→∞ and nmax→∞(9)fη=limI→∞⁡∑n=0∞anI×η-η0+IΔn.


## 3. 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 U and uniform temperature T∞ 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 a(t)=r01-βt, where r0 is a positive constant, β is the constant of expansion/contraction strength, and t 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, r, 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, (r,z), are given by(10)1r∂∂rrur+∂uz∂z=0,∂ur∂t+ur∂ur∂r+uz∂ur∂z=-1ρ∂p∂r+ν∂2ur∂r2+1r∂ur∂r+∂2ur∂z2-urr2,∂uz∂t+ur∂uz∂r+uz∂uz∂z=-1ρ∂p∂z+ν∂2uz∂r2+1r∂uz∂r+∂2uz∂z2,∂T∂t+ur∂ur∂r=α1r∂∂rr∂T∂r,where r and z are the polar coordinates in the radial and axial directions, respectively, ur and uz are the fluid velocity components in the radial and axial directions, respectively, and T is the fluid temperature. The function p 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(11)ur=-2νγr01-βt,uz=-4νzr021-βt,Ts=c01-βt+T∞atr=at,uz=0,T=T∞asr⟶∞,where Ts is the constant surface temperature and c0 is a positive constant.The similarity transformations which convert (10) into nonlinear ordinary differential equations are given by [8, 10](12)ur=-2νr01-βtfηη,uz=4νzr021-βtf′η,T=Ts-T∞θη-T∞+T∞,where f′(η)=df/dη and η is the similarity variable given by(13)η=rr0211-βt.In addition, it should be noted that f represents the dimensionless stream function and θ represents the normalized temperature. Applying the above similarity transformations, (10) and the boundary conditions (11) reduce to(14)ηf′′′+f′′+ff′′-f′2-Sηf′′+f′=0,(15)ηθ′′+θ′1+Prf-SPrη-SPrθ=0,subject to(16)f1=γ,f′1=-1,f′∞=0,θ1=1,θ∞=0,where S=r02β/4ν is the unsteadiness parameter representing the strength of contraction/expansion, γ=-r0U/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≤S≤0 at Pr=0.7. In this study, we will analyze the normalized skin friction coefficient, f′′(1), and the normalized heat transfer rate, -θ′(1).


## 4. 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 f, while (15) includes both θ and f. Therefore, it is more convenient to find the variable f from (14) and then solve (15). Furthermore, for the sake of simplicity, we render (14) to an initial-value problem; that is,(17)f′′′=1η-f′′-ff′′+f′2+Sηf′′+f′,with(18)f1=γ,f′1=-1,f′′1=λ,where λ must be chosen using the shooting method [10] so that the solution satisfies the boundary condition f′(∞)=0. Notice that by setting different initial values for λ in the shooting method, dual solutions are obtained.As mentioned in Section2, we start by expanding f(η) in power series around the initial point η0=1(19)f0η=∑n=0nmaxan0η-1n,and the derivatives, f0′, f0′′, and f0′′′ can be calculated simply by differentiating this series. Substituted in (14), the coefficients, an0, for n≥ 3, can be found recursively in terms of the initial conditions a00, a10, and a20 through the recursion relations. The first two recursion relations are given by(20)a30=16a102-2a201+a00-S+a10S,a40=124a102S-2-a00+2a202+a002+a003-2S-S+S2+a102a20+SS-2-a00.Recalculating f0(η), f0′(η), and f0′′(η) at η=1+Δ gives(21)a01=f0Δ,a11=f0′Δ,a21=f0′′Δ2.Now, a01, a11, and a21 play the role of the initial conditions for the next series expansion, where we expand the solution and its derivatives in power series around η0=1+Δ(22)f1η=∑n=0nmaxan1η-1+Δn.Resubstituting these power series expansions in the differential equation, we get the new recursion relations an1(a01,a11,a21). The next iterative step is to calculate f1(η) and its first two derivatives at η=1+2Δ which will give the initial conditions for the new power series. Repeating this process I times, the general forms of the first two recursion relations of (14) are found to be(23)a3I=161+IΔa1I2-2a2I1+a0I-S-SIΔ+a1IS,(24)a4I=1241+IΔ2a1I2-2-a0I+S+SIΔ+a1I2a2I1+IΔ+S-2-a0I+S+SIΔ+2a2IS-1-IΔ+S+SIΔ2+2+a0I2+a0I3-2S-2SIΔ,where(25)a0I=fI-1Δ,a1I=fI-1′Δ,a2I=fI-1′′Δ2.In summary, the IPS procedure can be reduced to the following algorithm:(26)fΔ=a0+a1Δ+a2Δ2+a3Δ3+a4Δ4+OΔ5,(27)a0=fΔ,a1=f′Δ,a2=f′′Δ2,where a3=a3(a0,a1,a2) and a4=a4(a0,a1,a2) 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 f(Δ), 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 I times with Δ=(η-1)/I.Similarly, for the energy equation(28)θ′′=1η-θ′1+Prf-SPrη+SPrθ,subject to(29)θ1=1,θ′1=σ,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(30)θ0η=∑n=0nmaxbn0η-1nand similarly for the derivatives, θ0′ and θ0′′. Substituted in (15), the coefficients, bn0, for n≥ 2, are found recursively in terms of the initial conditions b00 and b10 through the recursion relations. Employing the IPS method, the first two recursion relations are found to be (31)b2I=121+IΔb0IPrS+b1I-1-a0IPr+PrS+IΔPrS,b3I=161+IΔ2b0IPrS-2-a0IPr+S1+IΔPr+b1I2+a0I2Pr2-PrS+Pr2S2-PrSIΔ+2Pr2S2IΔ+Pr2S2I2Δ2-a1IPr1+IΔ+a0IPr3-2PrS1+IΔ,where(32)b0I=θI-1Δ,b1I=θI-1′Δ.


## 5. 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,I, nmax, and η∞. To achieve this target we focus on studying the following:(33)f′′′=1η-f′′-ff′′+f′2+Sηf′′+f′,η∈1,η∞,with(34)f1=γ,f′1=-1,f′′1=λ,where λ is updated using the shooting method to satisfy the condition f′(∞)=0. Notice that several recent studies such as [10–16] reported the existence of a critical value of S (named Sc) at which the problem has no solution (for S>Sc), only one solution (at S=Sc), and dual solutions (for S<Sc). 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 S=-1/γ; that is(35)fη=γe1-η/γ.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 I=1,2, and 3, together with the exact solution obtained by (35) when S=-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.Figure 1
Velocity profiles,f′(η), for the first solution; exact (solid) and approximate (dashed) solutions at different I with S=-1 and γ=1.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–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 f′(7)=-0.00247875 and f′(15)=-8.31529×10-7. However, we will only show the interval up to η∞=8 for the rest of the coming figures.Table 1
Progressing ofλ values with η∞ at γ=1 and S=-1.
η ∞ λ 7 1.06032130948772355259169445286450439527129433161564 8 1.02826614125492573507123370718399515004546832373446 9 1.01244316459533191876326159442255228857229463473780 10 1.00535010609275538041830857835707790239466505001516 11 1.00221279175507729615034859048196423158532888754168 12 1.00084874204304334646862685662321915732040012124886 13 1.00030312215822976659593816307972112761442861473174 14 1.00016671718702637162776598969384662018793573810246 15 1.00003031221582297665959381630797211276144286147317 16 1.00003031221582297665959381630797211276144286147317 17 1.00003031221582297665959381630797211276144286147317The upper bound of the error in the IPS method can be estimated as follows. At each iterative step an error of orderΔnmax+1 results from terminating the power series at nmax. This error will be magnified I times due to the iterative procedure. As a result, the upper bound of the error of the IPS method is(36)EIPS=IΔnmax+1.Table2 presents the CPU time in seconds, which is machine-dependent, versus the upper bound of the error, EIPS, of the IPS method for the first solution at γ=2 and S=-1, where nmax varies from 3 to 8.Table 2
The upper bound of the error for the first solution atγ=2 and S=-1 versus the CPU time at different values of nmax.
n m a x E I P S CPU time (seconds) 3 ≈ 1 0 - 5 0.0156001 4 ≈ 1 0 - 7 0.0312002 5 ≈ 1 0 - 9 0.0780005 6 ≈ 1 0 - 10 0.156001 7 ≈ 1 0 - 12 0.296402 8 ≈ 1 0 - 14 0.546004The 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 S=-1. The advantages of the present technique over the other one are notable.Figure 2
Error of the IPS method (a) and ERK4 (b) for the case ofγ=1 and S=-1. The error is defined as the difference between the numerical solution and the exact solution (35).
(a) (b)A further comparison is done in Figure3 on the normalized skin friction coefficient, f′′(1), as a function of S 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].Figure 3
Normalized skin friction coefficient,f′′(1), as a function of S for γ=0.It is worth mentioning here that our technique successfully showed its definite capability to exceed the machine precision which is10-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 specifications: model: Z470, processor: Core(TM) i5-2430M CPU @ 2.40 GHz, system type: 64-bit, and installed memory (RAM): 4.00 GB.


## 6. Results and Discussion

In this section, we discuss the effects of both suction parameter,γ, and unsteadiness parameter, S, on the velocity profile, f′(η), the normalized skin friction coefficient, f′′(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≤S≤-1.Figure4 shows the first and second solutions of the velocity profiles for γ=0,1,3,5,7 with a fixed unsteadiness parameter, S=-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 4
Velocity profiles,f′(η), for different values of γ at S=-2: first solution (a) and second solution (b).
(a) (b)Figure 5
Normalized skin friction coefficient,f′′(1), as a function of S for different values of γ: the first solution (solid) and second solution (dotted).Figure6 presents the temperature profiles of the fluid flow, θ(η), at S=-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 S 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.Figure 6
Temperature profiles,θ(η), for different values of γ at S=-2: first solution (a) and second solution (b).
(a) (b)Figure 7
Heat transfer rate,-θ′(1), as a function of S for different values of γ: the first solution (solid) and second solution (dotted).The variation of both the normalized skin friction coefficient,f′′(1), and the heat transfer rate, -θ′(1), as functions of S, 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 Sc in the S-domain at which the problem has no solution for S>Sc, only one solution at S=Sc, and dual solutions for S<Sc. Figure 5 shows that f′′(1) increases as γ increases which is due to the increase in the surface shear stress coefficient. Moreover, we observe that f′′(1) is decreasing with S. However, Figure 7 clearly shows that increasing γ will definitely increase the heat transfer rate while increasing S causes a decrease in the heat transfer rate.Figure8 displays the first and second solutions of the velocity profiles for S=-1,-2,-3,-4 with a fixed value of the suction parameter, γ=1. Generally speaking, the behaviour of f′(η) 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 8
Velocity profiles,f′(η), for different values of S at γ=1: first solution (a) and second solution (b).
(a) (b)Figure9 presents the temperature profiles of the fluid flow, θ(η), at γ=1 and S=-1,-2,-3,-4. Clearly, the increase in the unsteadiness parameter or the suction parameter leads to the same trend.Figure 9
Temperature profiles,θ(η), for different values of S at γ=1: first solution (a) and second solution (b).
(a) (b)Finally, we end up our discussion with Figure10 which presents an overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the γ-S domain for Pr=0.7 and 0≤γ≤2. The straight line in this figure represents the occurrence of unique solution of the problem.Figure 10
An overview of the solution for problems (14) and (15) subject to the boundary conditions (16) in the γ-S domain.


## 7. 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,S=-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 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.



---

*Source: 8269541-2018-03-12.xml*