# An Efficient Higher-Order Quasilinearization Method for Solving Nonlinear BVPs

Authors: Eman S. Alaidarous; Malik Zaka Ullah; Fayyaz Ahmad; A.S. Al-Fhaid
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

In this research paper, we present higher-order quasilinearization methods for the boundary value problems as well as
coupled boundary value problems. The construction of higher-order convergent methods depends on a decomposition method which is different from Adomain decomposition method (Motsa and Sibanda, 2013). The reported method is very general and can be extended to desired order of convergence for highly nonlinear differential equations and also computationally superior to proposed iterative method based on Adomain decomposition because our proposed iterative scheme avoids the calculations of Adomain polynomials and achieves the same computational order of convergence as authors
have claimed in Motsa and Sibanda, 2013. In order to check the validity and computational performance, the constructed iterative schemes are also successfully applied to bifurcation problems to calculate the values of critical parameters. The numerical performance is also tested for one-dimension Bratu and Frank-Kamenetzkii equations.


---

## Body

## 1. Introduction

In the very beginning, the Newton-Raphson method constructed for nonlinear single-variable algebraic equation whose efficiency index is optimal according to Kung-Traub conjecture [1], actually Newton-Raphson uses two function evaluations, and its computational order of convergence is two. Many authors [2–7] have made good effort to construct the iterative methods for algebraic equations. Similarly, for the system of algebraic equation, the version of Newton-Raphson gives us second-order convergence. In the literature, the efficiency index is only defined for single-variable algebraic equation on the basis of function evaluations. For multivariable case, if we consider system of n algebraic equations, then we require the computation of Jacobian inverse J-1 of dimension n×n and n function evaluations to perform one newton iteration and it is clearly evident the computational cost of matrix inversion is dominant over other binary operations to complete one newton iteration. The iterative methods which require only once the Jacobian inverse for the whole cycle of iterations are clearly efficient.When we talk about nonlinear boundary value problems, the quasilinear method (QLM) [8–11] is an iterative method which starts from initial guess for a boundary value problem (BVP) which is quadratically convergent. Initially, Bellman and Kalaba [12] proposed QLM and later Mandelzweig and coauthors [8–11] provide the second-order convergence proof for the BVPs. In [13], recently authors proposed higher-order quasilinearization method for single BVP as well as coupled BVPs. The original idea in [13] is to decompose the nonlinear operator as an infinite sum of Adomain [14] polynomials. The reported algorithm is efficient in the case of couple BVPs. The computation of Jacobian inverse is performed at the initial guess, but the calculation of Adomain polynomials is somehow difficult and also increases the computation cost of iterative scheme. To avoid the computational burden of Adomain polynomials, we use a different decomposition method for nonlinear operator which was actually introduced in [15]. Our proposed scheme uses only one calculation of Jacobian inverse J-1 and does not require any calculation of Adomain polynomials, and this fact increases its computational efficiency in comparison with [13]. The sequences of iteration schemes have convergence orders two, three, four, five, and so forth. The numerical stability and efficiency are tested over two problems, namely, one-dimensional Bratu problem [16–23]
(1)y′′(x)+λey(x)=0,y(0)=y(1)=0
and Frank-Kamenetzkii [24] boundary value problem
(2)y′′(x)+kxy′(x)+λey(x)=0,y′(0)=y(1)=0.
The Frank-Kamenetzkii BVP (2) has no solution if λ>2, unique solution if λ=2 and two solutions if λ<2. The closed form solution of (2) is reported in [24–26]. The solutions of (2) in [25] are given as
(3)c1=log[2(4-λ)±42(2-λ)],c2=log[4-λ±22(2-λ)2λ2],y=log[16e1c(2λ+ec1x2)2],y=log[16e1c(1+2λec2x2)2].
The closed form solution for Bratu equation [27] can be written as
(4)y(x)=-2log[cosh((x-1/2)(θ/2))cosh(θ/4)]θ=2λcosh(θ4).
The critical parameter λc for Bratu problem satisfies 4=2λcsinh(θc/4) and if λ<λc, λ=λc, and λ>λc, then there are two solutions, unique solution and no solution for (1). The numerical reported value of critical parameter is λc=3.51383071912516 [18].


## 2. Construction of Iterative Methods

### 2.1. Single Nonlinear Boundary Value Problem

Consider a nonlinear ordinary differential equation(5)L(y)+f(y)=0,
where L is a linear derivative operator; for Bratu and Frank-Kamenetzkii problems linear operators are L[y]=d2y/dx2 and L[y]=d2y/dx2+(k/x)(dy/dx), and f is any nonlinear function of y. Let γ be an initial guess (satisfying the boundary conditions) for the solution y of (5). By expanding f around γ, we obtain
(6)f(y-γ+γ)=f(γ)+(y-γ)fy(γ)+g(y),g(y)=f(y)-f(γ)-(y-γ)fy(γ),hhhhhhhhwherefy(γ)=∂f(y)∂y|y=γ,M(y)=y-f(y)fy(γ)-γ+f(γ)fy(γ),hhhh1whereM(y)=-g(y)fy(γ),(7)M(y)+c=y-f(y)fy(γ),wherec=γ-f(γ)fy(γ),(8)L(y)+f(γ)+(y-γ)fy(γ)+g(y)=0,(9)L(y)fy(γ)+f(γ)fy(γ)+(y-γ)=M(y),(10)L(y)fy(γ)+y=c+M(y).
Suppose, that we can decompose the solution y into infinite series sum
(11)y=∑i=0∞yi.
Further, we obtain the decomposition of nonlinear operator M(y) as follows:
(12)M(y)=M(∑i=0∞yi)=M(y0)+(M(y0+y1)-M(y0))+(M(y0+y1+y2)-M(y0+y1))+⋯(13)M(y0+y1+⋯+yn+1)=M(y0)+(M(y0+y1)-M(y0))+(M(y0+y1+y2)-M(y0+y1))+⋯+(M(y0+y1+⋯+yn+1)-M(y0+y1+⋯+yn)).
By substituting (11) in (10), we get
(14)L(y0+y1+y2+⋯)fy(γ)+(y0+y1+y2+⋯)=c+M(∑i=0∞yi).
By using (13), we obtain
(15)L(y0+y1+y2+⋯)fy(γ)+(y0+y1+y2+⋯)=c+M(y0)+(M(y0+y1)-M(y0))+⋯.
By comparing left and right sides in (15), we have
(16)L(y0)fy(γ)+y0=c,(17)L(y1)fy(γ)+y1=M(y0),L(y2)fy(γ)+y2=M(y0+y1)-M(y0)⋮(18)L(yn+1)fy(γ)+yn+1=M(y0+y1+⋯+yn)-M(y0+y1+⋯+yn-1).If we approximate solution(19)y≈Yn+1=y0+y1+⋯+yn+1.
By adding (16) to (18), we get
(20)L(y0+y1+⋯+yn+1)fy(γ)+y0+y1+⋯+yn+1=c+M(y0+y1+⋯+yn),L(y0+y1+⋯+yn+1)fy(γ)+y0+y1+⋯+yn+1=y0+y1+⋯+yn-f(y0+y1+⋯+yn)fy(γ),(21)L(Yn+1)+fy(γ)Yn+1=fy(γ)Yn-f(Yn).
Equations (16) and (21) give
(22)L(Y0)+Y0=fy(γ)γ-f(γ),whereY0=y0,L(Yn+1)+fy(γ)Yn+1=fy(γ)Yn-f(Yn),n=0,1,….After renaming the variables, we obtain the following iterative schemes.Scheme
S
0. Consider
(23)L(y1)+fy(y0)y1=fy(y0)y0-f(y0).
Note that S0 corresponds to the QLM scheme which is quadratically convergent.Scheme
S
1
. Consider
(24)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1).Scheme
S
2
. Consider
(25)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1),L(y3)+fy(y0)y3=fy(y0)y2-f(y2).Scheme
S
n
. Consider
(26)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1),L(y3)+fy(y0)y3=fy(y0)y2-f(y2),⋮L(yn+1)+fy(y0)yn+1=fy(y0)yn-f(yn).
By calculating computational order of convergence, we show that the order of convergence of Sn is n+2.Consider the following nonlinear coupled boundary value problem:(27)L1(z(1))+f1(z)=0,L2(z(2))+f2(z)=0,L3(z(3))+f3(z)=0,⋮Ln(z(n))+fn(z)=0,
where z=[z(1),z(2),z(3),…,z(n)]T and L1,L2,L3,…,Ln are linear derivative operators. Equation (27) can be rewritten as
(28)L(z)+f(z)=0,
where L=[L1,L2,L3,…,Ln]T and f(z)=[f1(z),f2(z),…,fn(z)]T. Let γ=[γ1,γ2,…,γn] be an initial guess, which satisfies boundary conditions, for problem (28). Taylor’s expansion of f around γ is
(29)h1f(z)=f(γ)+J(z-γ)+g(z),whereJisaJacobianJ=∂f∂z|z=γ,g(z)=f(z)-f(γ)-J(z-γ),-J-1g(z)=-J-1f(z)+J-1f(γ)+z-γ,(30)M(z)+c=z-J-1f(z),wherec=γ-J-1f(γ),M(z)+c=z-J-1f(z),wherecM(z)=-J-1g(z).
Equation (28) can be written as
(31)J-1L(z)+J-1f(z)=0,(30)implies(32)J-1L(z)+z=M(z)+c,
where z can be decomposed into infinite series sum
(33)z=∑i=0∞zi.
Substituting (33) in (32), we get
(34)J-1L(∑i=0∞zi)+∑i=0∞zi=M(∑i=0∞zi)+c.
By using the same decomposition for nonlinear operator M(z) for multivariable case which is given in (12), we obtain
(35)J-1L(z0)+z0=c,J-1L(z1)+z1=M(z0),J-1L(z2)+z2=M(z0+z1)-M(z0),J-1L(z3)+z3=M(z0+z1+z2)-M(z0+z1),⋮J-1L(zn+1)+zn+1=M(z0+z1+⋯+zn)-M(z0+z1+⋯+zn-1).If we approximate the solution by(36)z≈Zn+1=z0+z1+⋯+zn+1,
we get by adding (35)
(37)J-1L(Zn+1)+Zn+1=c+M(Zn).
From (32), we have
(38)J-1L(Zn+1)+Zn+1=Zn-J-1f(Zn),L(Zn+1)+JZn+1=JZn-f(Zn).
After renaming the variables, we get the following iterative schemes.Scheme
S
M
0
. Consider
(39)L(z1)+Jz1=Jz0-f(z0),whereJ=∂f∂z|z=z0.
Note that SM0 is the QLM scheme which is quadratically convergent.Scheme
S
M
1
. Consider
(40)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1).Scheme
S
M
2
. Consider
(41)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1),L(z3)+Jz3=Jz2-f(z2).Scheme
S
M
n
. Consider
(42)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1),L(z3)+Jz3=Jz2-f(z2),⋮L(zn+1)+Jzn+1=Jzn-f(zn).


### 2.2. Coupled Boundary Value with Many Variables

Consider the following nonlinear coupled boundary value problem with many variables:(43)L1(z(1))+f1(w)=0,L2(z(2))+f2(w)=0,L3(z(3))+f3(w)=0,⋮Ln(z(n))+fn(w)=0,
where w=[z(1),z(2),z(3),…,z(n),s(1),s(2),…,s(m)]T. The compact form of (43) is
(44)L(z)+f(w)=0,
where z=[z(1),z(2),z(3),…,z(n)]T, s=[s(1),s(2),…,s(m)]T, L=[L1(z(1)),L2(z(2)),…,Ln(z(n))]T, f(w)=[f1(w),f2(w),…,fn(w)]T. Let δ=[γT,ηT]T be an initial guess for (44) which satisfies the boundary conditions where γ=[γ1,γ2,…,γn]T and η=[η1,η2,…,ηm]T. The expansion of f around δ is
(45)f(w)=f(δ)+J(w-δ)+g(w),f(w)=f(δ)+J1(z-γ)+J2(s-η)+g(w),hhhhhhhhhwhereJ1=∂f∂z|z=γ,J2=∂f∂s|s=η,J1-1f(w)=J1-1f(δ)+z-γ+J1-1J2(s-η)+J1-1g(w),M(w)=-J1-1f(w)+J1-1f(δ)+z-γ+J1-1J2(s-η),hhhhhhhhhhhhhhhhhhhwhereM(w)=-J1-1g(w),M(w)+γ+J1-1J2η-J1-1f(δ)=-J1-1f(w)+z+J1-1J2s,(46)M(w)+c=z+J1-1J2s-J1-1f(w),wherec=γ+J1-1J2η-J1-1f(δ),(47)f(w)=J1z+J2s-J1(M(w)+c).
From (44), we obtain
(48)L(z)+J1z+J2s=J1c+J1M(w).
We decompose the solution into infinite series sum
(49)w=∑i=0∞wi.
Equation (48) implies
(50)L(∑i=0∞zi)+J1(∑i=0∞zi)+J2(∑i=0∞si)=J1c+J1M(∑i=0∞wi),L(∑i=0∞zi)+J1(∑i=0∞zi)+J2(∑i=0∞si)=J1c+J1(M(w0)+(M(w0+w1)-M(w0))hhhhhhh1+(M(w0+w1+w2)hhhhhhhhhhhh-M(w0+w1))+⋯(M(w0+w1)-M(w0))),(51)L(z0)+J1z0+J2s0=J1c,(52)L(z0)+J1z0+J2s0=J1γ+J2η-f(δ),L(z0)+Jw0=Jδ-f(δ),whereJ=∂f∂w|w=δ,(53)L(z1)+Jw1=J1M(w0),zn′′(x)+kxzn′(x)+J0zn(x)=J0yn(x)-f(yn(x)),L(z2)+Jw2=J1(M(w0+w1)-M(w0)),⋮L(zn+1)+Jwn+1=J1(M(w0+w1+w2+⋯+wn)-M(w0+w1+w2+⋯+wn-1)).
By adding (51) and (53), we get
(54)L(z0+z1+⋯+zn+1)+J(w0+w1+⋯+wn+1)=J1(c+M(w0+w1+⋯+wn)).
We denote
(55)w≈Wn+1=w0+w1+⋯+wn+1,z≈Zn+1=z0+z1+⋯+zn+1.s≈Sn+1=s0+s1+⋯+sn+1.
From (54) and (46),
(56)L(Zn+1)+JWn+1=J1(c+M(Wn)),L(Zn+1)+JWn+1=J1Zn+J2Sn-f(Wn),L(Zn+1)+JWn+1=JWn-f(Wn).
After renaming the variables, we get the following iterative schemes.Scheme
S
W
0
. Consider
(57)L(z1)+Jw1=Jw0-f(w0),whereJ=∂f∂w|w=w0.Scheme
S
W
1
. Consider
(58)L(z1)+Jw1=Jw0-f(w0),L(w2)+Jw2=Jw1-f(w1).Scheme
S
W
2
. Consider
(59)L(z1)+Jw1=Jw0-f(w0),L(z2)+Jw2=Jw1-f(w1),L(z3)+Jw3=Jw2-f(w2).Scheme
S
W
n
. Consider
(60)L(z1)+Jw1=Jw0-f(w0),L(z2)+Jw2=Jw1-f(w1),L(z3)+Jw3=Jw2-f(w2),⋮L(zn+1)+Jwn+1=Jwn-f(wn).


## 2.1. Single Nonlinear Boundary Value Problem

Consider a nonlinear ordinary differential equation(5)L(y)+f(y)=0,
where L is a linear derivative operator; for Bratu and Frank-Kamenetzkii problems linear operators are L[y]=d2y/dx2 and L[y]=d2y/dx2+(k/x)(dy/dx), and f is any nonlinear function of y. Let γ be an initial guess (satisfying the boundary conditions) for the solution y of (5). By expanding f around γ, we obtain
(6)f(y-γ+γ)=f(γ)+(y-γ)fy(γ)+g(y),g(y)=f(y)-f(γ)-(y-γ)fy(γ),hhhhhhhhwherefy(γ)=∂f(y)∂y|y=γ,M(y)=y-f(y)fy(γ)-γ+f(γ)fy(γ),hhhh1whereM(y)=-g(y)fy(γ),(7)M(y)+c=y-f(y)fy(γ),wherec=γ-f(γ)fy(γ),(8)L(y)+f(γ)+(y-γ)fy(γ)+g(y)=0,(9)L(y)fy(γ)+f(γ)fy(γ)+(y-γ)=M(y),(10)L(y)fy(γ)+y=c+M(y).
Suppose, that we can decompose the solution y into infinite series sum
(11)y=∑i=0∞yi.
Further, we obtain the decomposition of nonlinear operator M(y) as follows:
(12)M(y)=M(∑i=0∞yi)=M(y0)+(M(y0+y1)-M(y0))+(M(y0+y1+y2)-M(y0+y1))+⋯(13)M(y0+y1+⋯+yn+1)=M(y0)+(M(y0+y1)-M(y0))+(M(y0+y1+y2)-M(y0+y1))+⋯+(M(y0+y1+⋯+yn+1)-M(y0+y1+⋯+yn)).
By substituting (11) in (10), we get
(14)L(y0+y1+y2+⋯)fy(γ)+(y0+y1+y2+⋯)=c+M(∑i=0∞yi).
By using (13), we obtain
(15)L(y0+y1+y2+⋯)fy(γ)+(y0+y1+y2+⋯)=c+M(y0)+(M(y0+y1)-M(y0))+⋯.
By comparing left and right sides in (15), we have
(16)L(y0)fy(γ)+y0=c,(17)L(y1)fy(γ)+y1=M(y0),L(y2)fy(γ)+y2=M(y0+y1)-M(y0)⋮(18)L(yn+1)fy(γ)+yn+1=M(y0+y1+⋯+yn)-M(y0+y1+⋯+yn-1).If we approximate solution(19)y≈Yn+1=y0+y1+⋯+yn+1.
By adding (16) to (18), we get
(20)L(y0+y1+⋯+yn+1)fy(γ)+y0+y1+⋯+yn+1=c+M(y0+y1+⋯+yn),L(y0+y1+⋯+yn+1)fy(γ)+y0+y1+⋯+yn+1=y0+y1+⋯+yn-f(y0+y1+⋯+yn)fy(γ),(21)L(Yn+1)+fy(γ)Yn+1=fy(γ)Yn-f(Yn).
Equations (16) and (21) give
(22)L(Y0)+Y0=fy(γ)γ-f(γ),whereY0=y0,L(Yn+1)+fy(γ)Yn+1=fy(γ)Yn-f(Yn),n=0,1,….After renaming the variables, we obtain the following iterative schemes.Scheme
S
0. Consider
(23)L(y1)+fy(y0)y1=fy(y0)y0-f(y0).
Note that S0 corresponds to the QLM scheme which is quadratically convergent.Scheme
S
1
. Consider
(24)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1).Scheme
S
2
. Consider
(25)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1),L(y3)+fy(y0)y3=fy(y0)y2-f(y2).Scheme
S
n
. Consider
(26)L(y1)+fy(y0)y1=fy(y0)y0-f(y0),L(y2)+fy(y0)y2=fy(y0)y1-f(y1),L(y3)+fy(y0)y3=fy(y0)y2-f(y2),⋮L(yn+1)+fy(y0)yn+1=fy(y0)yn-f(yn).
By calculating computational order of convergence, we show that the order of convergence of Sn is n+2.Consider the following nonlinear coupled boundary value problem:(27)L1(z(1))+f1(z)=0,L2(z(2))+f2(z)=0,L3(z(3))+f3(z)=0,⋮Ln(z(n))+fn(z)=0,
where z=[z(1),z(2),z(3),…,z(n)]T and L1,L2,L3,…,Ln are linear derivative operators. Equation (27) can be rewritten as
(28)L(z)+f(z)=0,
where L=[L1,L2,L3,…,Ln]T and f(z)=[f1(z),f2(z),…,fn(z)]T. Let γ=[γ1,γ2,…,γn] be an initial guess, which satisfies boundary conditions, for problem (28). Taylor’s expansion of f around γ is
(29)h1f(z)=f(γ)+J(z-γ)+g(z),whereJisaJacobianJ=∂f∂z|z=γ,g(z)=f(z)-f(γ)-J(z-γ),-J-1g(z)=-J-1f(z)+J-1f(γ)+z-γ,(30)M(z)+c=z-J-1f(z),wherec=γ-J-1f(γ),M(z)+c=z-J-1f(z),wherecM(z)=-J-1g(z).
Equation (28) can be written as
(31)J-1L(z)+J-1f(z)=0,(30)implies(32)J-1L(z)+z=M(z)+c,
where z can be decomposed into infinite series sum
(33)z=∑i=0∞zi.
Substituting (33) in (32), we get
(34)J-1L(∑i=0∞zi)+∑i=0∞zi=M(∑i=0∞zi)+c.
By using the same decomposition for nonlinear operator M(z) for multivariable case which is given in (12), we obtain
(35)J-1L(z0)+z0=c,J-1L(z1)+z1=M(z0),J-1L(z2)+z2=M(z0+z1)-M(z0),J-1L(z3)+z3=M(z0+z1+z2)-M(z0+z1),⋮J-1L(zn+1)+zn+1=M(z0+z1+⋯+zn)-M(z0+z1+⋯+zn-1).If we approximate the solution by(36)z≈Zn+1=z0+z1+⋯+zn+1,
we get by adding (35)
(37)J-1L(Zn+1)+Zn+1=c+M(Zn).
From (32), we have
(38)J-1L(Zn+1)+Zn+1=Zn-J-1f(Zn),L(Zn+1)+JZn+1=JZn-f(Zn).
After renaming the variables, we get the following iterative schemes.Scheme
S
M
0
. Consider
(39)L(z1)+Jz1=Jz0-f(z0),whereJ=∂f∂z|z=z0.
Note that SM0 is the QLM scheme which is quadratically convergent.Scheme
S
M
1
. Consider
(40)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1).Scheme
S
M
2
. Consider
(41)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1),L(z3)+Jz3=Jz2-f(z2).Scheme
S
M
n
. Consider
(42)L(z1)+Jz1=Jz0-f(z0),L(z2)+Jz2=Jz1-f(z1),L(z3)+Jz3=Jz2-f(z2),⋮L(zn+1)+Jzn+1=Jzn-f(zn).


## 2.2. Coupled Boundary Value with Many Variables

Consider the following nonlinear coupled boundary value problem with many variables:(43)L1(z(1))+f1(w)=0,L2(z(2))+f2(w)=0,L3(z(3))+f3(w)=0,⋮Ln(z(n))+fn(w)=0,
where w=[z(1),z(2),z(3),…,z(n),s(1),s(2),…,s(m)]T. The compact form of (43) is
(44)L(z)+f(w)=0,
where z=[z(1),z(2),z(3),…,z(n)]T, s=[s(1),s(2),…,s(m)]T, L=[L1(z(1)),L2(z(2)),…,Ln(z(n))]T, f(w)=[f1(w),f2(w),…,fn(w)]T. Let δ=[γT,ηT]T be an initial guess for (44) which satisfies the boundary conditions where γ=[γ1,γ2,…,γn]T and η=[η1,η2,…,ηm]T. The expansion of f around δ is
(45)f(w)=f(δ)+J(w-δ)+g(w),f(w)=f(δ)+J1(z-γ)+J2(s-η)+g(w),hhhhhhhhhwhereJ1=∂f∂z|z=γ,J2=∂f∂s|s=η,J1-1f(w)=J1-1f(δ)+z-γ+J1-1J2(s-η)+J1-1g(w),M(w)=-J1-1f(w)+J1-1f(δ)+z-γ+J1-1J2(s-η),hhhhhhhhhhhhhhhhhhhwhereM(w)=-J1-1g(w),M(w)+γ+J1-1J2η-J1-1f(δ)=-J1-1f(w)+z+J1-1J2s,(46)M(w)+c=z+J1-1J2s-J1-1f(w),wherec=γ+J1-1J2η-J1-1f(δ),(47)f(w)=J1z+J2s-J1(M(w)+c).
From (44), we obtain
(48)L(z)+J1z+J2s=J1c+J1M(w).
We decompose the solution into infinite series sum
(49)w=∑i=0∞wi.
Equation (48) implies
(50)L(∑i=0∞zi)+J1(∑i=0∞zi)+J2(∑i=0∞si)=J1c+J1M(∑i=0∞wi),L(∑i=0∞zi)+J1(∑i=0∞zi)+J2(∑i=0∞si)=J1c+J1(M(w0)+(M(w0+w1)-M(w0))hhhhhhh1+(M(w0+w1+w2)hhhhhhhhhhhh-M(w0+w1))+⋯(M(w0+w1)-M(w0))),(51)L(z0)+J1z0+J2s0=J1c,(52)L(z0)+J1z0+J2s0=J1γ+J2η-f(δ),L(z0)+Jw0=Jδ-f(δ),whereJ=∂f∂w|w=δ,(53)L(z1)+Jw1=J1M(w0),zn′′(x)+kxzn′(x)+J0zn(x)=J0yn(x)-f(yn(x)),L(z2)+Jw2=J1(M(w0+w1)-M(w0)),⋮L(zn+1)+Jwn+1=J1(M(w0+w1+w2+⋯+wn)-M(w0+w1+w2+⋯+wn-1)).
By adding (51) and (53), we get
(54)L(z0+z1+⋯+zn+1)+J(w0+w1+⋯+wn+1)=J1(c+M(w0+w1+⋯+wn)).
We denote
(55)w≈Wn+1=w0+w1+⋯+wn+1,z≈Zn+1=z0+z1+⋯+zn+1.s≈Sn+1=s0+s1+⋯+sn+1.
From (54) and (46),
(56)L(Zn+1)+JWn+1=J1(c+M(Wn)),L(Zn+1)+JWn+1=J1Zn+J2Sn-f(Wn),L(Zn+1)+JWn+1=JWn-f(Wn).
After renaming the variables, we get the following iterative schemes.Scheme
S
W
0
. Consider
(57)L(z1)+Jw1=Jw0-f(w0),whereJ=∂f∂w|w=w0.Scheme
S
W
1
. Consider
(58)L(z1)+Jw1=Jw0-f(w0),L(w2)+Jw2=Jw1-f(w1).Scheme
S
W
2
. Consider
(59)L(z1)+Jw1=Jw0-f(w0),L(z2)+Jw2=Jw1-f(w1),L(z3)+Jw3=Jw2-f(w2).Scheme
S
W
n
. Consider
(60)L(z1)+Jw1=Jw0-f(w0),L(z2)+Jw2=Jw1-f(w1),L(z3)+Jw3=Jw2-f(w2),⋮L(zn+1)+Jwn+1=Jwn-f(wn).


## 3. Numerical Results and Rate of Convergence

In all numerical experimentation, we use Chebyshev pseudospectral methods (for more details, see [13]). In order to show the rate of convergence, we require the definition of computational order of convergence, The computational order of convergence can be approximated by [28]
(61)COC≈ln∥(yn+1-α)(yn-α)-1∥∞ln∥(yn-α)(yn-1-α)-1∥∞,
where yn-1, yn, and yn+1 are successive iterations closer to the solution α(x) of boundary value problem and ∥vq∥∞=max{|vq(x1)|,|vq(x2)|,…,|vq(xN)|} such that {x1,x2,…,xN} defines the partition of domain of the BVP. For BVPs with infinite domain, for instance, [a,∞), (-∞,a], or (-∞,∞), one could replace infinity by a suitable large number to make the domain compact. The iterative scheme S0 for Bratu problem with boundary conditions is
(62)yn+1′′(x)+λeyn(x)yn+1(x)=λeyn(x)yn(x)-λeyn(x),yn+1(0)=yn+1(1)=0,
and y0(x)=x(1-x) is an initial guess for (62). Similarly, we can obtain the scheme S1 as follows:
(63)zn′′(x)+λeyn(x)zn(x)=λeyn(x)yn(x)-λeyn(x),yn+1′′(x)+λeyn(x)yn+1(x)=λeyn(x)zn(x)-λezn(x).
The S2 and S3 schemes are
(64)zn′′(x)+λeyn(x)zn(x)=λeyn(x)yn(x)-λeyn(x),wn′′(x)+λeyn(x)wn(x)=λeyn(x)zn(x)-λezn(x),yn+1′′(x)+λeyn(x)yn+1(x)=λeyn(x)wn(x)-λewn(x),(65)zn′′(x)+λeyn(x)zn(x)=λeyn(x)yn(x)-λeyn(x),wn′′(x)+λeyn(x)wn(x)=λeyn(x)zn(x)-λezn(x),pn′′(x)+λeyn(x)pn(x)=λeyn(x)wn(x)-λewn(x),yn+1′′(x)+λeyn(x)yn+1(x)=λeyn(x)pn(x)-λepn(x),
respectively. Tables 1 and 2 show the infinity norms of error and rates of convergence for (62), (63), (64), and (65) for λ=1, λ=3 and N runs over 50, 100, 150, and 200. We denote f(y)=λey(x) and fy(y)=∂f(y)/∂y=λey(x) and J0=fy(y0(x)). The construction of S3 scheme for (2) is given below and others are similar. Consider
(66)zn′′(x)+kxzn′(x)+J0zn(x)=J0yn(x)-f(yn(x)),wn′′(x)+kxwn′(x)+J0wn(x)=J0zn(x)-f(zn(x)),pn′′(x)+kxpn′(x)+J0pn(x)=J0wn(x)-f(wn(x)),yn+1′′(x)+kxyn+1′(x)+J0yn+1(x)=J0pn(x)-f(pn(x)),yn+1′(0)=yn+1(1)=0.
The initial guess y0(x)=1-x2 for Frank-Kamenetzkii problem is used to start the proposed schemes. The infinity norms of error and rates of convergence for Frank-Kamenetzkii problem (λ=1,k=1) are depicted in Table 4. In order to plot bifurcation diagram for the one-dimensional Bratu problem and the Frank-Kamenetzkii problem, we rewrite the both problems as follows:
(67)u′′(x)+λeu(x)=0,u(0)=u(1)=0,u′′(x)+kxu′(x)+ηeu(x)=0,u′(0)=u(1)=0.
If λ<λc(=3.51383071912516), λ=λc, and λ>λc, then the Bratu problem has two solutions, unique solution and no solution, respectively, and similarly, for the Frank-Kamenetzkii problem statement is valid if k=1 and ηc=2. We define B1-problem and FK1-problem
(68)u′′(x)+λeu(x)=0,u(0)=0,u′(0)=α,u(1)=0,(69)u′′(x)+kxu′(x)+λeu(x)=0,u′(0)=0,u(0)=α,u(1)=0,
respectively. The iterative forms of (68) and (69) are
(70)un+1′′(x)+λneun(x)un+1(x)+eun(x)λn+1=λneun(x)un(x),u(0)=0,u′(0)=α1,u(1)=0,un+1′′(x)+kxun+1′(x)+ηneun(x)un+1(x)+eun(x)ηn+1=ηneun(x)un(x),u(0)=α2,u′(0)=0,u(1)=0.
The bifurcation diagrams are shown in Figures 1 and 2. The calculation of critical parameters for the B1-problem and FK1-problem is performed by using proposed system of equation in [13] and Tables 2, 5, and 6 show numerical results of different iterative schemes for B1-problem and FK1-problem.Table 1
Convergence order and absolute error for Bratu problem forλ=1 and different values of N.
Iter.
N
50
100
150
200
S
0
1
0.000764899
0.000765714
0.000765861
0.000765913
2
3.22437
e
-
8
3.228
e
-
8
3.22866
e
-
8
3.22889
e
-
8
3
5.67396
e
-
17
5.68041
e
-
17
5.68158
e
-
17
5.68198
e
-
17
4
1.75415
e
-
34
1.75614
e
-
34
1.7565
e
-
34
1.75663
e
-
34
5
7.71203
e
-
49
1.67808
e
-
69
1.67843
e
-
69
1.67855
e
-
69
Rate
2.03003
2.03014
2.03016
2.03017
2.00097
2.00098
2.00098
2.00098
2.00008
2.00008
2.00008
2.00008
S
1
1
0.000010196
0.0000102072
0.0000102092
0.0000102099
2
6.3566
e
-
18
6.36381
e
-
18
6.36511
e
-
18
6.36557
e
-
18
3
7.71203
e
-
49
1.52735
e
-
54
1.52766
e
-
54
1.52777
e
-
54
4
7.71203
e
-
49
1.14254
e
-
95
3.04915
e
-
142
2.11033
e
-
164
Rate
3.02829
3.0284
3.02842
3.02842
2.53302
3.00035
3.00035
3.00035
S
2
1
1.35077
e
-
7
1.35226
e
-
7
1.35253
e
-
7
1.35262
e
-
7
2
2.18092
e
-
31
2.1834
e
-
31
2.18384
e
-
31
2.184
e
-
31
3
7.71203
e
-
49
1.1425376
e
-
95
1.46619
e
-
126
1.46629
e
-
126
Rate
4.02691
4.02701
4.02702
4.02703
S
3
1
1.7892
e
-
9
1.79116
e
-
9
1.79152
e
-
9
1.79164
e
-
9
2
1.79244
e
-
48
1.31437
e
-
48
1.31464
e
-
48
1.31473
e
-
48
3
7.71203
e
-
49
1.14254
e
-
95
3.04915
e
-
142
1.03327
e
-
188
Rate
5.00878
5.02624
5.02626
5.02627Table 2
Convergence order for Bratu problem forλ=3 and different values of N.
Iter.
N
50
100
150
S
0
1
0.0458144
0.045864
0.045873
2
0.000954118
0.000955184
0.000955377
3
4.40809
e
-
7
4.41304
e
-
7
4.41393
e
-
7
4
9.41641
e
-
14
9.42698
e
-
14
9.4289
e
-
14
5
4.29662
e
-
27
4.30146
e
-
27
4.30234
e
-
27
Rate
1.80854
1.80867
1.80869
1.98368
1.98369
1.98369
1.9999
1.9999
1.9999
S
1
1
0.0104901
0.0105016
0.0105037
2
5.27186
e
-
7
5.27777
e
-
7
5.27884
e
-
7
3
6.88292
e
-
20
6.89065
e
-
20
6.89205
e
-
20
4
1.80799
e
-
32
1.53308
e
-
58
1.503335
e
-
58
Rate
2.73823
2.73837
2.73839
2.99715
2.99716
2.99716
S
2
1
0.00253291
0.0025357
0.00253621
2
1.86441
e
-
11
1.86651
e
-
11
1.86688
e
-
11
3
1.80799
e
-
32
5.50708
e
-
44
5.5082
e
-
44
Rate
3.71869
3.71882
3.71885
S
3
1
0.000618859
0.000619541
0.000619664
2
4.0232
e
-
17
4.02771
e
-
17
4.02853
e
-
17
3
1.80799
e
-
32
4.68
e
-
63
4.65698
e
-
83
Rate
4.71115
4.71128
4.7113Figure 1
Bifurcation diagram for Bratu problem (7).Figure 2
Bifurcation diagram for Frank-Kamenetzkii problem (8).The authors of [13] are pioneer to talk about higher-order iterative quasilinearization method (QLM). Tables 1, 2, and 4 confirm the orders of convergence of their respective iterative schemes for the different values of parameters under a different range of grid points for Chebyshev pseudospectral method. We make all the calculation for Tables 1, 2, and 4 in Mathematica (MinPrecision = 200). The scheme-0 in [13] and S0 in this paper are the same because both are QLM and the calculated results should be same but unfortunately this is not the case. For the QLM, the infinity norm of error for Bratu problem is (Scheme-0, N=50, 6.66e−49), (Scheme-2, N=100, 6.34e−49), (Scheme-2, N=150, 6.34e−49), and (Scheme-3, N=200, 6.34e−49) in Table 1 [13] and in this article, in Table 1 (S0, N=50, 7.7e−49), (S1, N=100, 1.67e−69), (S2, N=150, 1.67e−69), (S3, N=200, 1.67e−69) and in other cases, our results show better reduction in error as compared to [13] for Table 1 as well as for Table 2 for Bratu problem. For the case of Frank-Kamenetzkii problem, Table 4 results are comparable with the reported results in Table 4 [13] and especially for QLM case almost results are the same. The Table 3 shows the computation of critical parameter λc for Bratu problem, and again the results produced in this article and in [13] for QLM are not the same which in fact should be the same. The results presented by other authors are surprisingly superior for QLM and for other schemes. Notice that we use Matlab to compute critical parameters. Tables 5 and 6 of [13] show better performance and our results are comparable with them. It is noticed that in [13], the constructed matrix systems (33) and (34) for the Bratu and the Frank-Kamenetzkii problems have not properly implemented for all boundary conditions. It is also noted that in some of the cases, if we increase the grid points by keeping the same scheme, there is an improvement in the accuracy of the calculated results, but in some cases this is not valid rule. It is also very clear from Tables 1, 2, and 4, that we ensure the convergence order which was the claim.Table 3
Iterative convergence for different schemes for the critical value ofλc for Bratu problem, N=200.
Iter.
S
0
S
1
S
2
Criticalλc
0
2.94303552937154
2.94303552937154
2.94303552937154
5
3.63644753861909
3.51376131743441
3.95335426360815
10
3.51383579457353
3.51383064747221
3.51383072002551
15
3.5138308350082
3.5138307190047
3.51383071919088
Absolute errors
0
0.570795189753462
0.570795189753462
0.570795189753462
5
0.122616819494092
6.94016905877781
e
-
005
0.439523544483147
10
5.07544853212138
e
-
006
7.16527859268012
e
-
008
9.00510332968452
e
-
010
15
1.15883229678815
e
-
007
1.20208731857474
e
-
010
6.58797461028371
e
-
011Table 4
Convergence order for Frank-Kamenetzkii problem fork=1, λ=1 and different values of N.
Iter.
N
50
100
150
S
0
1
0.114557
0.114557
0.114557
2
0.00125407
0.00125407
0.00125407
3
1.58064
e
-
7
1.58064
e
-
7
1.58064
e
-
7
4
2.48121
e
-
15
2.48121
e
-
15
2.48121
e
-
15
5
6.09322
e
-
31
6.09322
e
-
31
6.09322
e
-
31
Rate
2.52799
2.52799
2.52799
1.98883
1.98883
1.98883
2.00133
2.00133
2.00133
S
1
1
0.0329952
0.0329952
0.0329952
2
7.69119
e
-
7
7.69119
e
-
7
7.69119
e
-
7
3
8.94424
e
-
21
8.94424
e
-
21
8.94424
e
-
21
4
9.38343
e
-
52
1.39988
e
-
62
1.39988
e
-
62
Rate
3.51966
3.51966
3.51966
3.008
3.008
3.008
S
2
1
0.00877215
0.00877215
0.00877215
2
2.38344
e
-
11
2.38344
e
-
11
2.38344
e
-
11
3
1.25899
e
-
45
1.25899
e
-
45
1.25899
e
-
45
Rate
4.52861
4.52861
4.52861
S
3
1
0.00237755
0.00237755
0.00237755
2
6.27393
e
-
17
6.27393
e
-
17
6.27393
e
-
17
3
9.38343
e
-
52
7.50684
e
-
85
7.50684
e
-
85
Rate
5.52316
5.52316
5.52316Table 5
Iterative convergence for different schemes for the critical value ofλc for Frank-Kamenetzkii problem, N=200.
Iter.
S
0
S
1
S
2
Criticalλc
0
1.96202368624769
1.96202368624769
1.96202368624769
1
2.00031964646299
1.99933495330998
1.99987632667614
2
2.00025140688174
2.00000007821613
1.99999999897116
3
2.00000003077048
2.00000000000462
1.99999999999956
4
2.00000001017936
1.99999999999926
1.99999999999949
Absolute errors
0
0.0379763137523081
0.0379763137523081
0.0379763137523081
1
0.000319646462989098
0.000665046690023985
0.000123673323863427
2
0.000251406881741634
7.82161291290606
e
-
008
1.02883590535896
e
-
009
3
3.07704830504463
e
-
008
4.6216364069096
e
-
012
4.42312853010662
e
-
013
4
1.0179361797924
e
-
008
7.41184891239755
e
-
013
5.14255305006373
e
-
013Table 6
Iterative convergence for different schemes for the critical value ofαc for Frank problem, N=200.
Iter.
S
0
S
1
S
2
Criticalαc
0
1.33333333333333
1.33333333333333
1.33333333333333
1
1.37022694167246
1.38416475416079
1.38597531980146
2
1.38643019839731
1.38630914733076
1.38629436351113
3
1.3864007758502
1.38629428714795
1.38629438822645
4
1.38629414907138
1.38629433591242
1.38629439729082
Absolute errors
0
0.0529610277866668
0.0529610277866668
0.0529610277866668
1
0.0160674194475423
0.0021296069592109
0.00031904131854454
2
0.000135837277305928
1.47862107602315
e
-
005
2.39112551980725
e
-
009
3
0.000106414730195015
7.3972053149518
e
-
008
2.71064524159925
e
-
008
4
2.12048618219995
e
-
007
2.52075842421817
e
-
008
3.61708230034452
e
-
008


## 4. Conclusions

An efficient method is presented in this article to get arbitrary higher-order iterative schemes. The construction of iterative schemes is very simple and straightforward. In [13], authors used the idea to decompose the nonlinear operator by using Adomain decomposition method (ADM) which requires the calculation of Adomain polynomials and in a result, the computational cost will be high. In our case, there is no need to calculate any extra polynomial in order to enhance the order of convergence of an iterative method and both methods require only one inversion of Jacobian. Our presented method is computationally efficient because it does not require to construct any Adomain polynomial and hence implementation is also simple to achieve the same order of convergence. We think that the results presented in [13] and in this article should be the same for the case of QLM, but it is not the case and we do not know the implementation of the iterative schemes in [13].



---

*Source: 259371-2013-11-28.xml*