On the Method of Inverse Mapping for Solutions of Coupled Systems of Nonlinear Differential Equations Arising in Nanofluid Flow, Heat and Mass Transfer

Consider a steady, incompressible, laminar, two-dimensional boundary layer flow of a nanofluid at a vertical wall coincide with the plane y = 0, the flow being confined to y > 0 (see Figure 1). Two equal and opposite forces are introduced along the x-axis so that the wall is stretched while keeping the origin fixed. The sheet is then stretched with a velocity u_w = axⁿ where a is a constant, n is a nonlinear stretching parameter and x is the coordinate measured along the stretching surface. We make following assumptions:

Fig. 1

the pressure gradient and external forces are neglected

Under these assumptions, the basic equations for the conservation of mass, momentum, thermal energy and nanoparticles of the nanofluid can be written in Cartesian coordinates x and y as ( for details see Rana and Bhargava [1])

∂<!-- ? -->u∂<!-- ? -->x+∂<!-- ? -->v∂<!-- ? -->y=0,$$\begin{array}{} \displaystyle \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, \end{array}$$(1)

u∂<!-- ? -->u∂<!-- ? -->x+v∂<!-- ? -->u∂<!-- ? -->y=v∂<!-- ? -->2u∂<!-- ? -->y2,$$\begin{array}{} \displaystyle u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=v\frac{\partial^2 u}{\partial y^2}, \end{array}$$(2)

u∂<!-- ? -->T∂<!-- ? -->x+v∂<!-- ? -->T∂<!-- ? -->y=α<!-- a -->m∇<!-- ? -->2T+τ<!-- t -->DB∂<!-- ? -->C∂<!-- ? -->y∂<!-- ? -->T∂<!-- ? -->y+DTT∞<!-- 8 -->(∂<!-- ? -->T∂<!-- ? -->y)2,$$\begin{array}{} \displaystyle u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha_m \nabla^2 T+\tau\left[D_B\frac{\partial C}{\partial y}\frac{\partial T}{\partial y}+\frac{D_T}{T_{\infty}}(\frac{\partial T}{\partial y})^2\right], \end{array}$$(3)

u∂<!-- ? -->C∂<!-- ? -->x+v∂<!-- ? -->C∂<!-- ? -->y=DB∂<!-- ? -->2C∂<!-- ? -->y2+(DTT∞<!-- 8 -->)∂<!-- ? -->2T∂<!-- ? -->y2,$$\begin{array}{} \displaystyle u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D_B\frac{\partial^2 C}{\partial y^2}+(\frac{D_T}{T_{\infty}})\frac{\partial^2 T}{\partial y^2}, \end{array}$$(4)

where

α<!-- a -->m=km(ρ<!-- ? -->c)f,τ<!-- t -->=(ρ<!-- ? -->c)p(ρ<!-- ? -->c)f.$$\begin{array}{} \displaystyle \alpha_m=\frac{k_m}{(\rho c)_f},\quad \tau=\frac{(\rho c)_p}{(\rho c)_f}. \end{array}$$(5)

The boundary conditions for the problem are

v=0,uw=axn,T=Tw,C=Cwaty=0,$$\begin{array}{} \displaystyle v=0,\quad u_w=ax^n,\quad T=T_w,\quad C=C_w \quad \text{at}\quad y=0, \end{array}$$(6)

u=v=0,T=T∞<!-- 8 -->,C=C∞<!-- 8 -->asy→<!-- ? -->∞<!-- 8 -->.$$\begin{array}{} \displaystyle u=v=0,\quad T=T_{\infty},\quad C=C_{\infty}\quad \text{as}\quad y\rightarrow \infty. \end{array}$$(7)

Here u and v are the velocity in the x and y directions, ρ_f is the density of the base fluid, α_m is the thermal diffusivity, ν is the kinematic viscosity, a is a positive constant, D_B is the Brownian coefficient, D_T is the thermophoretic diffusion coefficient, τ is the ratio between the effective heat capacity of the nanoparticle material and heat capacity of the fluid, c is the volumetric volume expansion coefficient and ρ_p is the density of the nanoparticles.

Defining the new variables

η<!-- ? -->=ya(n+1)2vxn−<!-- - -->12,u=axnf′(η<!-- ? -->),v=−<!-- - -->av(n+1)2xn−<!-- - -->12(f+(n−<!-- - -->1n+1)η<!-- ? -->f′),$$\begin{array}{} \displaystyle \eta=y\sqrt{\frac{a(n+1)}{2v}}x^{\frac{n-1}{2}}, \quad u=ax^n f'(\eta),\quad v=-\sqrt{\frac{av(n+1)}{2}}x^{\frac{n-1}{2}}(f+(\frac{n-1}{n+1})\eta f'), \end{array}$$(8)

θ<!-- ? -->(η<!-- ? -->)=T−<!-- - -->T∞<!-- 8 -->Tw−<!-- - -->T∞<!-- 8 -->,ϕ<!-- ? -->(η<!-- ? -->)=C−<!-- - -->C∞<!-- 8 -->Cw−<!-- - -->C∞<!-- 8 -->,$$\begin{array}{} \displaystyle \theta(\eta)=\frac{T-T_{\infty}}{T_w-T_{\infty}},\quad \phi(\eta)=\frac{C-C_{\infty}}{C_w-C_{\infty}}, \end{array}$$(9)

and substituting in (1)-(4), we obtained

f‴+ff″−<!-- - -->2nn+1f′2=0,$$\begin{array}{} \displaystyle f'''+ff''-\left(\frac{2n}{n+1}\right)f'^2=0, \end{array}$$(10)

1Prθ<!-- ? -->″+fθ<!-- ? -->′+Nbθ<!-- ? -->′ϕ<!-- ? -->′+Nt(θ<!-- ? -->′)2=0,$$\begin{array}{} \displaystyle \frac{1}{Pr}\theta''+f\theta '+Nb\theta' \phi '+Nt(\theta ')^2=0, \end{array}$$(11)

ϕ<!-- ? -->″+12Lefϕ<!-- ? -->′+NtNbθ<!-- ? -->″=0,$$\begin{array}{} \displaystyle \phi ''+\frac{1}{2}Le f\phi '+\frac{Nt}{Nb}\theta ''=0, \end{array}$$(12)

with boundary conditions,

atη<!-- ? -->=0,f=0,f′=1,θ<!-- ? -->=1,ϕ<!-- ? -->=1,$$\begin{array}{} \displaystyle \text{ at} \quad \eta=0,\quad f=0,\quad f'=1,\quad \theta =1,\quad \phi=1, \end{array}$$(13)

asη<!-- ? -->→<!-- ? -->∞<!-- 8 -->,f′=0,θ<!-- ? -->=0,ϕ<!-- ? -->=0.$$\begin{array}{} \displaystyle \text{as}\quad \eta\rightarrow \infty, \quad f'=0,\quad \theta =0,\quad \phi=0. \end{array}$$(14)

The key thermophysical parameters are defined by:

Pr=vα<!-- a -->,Le=vDB,Nb=(ρ<!-- ? -->c)pDB(Cw−<!-- - -->C∞<!-- 8 -->)(ρ<!-- ? -->c)fv,Nt=(ρ<!-- ? -->c)pDT(Tw−<!-- - -->T∞<!-- 8 -->)(ρ<!-- ? -->c)fvT∞<!-- 8 -->.$$\begin{array}{} \displaystyle Pr=\frac{v}{\alpha},\quad Le=\frac{v}{D_B}, \quad Nb=\frac{(\rho c)_pD_B(C_w-C_{\infty})}{(\rho c)_fv}, \quad Nt=\frac{(\rho c)_pD_T(T_w-T_{\infty})}{(\rho c)_fvT_{\infty}}. \end{array}$$(15)

Here Pr,Le,Nb, and Nt denote the Prandtl number, the lewis number, the Brownian motion parameter and the thermophoresis parameter respectively.

In the present paper, we study the nonliner system analyticaly through the Optimam Homotopy Analysis Method by directly defining an inverse mapping 𝒯, i.e. without calculating any inverse operator. This method was intoduced by Liao [2] for a single differentiall equation. Vajravelu et al. [3] extended it to solve coupled systems. Here, we extend the method to a system of three nonlinear diferential equations using a common inverse linear mapping and approximated f(η), θ(η) and ϕ(η).

In this section, we discuss the set up of the problem using the details of OHAM ( see [4], [5] for more details) and MDDiM for the nonlinear system. First, we discuss the space that solution and base functions come from and then we derive deformation equations that we are trying to solve (nonlinear system). Finally, we use MDDiM to solve these deformation equations by introducing an appropriate inverse linear map 𝒯.

Define three nonlinear operators

N1[f(η<!-- ? -->),θ<!-- ? -->(η<!-- ? -->),ϕ<!-- ? -->(η<!-- ? -->)]=f‴+ff″−<!-- - -->2nn+1f′2,$$\begin{array}{} \displaystyle N_1[f(\eta),\theta(\eta), \phi(\eta)] =f'''+ff''-\left(\frac{2n}{n+1}\right)f'^2, \end{array}$$(16)

N2[f(η<!-- ? -->),θ<!-- ? -->(η<!-- ? -->),ϕ<!-- ? -->(η<!-- ? -->)]=1Prθ<!-- ? -->″+fθ<!-- ? -->′+Nbθ<!-- ? -->′ϕ<!-- ? -->′+Nt(θ<!-- ? -->′)2,$$\begin{array}{} \displaystyle N_2[f(\eta),\theta(\eta), \phi(\eta)] =\frac{1}{Pr}\theta''+f\theta '+Nb\theta' \phi '+Nt(\theta ')^2, \end{array}$$(17)

N3[f(η<!-- ? -->),θ<!-- ? -->(η<!-- ? -->),ϕ<!-- ? -->(η<!-- ? -->)]=ϕ<!-- ? -->″+12Lefϕ<!-- ? -->′+NtNbθ<!-- ? -->″=0,$$\begin{array}{} \displaystyle N_3[f(\eta),\theta(\eta), \phi(\eta)] =\phi ''+\frac{1}{2}Le f\phi '+\frac{Nt}{Nb}\theta ''=0, \end{array}$$(18)

so that N₁[f(η), θ(η), ϕ(η)] = 0, N₂[f(η), θ(η), ϕ(η)] = 0 and N₃[f(η), θ(η), ϕ(η)] = 0 give the original system (10)-(12). Take complete set of an infinite number of base functions that are linearly independent

S∞<!-- 8 -->=1,e−<!-- - -->δ<!-- d -->η<!-- ? -->,e−<!-- - -->2δ<!-- d -->η<!-- ? -->,…,$$\begin{array}{} \displaystyle S_{\infty} = \left\lbrace 1, e^{-\delta\eta}, e^{-2\delta\eta}, \ldots \right\rbrace , \end{array}$$(19)

and define the space of functions that is their linear combinations to be

V=∑<!-- ? -->k=0∞<!-- 8 -->ake−<!-- - -->kδ<!-- d -->η<!-- ? -->|ak∈<!-- ? -->R.$$\begin{array}{} \displaystyle V = \left\lbrace \sum_{k=0}^\infty a_k e^{-k \delta\eta} \big| a_k \in \mathbb{R} \right\rbrace . \end{array}$$(20)

That is, V is the solution and base space for f(η), θ(η) and ϕ(η).

Let

S∗<!-- * -->=1,e−<!-- - -->δ<!-- d -->η<!-- ? -->.$$\begin{array}{} \displaystyle S^* = \left\lbrace 1, e^{-\delta\eta}\right\rbrace . \end{array}$$(21)

denote a set, consists of first 2 members of S_∞. Next, form the space of functions taking their linear combinations

V∗<!-- * -->=a0+a1e−<!-- - -->δ<!-- d -->η<!-- ? -->|a0,a1∈<!-- ? -->R.$$\begin{array}{} \displaystyle V^* = \left\lbrace a_0 + a_1e^{-\delta\eta} | a_0, a_1\in \mathbb{R} \right\rbrace. \end{array}$$(22)

Then the primary solutions, or our initial guesses, μ(η) ∈ V^∗ have the form

μ(η<!-- ? -->)=∑<!-- ? -->j=01aje−<!-- - -->δ<!-- d -->η<!-- ? -->.$$\begin{array}{} \displaystyle \mu(\eta)=\sum_{j=0}^1 a_je^{-\delta\eta}. \end{array}$$(23)

Write

S^<!-- ^ -->=e−<!-- - -->2δ<!-- d -->η<!-- ? -->,e−<!-- - -->3δ<!-- d -->η<!-- ? -->,…,$$\begin{array}{} \displaystyle \widehat{S}= \left\lbrace e^{-2\delta\eta}, e^{-3\delta\eta}, \ldots \right\rbrace , \end{array}$$(24)

and define

V^<!-- ^ -->=∑<!-- ? -->k=2∞<!-- 8 -->ake−<!-- - -->kδ<!-- d -->η<!-- ? -->|ak∈<!-- ? -->R.$$\begin{array}{} \displaystyle \widehat{V} = \left\lbrace \sum_{k=2}^\infty a_ke^{-k\delta\eta}\big| a_k \in \mathbb{R} \right\rbrace . \end{array}$$(25)

Obviously, V = V̂ ∪ V^∗.

Next, define

SR=ψ<!-- ? -->1(η<!-- ? -->),ψ<!-- ? -->2(η<!-- ? -->),…,$$\begin{array}{} \displaystyle S_R = \left\lbrace \psi_1(\eta), \psi_2(\eta), \ldots \right\rbrace , \end{array}$$(26)

which is an infinite set of base functions that are linearly independent, and set of linear combinations of functions from S_R

U=∑<!-- ? -->k=1∞<!-- 8 -->ckψ<!-- ? -->k(η<!-- ? -->)|ck∈<!-- ? -->R.$$\begin{array}{} \displaystyle U = \left\lbrace \sum_{k=1}^\infty c_k \psi _k (\eta) \big| c_k \in \mathbb{R} \right\rbrace . \end{array}$$(27)

Assuming that N₁[f(η), θ(η), ϕ(η)], N₂[f(η), θ(η), ϕ(η)], N₃[f(η), θ(η), ϕ(η)] ∈ U, then N₁, N₂, N₃ : V → U.

Optimal Homotopy Analysis Method allows us to obtain approximate series solutions to wide variety of nonlinear systems. Define three homotopies of operators ℋ₁, ℋ₂ and ℋ₃

0≡<!-- = -->H1(f,θ<!-- ? -->,ϕ<!-- ? -->,q)=(1−<!-- - -->q)L1[f]−<!-- - -->c0qN1[f,θ<!-- ? -->,ϕ<!-- ? -->],$$\begin{array}{} \displaystyle 0 \equiv \mathscr{H}_1(f,\theta,\phi, q) = (1-q)L_1[f] - c_0qN_1[f,\theta, \phi], \end{array}$$(28)

0≡<!-- = -->H2(f,θ<!-- ? -->,ϕ<!-- ? -->,q)=(1−<!-- - -->q)L2[θ<!-- ? -->]−<!-- - -->c1qN2[f,θ<!-- ? -->,ϕ<!-- ? -->],$$\begin{array}{} \displaystyle 0 \equiv \mathscr{H}_2(f,\theta,\phi, q) = (1-q)L_2[\theta] - c_1qN_2[f,\theta, \phi], \end{array}$$(29)

0≡<!-- = -->H3(f,θ<!-- ? -->,ϕ<!-- ? -->,q)=(1−<!-- - -->q)L3[ϕ<!-- ? -->]−<!-- - -->c2qN3[f,θ<!-- ? -->,ϕ<!-- ? -->],$$\begin{array}{} \displaystyle 0 \equiv \mathscr{H}_3(f,\theta,\phi, q) = (1-q)L_3[\phi] - c_2qN_3[f,\theta, \phi], \end{array}$$(30)

through the homotopy embedding parameter q ∈ [0, 1], between nonlinear operators N₁, N₂, N₃ and an auxiliary linear operators L₁, L₂, L₃. Here, c₀, c₁, c₂ ≠ 0 are the converge control parameters which will be used to optimize the function approximations in the next section. In the frame of OHAM, the series solution of f, θ and ϕ is given by

f(η<!-- ? -->)=f0(η<!-- ? -->)+∑<!-- ? -->k=1∞<!-- 8 -->fk(η<!-- ? -->)qk,$$\begin{array}{} \displaystyle f(\eta) = f_0(\eta) + \sum_{k=1}^\infty f_k(\eta)q^k, \end{array}$$(31)

θ<!-- ? -->(η<!-- ? -->)=θ<!-- ? -->0(η<!-- ? -->)+∑<!-- ? -->k=1∞<!-- 8 -->θ<!-- ? -->k(η<!-- ? -->)qk,$$\begin{array}{} \displaystyle \theta(\eta) = \theta_0(\eta) + \sum_{k=1}^\infty \theta_k(\eta)q^k, \end{array}$$(32)

ϕ<!-- ? -->(η<!-- ? -->)=ϕ<!-- ? -->0(η<!-- ? -->)+∑<!-- ? -->k=1∞<!-- 8 -->ϕ<!-- ? -->k(η<!-- ? -->)qk,$$\begin{array}{} \displaystyle \phi(\eta) = \phi_0(\eta) + \sum_{k=1}^\infty \phi_k(\eta)q^k, \end{array}$$(33)

where f₀(η), θ₀(η) and ϕ₀(η) are initial guesses that satisfy boundary conditions (13)-(14) and belong to the set V.

It is clear that when q = 0 in the homotopies (28)-(30), they become L₁[f] = 0, L₂[θ] = 0 and L₃[ϕ] = 0; but for q = 1, the original nonlinear differential equations N₁[f, θ, ϕ] = 0, N₂[f, θ, ϕ] = 0 and N₃[f, θ, ϕ] = 0 are recovered. In addition, when q = 1 in the expansions (31)-(33), the solutions f, θ and ϕ are a sum of the components f₀, f₁, f₂, …, θ₀, θ₁, θ₂, … and ϕ₀, ϕ₁, ϕ₂, …. Substituting (31)-(33) in to the first homotopy (28), we get the deformation equations

L1[f0(η<!-- ? -->)]=0,f0(0)=0,f0′(0)=1f0′→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->,$$\begin{array}{} \displaystyle L_1[f_0(\eta)] = 0, \quad f_0(0)=0,\quad f_0'(0)=1 \quad f_0' \to 0 \text{ as } \eta \to \infty, \end{array}$$(34)

and for k ≥ 1 we have

L1[fk(η<!-- ? -->)]=χ<!-- ? -->kL1[fk−<!-- - -->1(η<!-- ? -->)]+c0Dk−<!-- - -->11(η<!-- ? -->),fk(0)=0,fk′(0)=0,fk′→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->,$$\begin{array}{} \displaystyle L_1[f_k(\eta)] = \chi_k L_1[f_{k-1}(\eta)]+c_0 \mathscr{D}_{k-1}^1(\eta), \quad f_k(0)=0,\quad f_k'(0)=0, \quad f_k' \to 0 \text{ as } \eta \to \infty, \end{array}$$(35)

where

χ<!-- ? -->k=0,k≤<!-- = -->1,1,k≥<!-- = -->1.$$\begin{array}{} \displaystyle \chi_k = \begin{cases} 0, \quad k \leq 1, \\ 1, \quad k \geq 1. \end{cases} \end{array}$$(36)

Here Dkξ<!-- ? -->$\begin{array}{} \displaystyle \mathscr{D}_k^\xi \end{array}$, for ξ = 1, 2, 3, is the homotopy derivative defined to be

Dk−<!-- - -->1ξ<!-- ? -->(η<!-- ? -->)=1(k−<!-- - -->1)!∂<!-- ? -->k−<!-- - -->1∂<!-- ? -->qk−<!-- - -->1Nξ<!-- ? -->∑<!-- ? -->j=0∞<!-- 8 -->fj(η<!-- ? -->)qj,∑<!-- ? -->j=0∞<!-- 8 -->θ<!-- ? -->j(η<!-- ? -->)qj|q=0.$$\begin{array}{} \displaystyle \mathscr{D}_{k-1}^\xi(\eta) = \frac{1}{(k-1)!}\left(\frac{\partial^{k-1}}{\partial q^{k-1}} N_\xi \left[\sum_{j=0}^\infty f_j(\eta)q^j, \sum_{j=0}^\infty \theta_j (\eta)q^j \right]\right) \Big| _{q=0}. \end{array}$$(37)

Similarly, substituting (31)-(33) into (29) and (30) obtained:

L2[θ<!-- ? -->0(η<!-- ? -->)]=0,θ<!-- ? -->0(0)=1,θ<!-- ? -->0→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->,$$\begin{array}{} \displaystyle L_2[\theta_0(\eta)] = 0, \quad \theta_0 (0) = 1, \quad \theta_0 \to 0 \text{ as } \eta \to \infty , \end{array}$$(38)

L3[ϕ<!-- ? -->0(η<!-- ? -->)]=0,ϕ<!-- ? -->0(0)=1,ϕ<!-- ? -->0→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->,$$\begin{array}{} \displaystyle L_3[\phi_0(\eta)] = 0, \quad \phi_0 (0) = 1, \quad \phi_0 \to 0 \text{ as } \eta \to \infty , \end{array}$$(39)

and for k ≥ 1

L2[θ<!-- ? -->k(η<!-- ? -->)]=χ<!-- ? -->kL2[θ<!-- ? -->k−<!-- - -->1(η<!-- ? -->)]+c1Dk−<!-- - -->12(η<!-- ? -->),θ<!-- ? -->k(0)=0,θ<!-- ? -->k→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->.$$\begin{array}{} \displaystyle L_2[\theta_k(\eta)] = \chi_k L_2[\theta_{k-1}(\eta)]+c_1 \mathscr{D}_{k-1}^2(\eta), \quad \theta_k (0)=0, \quad \theta_k \to 0 \text{ as } \eta \to \infty. \end{array}$$(40)

L3[ϕ<!-- ? -->k(η<!-- ? -->)]=χ<!-- ? -->kL3[ϕ<!-- ? -->k−<!-- - -->1(η<!-- ? -->)]+c1Dk−<!-- - -->13(η<!-- ? -->),ϕ<!-- ? -->k(0)=0,ϕ<!-- ? -->k→<!-- ? -->0 as η<!-- ? -->→<!-- ? -->∞<!-- 8 -->.$$\begin{array}{} \displaystyle L_3[\phi_k(\eta)] = \chi_k L_3[\phi_{k-1}(\eta)]+c_1\mathscr{D}_{k-1}^3(\eta), \quad \phi_k (0)=0, \quad \phi_k \to 0 \text{ as } \eta \to \infty. \end{array}$$(41)

Using Liao’s Method of Directly Defined Inverses, the deformation equations (35) and (40)-(41) are

fk(η<!-- ? -->)=χ<!-- ? -->kfk−<!-- - -->1(η<!-- ? -->)+c0JDk−<!-- - -->11(η<!-- ? -->)+ak,1e−<!-- - -->δ<!-- d -->η<!-- ? -->+ak,0,$$\begin{array}{} \displaystyle f_k(\eta) = \chi_k f_{k-1}(\eta) + c_0 \mathscr{J}\left[\mathscr{D}_{k-1}^1(\eta)\right] + a_{k,1}e^{-\delta\eta}+a_{k,0}, \end{array}$$(42)

θ<!-- ? -->k(η<!-- ? -->)=χ<!-- ? -->kθ<!-- ? -->k−<!-- - -->1(η<!-- ? -->)+c1JDk−<!-- - -->12(η<!-- ? -->)+bk,1e−<!-- - -->δ<!-- d -->η<!-- ? -->+bk,0,$$\begin{array}{} \displaystyle \theta_k(\eta) = \chi_k \theta_{k-1}(\eta) + c_1\mathscr{J}\left[\mathscr{D}_{k-1}^2(\eta)\right] + b_{k,1}e^{-\delta\eta}+b_{k,0}, \end{array}$$(43)

ϕ<!-- ? -->k(η<!-- ? -->)=χ<!-- ? -->kϕ<!-- ? -->k−<!-- - -->1(η<!-- ? -->)+c1JDk−<!-- - -->13(η<!-- ? -->)+ck,1e−<!-- - -->δ<!-- d -->η<!-- ? -->+ck,0.$$\begin{array}{} \displaystyle \phi_k(\eta) = \chi_k \phi_{k-1}(\eta) + c_1\mathscr{J}\left[\mathscr{D}_{k-1}^3(\eta)\right] + c_{k,1}e^{-\delta\eta}+c_{k,0}. \end{array}$$(44)

The benifit of the Optimal Homotopy Analysis Mehod is that it has a great freedom to choose the auxillary linear operators L₁, L₂ and L₃ and initial guesses f₀(η), θ₀(η), ϕ₀(η). After auxillary linear operator and initial guesses are properly choosen we are free to determine how many terms f_k, θ_k, ϕ_k ∈ V we want and can do iteratively. There has been great success in solving systems of nonlinear differential equations using OHAM (see [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17]).

The only drawback of this homotopy analysis method is spending a lot of CPU time. First we choose auxiliary linear operators, and then solving the linear higher order deformation equation only to find out the inverse operators and applying them to quickly-growing expressions. However, in the latest innovation of Liao we have the freedom to directly define inverse operator by completely neglecting the linear operator. So, using this novel method can solve higher order deformation equations quickly and it’s unnecessary to calculate inverse linear operators.

In our work the inversely defined mapping, 𝒯, is the same for all three equations. But different directly defined inverses could be chosen if a different structure for the solutions f, θ and ϕ is required.

Define 𝒯 : U → V by

Je−<!-- - -->kδ<!-- d -->η<!-- ? -->=e−<!-- - -->kδ<!-- d -->η<!-- ? -->Ak3+k,$$\begin{array}{} \displaystyle \mathscr{J}\left[e^{-k \delta\eta}\right] = \frac{e^{-k\delta \eta }}{Ak^3+k }, \end{array}$$(45)

where A, δ are parameters which will be used to optimize the square residual error functions.

The appropriate solutions for the system (10)-(12) with boundary conditions (13)-(14) are obtained using MDDiM. Further used error analysis to get a general idea of how good the approximations are.

Define three term approximation f̂, θ̂ and ϕ̂ which is sum of the first three solutions to the deformation equations. If they are exact, then they solve system (10)-(12), i.e., if N₁[f̂, θ̂, ϕ̂] = 0, N₂[f̂, θ̂, ϕ̂] = 0 and N₃[f̂, θ̂, ϕ̂] = 0, then the three term approximations are exact solutions. If not N₁[f̂, θ̂, ϕ̂], N₂[f̂, θ̂, ϕ̂] and N₃[f̂, θ̂, ϕ̂] become residual error functions that can be evaluated at any point η in the domain of the problem. Taking square of the L²-norm of error functions and setting converge control parameters to be c₀ = c₂ = c₃ define square residual error functions

Eξ<!-- ? -->(Le,Nb,Pr,Nt,n,A,c0,δ<!-- d -->)=∫<!-- ? -->0∞<!-- 8 -->Nξ<!-- ? -->f^<!-- ^ -->(η<!-- ? -->),θ<!-- ? -->^<!-- ^ -->(η<!-- ? -->),ϕ<!-- ? -->^<!-- ^ -->(η<!-- ? -->)2dη<!-- ? -->,$$\begin{array}{} \displaystyle E_{\xi}(Le, Nb, Pr, Nt, n,A, c_0, \delta) = \int_0^\infty \left(N_\xi\left[\widehat{f}(\eta), \widehat{\theta}(\eta), \widehat{\phi}(\eta)\right]\right)^2 d\eta , \end{array}$$(46)

for ξ = 1, 2, 3. Since we have three error functions we will take affine combination of them as

E(Le,Nb,Pr,Nt,n,A,c0,δ<!-- d -->)=∑<!-- ? -->ξ<!-- ? -->=13Eξ<!-- ? -->(Le,Nb,Pr,Nt,n,A,c0,δ<!-- d -->).$$\begin{array}{} \displaystyle E(Le, Nb, Pr, Nt, n,A, c_0, \delta) = \sum_{\xi=1}^3E_{\xi}(Le, Nb, Pr, Nt, n,A, c_0, \delta). \end{array}$$(47)

But in practice the evaluation of E_ξ(Le,Nb,Pr,Nt,n,A,c₀,δ) is much time consuming so instead of exact residual error we use average residual error defined as

E^<!-- ^ -->ξ<!-- ? -->(Le,Nb,Pr,Nt,n,A,c0,δ<!-- d -->)=1M+1∑<!-- ? -->j=0MNξ<!-- ? -->f^<!-- ^ -->(j),θ<!-- ? -->^<!-- ^ -->(j),ϕ<!-- ? -->^<!-- ^ -->(j)2.$$\begin{array}{} \displaystyle \widehat{E}_{\xi}(Le, Nb, Pr, Nt, n,A, c_0, \delta) =\frac{1}{M+1} \sum_{j=0}^M \left(N_\xi\left[\widehat{f}(j), \widehat{\theta}(j), \widehat{\phi}(j)\right]\right)^2. \end{array}$$(48)

Now, we minimize error functions with respect to A, c₀, δ and obtained optimal values of A, c₀, δ. Substituting those values in to f̂, θ̂ and ϕ̂ we get three term approximation solution to the system (10)-(12) which satisfies the conditions (13)-(14).

We start with initial guesses f₀(η), θ₀(η) and ϕ₀(η) that satisfy the boundary conditions (13)-(14), respectively. We choose

f0(η<!-- ? -->)=1δ<!-- d -->−<!-- - -->1δ<!-- d -->e−<!-- - -->δ<!-- d -->η<!-- ? -->,$$\begin{array}{} \displaystyle f_0(\eta)=\frac{1}{\delta}-\frac{1}{\delta}e^{-\delta\eta}, \end{array}$$(49)

θ<!-- ? -->0(η<!-- ? -->)=e−<!-- - -->δ<!-- d -->η<!-- ? -->,$$\begin{array}{} \displaystyle \theta_0(\eta)=e^{-\delta\eta}, \end{array}$$(50)

and

ϕ<!-- ? -->0(η<!-- ? -->)=e−<!-- - -->δ<!-- d -->η<!-- ? -->.$$\begin{array}{} \displaystyle \phi_0(\eta)=e^{-\delta\eta}. \end{array}$$(51)

Now, using the deformation equations (42)-(44) to find f₁(η), θ₁(η) and ϕ₁(η), they are

f1(η<!-- ? -->)=−<!-- - -->12c0(n−<!-- - -->1)(n+1)(4A+1)+c0(n−<!-- - -->1)(n+1)(4A+1)e−<!-- - -->δ<!-- d -->η<!-- ? -->−<!-- - -->12c0(n−<!-- - -->1)(n+1)(4A+1)e−<!-- - -->2δ<!-- d -->η<!-- ? -->,$$\begin{array}{} \displaystyle f_1(\eta)=-\frac{1}{2}\frac{c_0(n-1)}{(n+1)(4A+1)}+\frac{c_0(n-1)}{(n+1)(4A+1)}e^{-\delta\eta}-\frac{1}{2}\frac{c_0(n-1)}{(n+1)(4A+1)}e^{-2\delta\eta}, \end{array}$$(52)

θ<!-- ? -->1(η<!-- ? -->)−<!-- - -->12(1+Nt⋅δ<!-- d -->2+Nb⋅δ<!-- d -->2)c0⋅δ<!-- d -->24A+1e−<!-- - -->δ<!-- d -->η<!-- ? -->+12(1+Nt⋅δ<!-- d -->2+Nb⋅δ<!-- d -->2)c0⋅δ<!-- d -->24A+1e−<!-- - -->2δ<!-- d -->η<!-- ? -->,$$\begin{array}{} \displaystyle \theta_1(\eta)-\frac{1}{2}\frac{(1+Nt\cdot\delta^2+Nb \cdot\delta^2)c_0\cdot\delta^2}{4A+1}e^{-\delta\eta}+\frac{1}{2}\frac{(1+Nt\cdot\delta^2+Nb \cdot\delta^2)c_0\cdot\delta^2}{4A+1}e^{-2\delta\eta}, \end{array}$$(53)

and

ϕ<!-- ? -->1(η<!-- ? -->)=−<!-- - -->14⋅Le⋅c04A+1e−<!-- - -->δ<!-- d -->η<!-- ? -->+14⋅Le⋅c04A+1e−<!-- - -->2δ<!-- d -->η<!-- ? -->.$$\begin{array}{} \displaystyle \phi_1(\eta)=-\frac{1}{4}\cdot\frac{Le \cdot c_0}{4A+1}e^{-\delta\eta}+\frac{1}{4}\cdot\frac{Le\cdot c_0}{4A+1}e^{-2\delta\eta}. \end{array}$$(54)

Using only three terms, let f̂(η) = f₀(η)+f₁(η)+f₂(η), θ̂(η) = θ₀(η)+θ₁(η)+θ₂(η) and ϕ̂(η) = ϕ₀(η)+ϕ₁(η)+ϕ₂(η), the sum of the square residual error function is given by

E(A,c0,δ<!-- d -->)=1500∑<!-- ? -->j=0499∑<!-- ? -->ξ<!-- ? -->=13Nξ<!-- ? -->f^<!-- ^ -->(j),θ<!-- ? -->^<!-- ^ -->(j),ϕ<!-- ? -->^<!-- ^ -->(j)2.$$\begin{array}{} \displaystyle E(A, c_0, \delta) =\frac{1}{500} \sum_{j=0}^{499} \left(\sum_{\xi=1}^3\left(N_\xi\left[\widehat{f}(j), \widehat{\theta}(j), \widehat{\phi}(j)\right]\right)\right)^2. \end{array}$$(55)

and it is a function of A, c₀ and δ with parameters Le,Nb,Pr,Nt and n in it.

Using three different sets of values for the parameters Le,Nb,Pr,Nt and n we found the sum of the square residual error E(A,c₀,δ) and are presented below.

The plot of the error functions E(A,c₀,δ) is given in Figures 2-4 for three schemes at their optimum A values.

Fig. 2

Fig. 3

Fig. 4

The plots of f̂(η) and f̂′(η) are presented in Figures 5-6, for parametric values in Table 1 for E₁(A,c₀,δ). In Figures 7-8 the plots of θ̂(η) and ϕ̂(η) are presented for parametric values in Table 1 for E(A,c₀,δ).

Fig. 5

Fig. 6

Fig. 7

Fig. 8

Table 1

Minimum of the squared residual error E(A,c₀,δ) for three different sets of parameters.

A very good validation of the present analytical results has been achieved with the numerical results as shown in Figure 9. Also, it is found that the squared residual error decreases as a function of the number of terms in the approximation series, as shown in Figure 10.

Fig. 9

Fig. 10

The skin friction at the surface |–f̂″(0)| as a function of the stretching parameter n is presented in Figure 11. It is found that skin friction decreases with an increase in stretching parameter. Figure 12 illustrated Nusselt number |–θ̂′(0)| as a function of Lewis number (Le) and Brownian motion parameter (Nb). It is found that Nusselt number decreases with increase Nt and Nb. Figure 13 illustrated Sherwood number |–ϕ̂′(0)| as a function of Nt, Nb and it is found that Sherwood number increases with increase Nt and but decreases with increasing Nb.

Fig. 11

Fig. 12

Fig. 13

Liao’s Directly Defining inverse Mapping method is extended to a system of three nonlinear diferential equations. Approximate series solutions for f(η), θ(η), and ϕ(η) are obtained. Also, illustrated dimensionless velocity (f(η)), dimentionless temperature (θ(η)) and dimensionless concentration (ϕ(η)) profiles for three set of parameters (see Figures 5-8) are presented. Further, analytical results are compared with the numerical results (see Figure 9) and studied convergence of analytical results (see Figure 10).

Since the inverse operator is directly defined, the series solutions are obtained with less CPU time. The freedom of choosing the inverse operator leads to obtaining less complicated terms for the approximation solution. Futher, the selected inverse linear operator leads to three term solution which is accurate up to five decimal places by optimizing square residual function with respect to A, δ, and c₀. Hence, we can conlude that MDDiM is not only easy to use, but also accurate. Theoretically, even if a smaller error was desired, it would just amount to computing more terms in the series by solving deformation equations. Furthermore, one can write an algorithm to iteration approach and truncate the approximate series solution at a given accuracy.

The idea is novel and is useful. This idea is not limited to a single nonlinear differential equation, but can be used for system of several equations. Also, it is important to note that finding an inverse operator that works well for the equation and it leads to an easily generated solution series. Hence, it is worth-while to investigate this inverse linear operator.