# Approximate Solution of Nonlinear System of BVP Arising in Fluid Flow Problem

Authors: A. K. Alomari; N. Ratib Anakira; A. Sami Bataineh; I. Hashim
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We extend for the first time the applicability of the Optimal Homotopy Asymptotic Method (OHAM) to find approximate solution of a system of two-point boundary-value problems (BVPs). The OHAM provides us with a very simple way to control and adjust the convergence of the series solution using the auxiliary constants which are optimally determined. Comparisons made show the effectiveness and reliability of the method.


---

## Body

## 1. Introduction

Many real-world problems can be modelled by nonlinear differential equations. For example, fluid flow problems can give rise to boundary-value problems (BVPs) or systems of BVPs with conditions specified at two or more different points. Finding a reliable method for solving BVPs is of great interest. Noor and Mohyud-Din [1–3] presented approximate solutions of some classes of BVPs by using the variational iteration method (VIM), homotopy perturbation method (HPM), and variational iteration decomposition method (VIDM). Herisanu et al. [4] developed the so-called Optimal Homotopy Asymptotic Method (OHAM) for solving nonlinear problems. OHAM provides us with a  very simple way to control and adjust the convergence of the series solution using the auxiliary constants which are optimally determined. Several promising applications of OHAM to problems in fluid dynamics have been presented [5–12]. Ali et al. [13, 14] solved several two-point and multipoint BVPs by OHAM. Very recently, Hashmi et al. [15] applied OHAM for finding the approximate solutions of a class of Volterra integral equations with weakly singular kernels.The laminar fully developed combined free and forced magnetoconvection in a vertical channel with symmetric and asymmetric boundary heatings in the presence of viscous and Joulean dissipations was studied by  Umavathi and Malashetty [16]. The mathematical model describing the channel flow problem is governed by a system of nonlinear BVPs. Umavathi and Malashetty [16] employed the classical perturbation technique to solve the system of BVPs. The aim of the present work is thus to propose an accurate approach to the channel flow problem using an analytical technique, namely, OHAM. The efficiency of the procedure is based on the construction and determination of the auxiliary functions combined with a convenient way to optimally control the convergence of the solution.


## 2. The Model Equation

The system of BVPs modelling the channel flow problem as given in [16] is
(1)d4udy4-M2d2udy2=M2GRBru2+GRBr(dudy)2,d2θdy2=-M2Bru2-Br(dudy)2,
subject to
(2)u(-14)=u(14)=0,d2udy2∣y=-1/4=-48+RTGR2,d2udy2∣y=1/4=-48-RTGR2,θ(-14)=-RT2,θ(14)=RT2,
where the parameter RT becomes one for asymmetric heating and zero for symmetric heating. The special case Br=0 was solved exactly by Umavathi and Malashetty [16], and the exact solutions are
(3)u=48M2(1-cosh(My)cosh(M/4))+2GRRTM2(y-sinh(My)4sinh(M/4)),θ=2RTy.
Furthermore, when GR=0, solutions of (1)-(2) become
(4)u=48M2(1-cosh(My)cosh(M/4)),θ=A(y2-116)+B[cosh(2My)-cosh(M2)]+C[cosh(My)-cosh(M4)]+2RTy,
where
(5)A=-1152BrM2,B=-576BrM4cosh2(M/4),C=-4608BrM4cosh(M/4).
We remark that the general case of both Br≠0 and GR≠0 is very difficult to solve exactly. For this case, Umavathi and Malashetty [16] have given the standard perturbation solutions by assumingε=BrGRto be the small parameter in the expansion.


## 3. Basic Idea of OHAM

Consider the following differential equations:(6)L(u(y))+g(y)+N(u(y))=0,B(u,dudy)=0,
where L is a linear operator, N is a nonlinear operator, u(y) is an unknown function, y denotes independent variable, g(y) is a known function, and B is a boundary operator.According to the basic idea of OHAM [4–6], we construct a homotopy h(v(y,p),p):R×[0,1]→R which satisfies
(7)(1-p)[L(v(y,p))+g(y)]=H(p)[L(v(y,p))+g(y)+N(v(y,p))],(8)B(v(y,p),∂v(y,p)∂y)=0,
where y∈R and p∈[0,1] is an embedding parameter, H(p) is a nonzero auxiliary function for p≠0, H(0)=0 and v(y,p) is an unknown function. Obviously, when p=0 and p=1 it holds that v(y,0)=u0(y) and v(y,1)=u(y), respectively. Thus, as  p varies from 0 to 1, the solution v(y,p) approaches from u(y) to u0(y),where u0(y) is obtained from (7) for p=0, and we have
(9)L(u0(y))+g(y)=0,B(u0,du0dy)=0.
Next, we choose auxiliary function H(p) in the form
(10)H(p)=pC1+p2C2+p3C3+⋯,
where C1,C2,C3,… are constants to be determined, and  H(p) can be expressed in many forms as reported in [4–7].To get an approximate solution, we expandv(y,p,Ci) in Taylor’s series about p in the following manner:
(11)v(y,p,Ci)=u0(y)+∑k=1∞‍uk(y,C1,C2,…,Ck)pk.
Substituting (11) into (7) and equating the coefficient of the like powers of p, we obtain the following linear equations. The zeroth-order problem is given by (9), and  the first- and second-order problems are given as
(12)L(u1(y))=C1N0(u0(y)),B(u1,du1dy)=0,L(u2(y))-L(u1(y))=C2N0(u0(y))+C1×[Lu1(y)+N1(u0(y),u1(y))],B(u2,du2dy)=0.
And the general governing equations for uk(y) are given as
(13)L(uk(y))-L(uk-1(y))=CkN0(u0(y))+∑i=1k-1‍Ci[L(uk-i(y))+Nk-i(u0(y),u1(y),…,uk-1(y))uuu1u1(y),…,uk-1(y))],B(uk,dukdy)=0,
where k=2,3,… and Nm(u0(y),u1(y),…,um(y)) is the coefficient of pm in the expansion of N(v(y,p)) about the embedding parameter p(14)N(v(y,p,Ci))=N0(u0(y))+∑m=1∞‍Nm(u0(y),u1(y),…,um(y))pm.
It has been observed that the convergence of the series (11) depends upon the auxiliary constants C1,C2,C3,…. If the series is convergent at  p=1,  one has
(15)v(y,Ci)=u0(y)+∑k=1∞‍uk(y,C1,C2,…,Ck).
The results of the mth-order approximations are
(16)u~(y,C1,C2,C3,…,Cm)=u0(y)+∑i=1m‍ui(y,C1,C2,…,Ci).
Substituting (16) into (6) it results the following residual:
(17)R(y,C1,C2,C3,…,Cm)=L(u~(y,C1,C2,C3,…,Cm))+g(y)+N(u~(y,C1,C2,C3,…,Cm)).
If R=0, then u~ will be the exact solution. Generally this does not happen, especially in nonlinear problems. In order to find the optimal values of Ci, i=1,2,3,…, we first construct the functional
(18)J(C1,C2,C3,…,Cm)=∫ab‍R2(y,C1,C2,C3,…,Cm)dy,
and then minimizing it, we have
(19)∂J∂C1=∂J∂C2=⋯=∂J∂Cm=0,
where a and b are in the domain of the problem. With these constants known, the approximate solution (of order m) is well determined.


### 3.1. Application of OHAM

In this section, we apply OHAM for solving the nonlinear system of two-point BVP (1)-(2). By applying the proposed method, the zeroth-order deformation equation is
(20)(1-p)L[u~(y,p)-u0(y)]=H(p,Ci)[N(u~(y,p))],(1-p)L[θ~(y,p)-θ0(y)]=H(p,Cj)[N(θ~(y,p),u~(y,p))],
subject to the boundary conditions
(21)u~(-14,p)=0,u~(14,p)=0,u~′′(-14,p)=-48+RTGR2,u~′′(14,p)=-48-RTGR2,θ~(-14,y)=-RT2,θ~(14,y)=RT2.
Using the framework of OHAM the mth-order(22)L1[u~-χmum-1]=R1,m(Ci,u→m-1),
where u→m-1={u0,u1,…,um-1},
(23)R1,m=∑i=0m-1‍Ci+1um-1-i(4)-M2GRBr∑i=0m-1‍um-1-i×∑j=0i‍Cj+1uj-i-GRBr∑i=0m-1‍um-1-i′∑j=0i‍Cj+1uj-1′,L2[θ~m-χmθm-1]=R2,m(Ci,θ→m-1),
where
(24)R2,m=∑i=0m-1‍Ci+1θm-1-i′′+M2Br∑i=0m-1‍um-1-i×∑j=0i‍Ci+1uj-i+Br∑i=0m-1‍um-1-i′∑j=0i‍Ci+1uj-1′.
Now the zeroth-order problem is
(25)u0(4)(y)=0,θ0′′(y)=0,
subject to the boundary conditions
(26)u0(-14)=0,u0(14)=0,u0′′(-14)=-48+RTGR2,u0′′(14)=-48-RTGR2,θ0(-14)=-RT2,θ0(14)=RT2.
The solutions are
(27)u0(y)=148(72-1152y2+yRTGR-16y3RTGR),θ0(y)=2yRT.
Now the first-order problem is
(28)u1(4)(y,C1)=-C1M2BrGRu02(y)-C1BrGR(u0′(y))2-C1M2u0′′(y)+C1u0(4)(y),θ1′′(y,C1)=C1M2Bru02(y)+C1Br(u0′(y))2+C1θ0′′(y),
subject to the boundary conditions
(29)u1(-14)=0,u1(14)=0,u1′′(-14)=0,u1′′(14)=0,θ1(-14)=0,θ1(14)=0.
The second-order problem is
(30)u2(4)(y,C1,C2)=-C2M2u0′′(y)+C2u0(4)(y)+u1(4)(y,C1)-C2M2BrGRu02(y)-2C1M2BrGRu0(y)u1(y,C1)-C2BrGR(u0′(y))2+C1u1(4)(y,C1)-C1M2u1′′(y,C1)-2C1BrGRu0′(y)u1′(y,C1),θ2′′(y,C1,C2)=C2M2Bru02(y)+2C1M2Bru0(y)u1(y,C1)+C2Br(u0′(y))2+2c1Bru0′(y)u1′(y,C1)+C2θ0′′(y)+(1+C1)θ1′′(y,C1),
subject to the boundary conditions
(31)u2(-14)=0,u2(14)=0,u2′′(-14)=0,u2′′(14)=0,θ2(-14)=0,θ2(14)=0.
The third-order problem  is
(32)u3(4)(y,C1,C2,C3)=-M2BrC3GRu02(y)-2M2BrC2GRu0(y)u1(y,C1)-C1M2BrGRu12(y,C1)+C3u0(4)(y)+C2u1(4)(y,C1)-C3BrGR(u0′(y))2-2C2BrGRu0′(y)-u1′(y,C1)-C1BrGR(u1′(y,C1))2-C1M2u2′′(y,C1,C2)-C3M2u0′′(y)-C2M2u1′′(y,C1)+C1u2(4)(y,C1,C2)-2C1M2BrGRu0(y)u2(y,C1,C2)+u2(4)(y,C1,C2)-2C1BrGRu0′(y)u2′(y,C1,C2),θ3′′(y,C1,C2,C3)=C3M2Bru02(y)+2C2M2Bru0(y)u1(y,C1)+C1M2Bru12(y,C1)+2C1M2Bru0(y)u2(y,C1,C2)+C3Br(u0′(y))2+2C2Bru0′(y)u1′(y,C1)+C1Br(u1′(y,C1))2+2C1Bru0′(y)u2′(y,C1,C2)+C3θ0′′(y)+C2θ1′′(y,C1)+θ2′′(y,C1,C2)+C1θ2′′(y,C1,C2),
subject to the boundary conditions
(33)u3(-14)=0,u3(14)=0,u3′′(-14)=0,u3′′(14)=0,θ3(-14)=0,θ3(14)=0.Using the solution of (25)–(32) we obtain the following four-term approximate solutions for u and θ by OHAM takingp=1:
(34)u~(y,C1,C2,C3)=u0(y)+u1(y,C1)+u2(y,C1,C2)+u3(y,C1,C2,C3),θ~(y,C1,C2,C3)=θ0(y)+θ1(y,C1)+θ2(y,C1,C2)+θ3(y,C1,C2,C3).
The explicit expressions for the individual terms of the approximate solutions are not given here for brevity. Taking the residual errors
(35)Ru~(y,C1,C2,C3)=u~(4)(y,C1,C2,C3)-M2u~′′(y,C1,C2,C3)-M2GRBru~(y,C1,C2,C3)-GRBr(u~′(y,C1,C2,C3))2,Rθ~(y,C1,C2,C3)=θ~′′(y,C1,C2,C3)+M2Br(u~(y,C1,C2,C3))2-Br(u~′(y,C1,C2,C3))2,
the optimal values of Ci’s can be obtained. Table 1 shows some optimal values of Ci for different values of GR and Br.Table 1
Optimal values ofCi  for the case M=2 and different values of GR and Br.
Br
GR
C
1
C
2
C
3
0
400
-
0.96436
-
0.00133
0.00005
8/100
100
-
1.11947
-
0.01693
-
0.00181
0
0
-
0.96436
-
0.00133
0.00005
8/500
500
-
1.19429
-
0.07113
-
0.01860
1
0
-
0.96439
-
0.00133
0.00005
8/100
-
100
-
0.80026
-
0.02050
0.00140
0
±
100
-
0.95976
-
0.00116
0.00003
0
±
500
-
0.95976
-
0.00116
0.00003
8/500
-
500
-
0.79213
-
0.02725
0.00383In Figure1 we compare our approximate four-term solutions (34) against the exact solutions (3) for the special case Br=0 and M=2 for several values of GR. The comparison of the special case GR=0 is shown in Figure 2 for M=2 and several values of Br. It is observed that our four-term OHAM solutions agree very well with the exact solutions. The general case of both Br≠0 and GR≠0 admits no explicit analytical solution. So, in Figures 3 and 4 we plot the four-term approximate OHAM solutions for several values of Br and GR in the case M=2 for both the asymmetric and symmetric heating conditions, respectively. The residual errors corresponding to selected cases of the solutions depicted in Figures 1 and 2 are presented in Figures 5(a) and 5(b), respectively. Finally, the residual errors for a selected case of Figure 3 are shown in Figure 6. Clearly, all the residual error plots suggest that the OHAM approximate solutions are accurate enough.Figure 1
Plots ofu and versus y in the case of asymmetric heating for different values of GR and Br=0, M=2.Figure 2
Plots ofθ versus y in the case of asymmetric heating for different values of Br and GR=0, M=2.Plots of (a)u and (b) θ versus y in the case of asymmetric heating for different values of Br and GR.
(a)
(b)Plots of (a)u and (b) θ versus y in the case of symmetric heating for different values of Br and GR.
(a)
(b)Plots of residual errors for (a)u in the case GR=400, Br=0 and (b) θ in the case GR=0, Br=0.5 of asymmetric heating and M=2.
(a)
(b)Plots of residual errors for (a)u and (b) θ in the case of asymmetric heating for GR=100, Br=8/100 and M=2.
(a)
(b)


## 3.1. Application of OHAM

In this section, we apply OHAM for solving the nonlinear system of two-point BVP (1)-(2). By applying the proposed method, the zeroth-order deformation equation is
(20)(1-p)L[u~(y,p)-u0(y)]=H(p,Ci)[N(u~(y,p))],(1-p)L[θ~(y,p)-θ0(y)]=H(p,Cj)[N(θ~(y,p),u~(y,p))],
subject to the boundary conditions
(21)u~(-14,p)=0,u~(14,p)=0,u~′′(-14,p)=-48+RTGR2,u~′′(14,p)=-48-RTGR2,θ~(-14,y)=-RT2,θ~(14,y)=RT2.
Using the framework of OHAM the mth-order(22)L1[u~-χmum-1]=R1,m(Ci,u→m-1),
where u→m-1={u0,u1,…,um-1},
(23)R1,m=∑i=0m-1‍Ci+1um-1-i(4)-M2GRBr∑i=0m-1‍um-1-i×∑j=0i‍Cj+1uj-i-GRBr∑i=0m-1‍um-1-i′∑j=0i‍Cj+1uj-1′,L2[θ~m-χmθm-1]=R2,m(Ci,θ→m-1),
where
(24)R2,m=∑i=0m-1‍Ci+1θm-1-i′′+M2Br∑i=0m-1‍um-1-i×∑j=0i‍Ci+1uj-i+Br∑i=0m-1‍um-1-i′∑j=0i‍Ci+1uj-1′.
Now the zeroth-order problem is
(25)u0(4)(y)=0,θ0′′(y)=0,
subject to the boundary conditions
(26)u0(-14)=0,u0(14)=0,u0′′(-14)=-48+RTGR2,u0′′(14)=-48-RTGR2,θ0(-14)=-RT2,θ0(14)=RT2.
The solutions are
(27)u0(y)=148(72-1152y2+yRTGR-16y3RTGR),θ0(y)=2yRT.
Now the first-order problem is
(28)u1(4)(y,C1)=-C1M2BrGRu02(y)-C1BrGR(u0′(y))2-C1M2u0′′(y)+C1u0(4)(y),θ1′′(y,C1)=C1M2Bru02(y)+C1Br(u0′(y))2+C1θ0′′(y),
subject to the boundary conditions
(29)u1(-14)=0,u1(14)=0,u1′′(-14)=0,u1′′(14)=0,θ1(-14)=0,θ1(14)=0.
The second-order problem is
(30)u2(4)(y,C1,C2)=-C2M2u0′′(y)+C2u0(4)(y)+u1(4)(y,C1)-C2M2BrGRu02(y)-2C1M2BrGRu0(y)u1(y,C1)-C2BrGR(u0′(y))2+C1u1(4)(y,C1)-C1M2u1′′(y,C1)-2C1BrGRu0′(y)u1′(y,C1),θ2′′(y,C1,C2)=C2M2Bru02(y)+2C1M2Bru0(y)u1(y,C1)+C2Br(u0′(y))2+2c1Bru0′(y)u1′(y,C1)+C2θ0′′(y)+(1+C1)θ1′′(y,C1),
subject to the boundary conditions
(31)u2(-14)=0,u2(14)=0,u2′′(-14)=0,u2′′(14)=0,θ2(-14)=0,θ2(14)=0.
The third-order problem  is
(32)u3(4)(y,C1,C2,C3)=-M2BrC3GRu02(y)-2M2BrC2GRu0(y)u1(y,C1)-C1M2BrGRu12(y,C1)+C3u0(4)(y)+C2u1(4)(y,C1)-C3BrGR(u0′(y))2-2C2BrGRu0′(y)-u1′(y,C1)-C1BrGR(u1′(y,C1))2-C1M2u2′′(y,C1,C2)-C3M2u0′′(y)-C2M2u1′′(y,C1)+C1u2(4)(y,C1,C2)-2C1M2BrGRu0(y)u2(y,C1,C2)+u2(4)(y,C1,C2)-2C1BrGRu0′(y)u2′(y,C1,C2),θ3′′(y,C1,C2,C3)=C3M2Bru02(y)+2C2M2Bru0(y)u1(y,C1)+C1M2Bru12(y,C1)+2C1M2Bru0(y)u2(y,C1,C2)+C3Br(u0′(y))2+2C2Bru0′(y)u1′(y,C1)+C1Br(u1′(y,C1))2+2C1Bru0′(y)u2′(y,C1,C2)+C3θ0′′(y)+C2θ1′′(y,C1)+θ2′′(y,C1,C2)+C1θ2′′(y,C1,C2),
subject to the boundary conditions
(33)u3(-14)=0,u3(14)=0,u3′′(-14)=0,u3′′(14)=0,θ3(-14)=0,θ3(14)=0.Using the solution of (25)–(32) we obtain the following four-term approximate solutions for u and θ by OHAM takingp=1:
(34)u~(y,C1,C2,C3)=u0(y)+u1(y,C1)+u2(y,C1,C2)+u3(y,C1,C2,C3),θ~(y,C1,C2,C3)=θ0(y)+θ1(y,C1)+θ2(y,C1,C2)+θ3(y,C1,C2,C3).
The explicit expressions for the individual terms of the approximate solutions are not given here for brevity. Taking the residual errors
(35)Ru~(y,C1,C2,C3)=u~(4)(y,C1,C2,C3)-M2u~′′(y,C1,C2,C3)-M2GRBru~(y,C1,C2,C3)-GRBr(u~′(y,C1,C2,C3))2,Rθ~(y,C1,C2,C3)=θ~′′(y,C1,C2,C3)+M2Br(u~(y,C1,C2,C3))2-Br(u~′(y,C1,C2,C3))2,
the optimal values of Ci’s can be obtained. Table 1 shows some optimal values of Ci for different values of GR and Br.Table 1
Optimal values ofCi  for the case M=2 and different values of GR and Br.
Br
GR
C
1
C
2
C
3
0
400
-
0.96436
-
0.00133
0.00005
8/100
100
-
1.11947
-
0.01693
-
0.00181
0
0
-
0.96436
-
0.00133
0.00005
8/500
500
-
1.19429
-
0.07113
-
0.01860
1
0
-
0.96439
-
0.00133
0.00005
8/100
-
100
-
0.80026
-
0.02050
0.00140
0
±
100
-
0.95976
-
0.00116
0.00003
0
±
500
-
0.95976
-
0.00116
0.00003
8/500
-
500
-
0.79213
-
0.02725
0.00383In Figure1 we compare our approximate four-term solutions (34) against the exact solutions (3) for the special case Br=0 and M=2 for several values of GR. The comparison of the special case GR=0 is shown in Figure 2 for M=2 and several values of Br. It is observed that our four-term OHAM solutions agree very well with the exact solutions. The general case of both Br≠0 and GR≠0 admits no explicit analytical solution. So, in Figures 3 and 4 we plot the four-term approximate OHAM solutions for several values of Br and GR in the case M=2 for both the asymmetric and symmetric heating conditions, respectively. The residual errors corresponding to selected cases of the solutions depicted in Figures 1 and 2 are presented in Figures 5(a) and 5(b), respectively. Finally, the residual errors for a selected case of Figure 3 are shown in Figure 6. Clearly, all the residual error plots suggest that the OHAM approximate solutions are accurate enough.Figure 1
Plots ofu and versus y in the case of asymmetric heating for different values of GR and Br=0, M=2.Figure 2
Plots ofθ versus y in the case of asymmetric heating for different values of Br and GR=0, M=2.Plots of (a)u and (b) θ versus y in the case of asymmetric heating for different values of Br and GR.
(a)
(b)Plots of (a)u and (b) θ versus y in the case of symmetric heating for different values of Br and GR.
(a)
(b)Plots of residual errors for (a)u in the case GR=400, Br=0 and (b) θ in the case GR=0, Br=0.5 of asymmetric heating and M=2.
(a)
(b)Plots of residual errors for (a)u and (b) θ in the case of asymmetric heating for GR=100, Br=8/100 and M=2.
(a)
(b)


## 4. Conclusion

In this paper we have extended the applicability of OHAM for the first time to solve a nonlinear system of two-point BVPs that arise in a fluid flow problem. OHAM is relatively simple to apply. It was shown that, with a few terms, the OHAM is capable of giving sufficient accuracy. OHAM can be a promising tool for solving strongly nonlinear systems of equations.



---

*Source: 136043-2013-05-15.xml*