# Three-Step Block Method for Solving Nonlinear Boundary Value Problems

Authors: Phang Pei See; Zanariah Abdul Majid; Mohamed Suleiman
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a three-step block method of Adam’s type to solve nonlinear second-order two-point boundary value problems of Dirichlet type and Neumann type directly. We also extend this method to solve the system of second-order boundary value problems which have the same or different two boundary conditions. The method will be implemented in predictor corrector mode and obtain the approximate solutions at three points simultaneously using variable step size strategy. The proposed block method will be adapted with multiple shooting techniques via the three-step iterative method. The boundary value problem will be solved without reducing to first-order equations. The numerical results are presented to demonstrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Boundary value problems (BVPs) arise in many areas of applied mathematics, for example, application to chemical reactor theory [1] and Bratu-type problem [2]. Recently many methods are available to solve BVPs such as Adomian decomposition method, variational iteration method, homotopy perturbation method, and modified homotopy perturbation method. The Adomian decomposition method and variational iteration method were developed by Singh and Kumar [3] and Lu [4], respectively. The homotopy perturbation method that has been introduced by Saadatmandi et al. [5] to solve second-order BVPs has less iteration compared to Adomian decomposition method. Asadi et al. [6] has extended the modified homotopy perturbation method to solve the nonlinear system of second-order BVPs. The quintic B-spline collocation method has been modified by Lang and Xu [7] to solve second-order BVPs. Srivastava et al. [8] developed a numerical algorithm based on the nonpolynomial quintic spline functions for the solution of second-order BVPs with engineering applications, while Ibraheem and Khalaf [9] have proposed a shooting neural networks algorithm for solving BVPs. Rahman et al. [10] solved numerically second-order BVPs by the Galerkin method. A smart nonstandard finite difference scheme has been proposed by Erdogan and Ozis [11] for solving second-order nonlinear BVPs. Besides that, Liu et al. [12] solve BVPs with Neumann type by using polynomial spline approach. The aim of this research is to propose a three-step block method to solve the BVPs directly using multiple shooting techniques with variable step size strategy.The present paper is organized as follows. In Section2, we present the derivation of the three-step block method. In Section 3, we show the analysis of the method including the order, consistency, and stability. In Section 4 we show the implementation of the multiple shooting techniques. The numerical results and the discussion are presented in Section 5. Finally, a conclusion is given in Section 6.


## 2. Derivation of Three-Step Block Method

The general two-point second-order BVP and the system of second-order BVP subject to two kinds of boundary condition which are Dirichlet type and Neumann type will be solved directly by three-step block method. The proposed method is the extended of block method proposed by Majid et al. [13] which uses the method to solve second-order ordinary differential equations.The two-point second-order boundary value problem is as follows:(1)y′′=f(x,y,y′),a≤x≤b.
Dirichlet boundary condition:
(2)y(a)=α,y(b)=β.
Neumann boundary condition:Type 1:(3)y′(a)=α,y′(b)=β.Type 2:(4)y(a)=α,y′(b)=β.The system of two-point second-order boundary value problem:
(5)y1′′=f1(x,y1,y1′,y2,y2′),y2′′=f2(x,y1,y1′,y2,y2′),a≤x≤b.
Dirichlet-Dirichlet boundary condition:
(6)y1(a)=α1,y1(b)=β1,y2(a)=α2,y2(b)=β2.
Neumann-Dirichlet boundary condition:
(7)y1(a)=α1,y1′(b)=β1,y2(a)=α2,y2(b)=β2.
We have divided the interval [a,b] into a series of blocks with each block containing three points as shown in Figure 1, where h is the step size and r is the ratio of the step size. Three approximate solutions, yn+1, yn+2, and yn+3 are simultaneously computed using the same back values by integrating (1) once and twice over the intervals [xn,xn+1], [xn,xn+2], and [xn,xn+3], respectively. Consider
(8)yn+1′=yn′+∫xnxn+1f(x,y,y′)dx,yn+1=yn+hyn′+∫xnxn+1(xn+1-x)f(x,y,y′)dx,yn+2′=yn′+∫xnxn+2f(x,y,y′)dx,yn+2=yn+2hyn′+∫xnxn+2(xn+2-x)f(x,y,y′)dx,yn+3′=yn′+∫xnxn+3f(x,y,y′)dx,yn+3=yn+3hyn′+∫xnxn+3(xn+3-x)f(x,y,y′)dx.
The method is derived by replacing the function f(x,y,y′) in (8) with Lagrange interpolation polynomial where five interpolating points are involved and will be implemented using the variable step size strategy. The choices of the next step size will be restricted to half, double, or the same as the current step size. When the next step size is doubled, the ratio r is 0.5 and while the step size remains constant, r is 1. If step fails, the current step size is half the previous step size and the ratio r is 2. Then the approximate solutions in the block will be recalculated. The corrector formulae for r=1 are as follows:
(9)yn+1′=yn′+h720×(-19fn-1+346fn+456fn+1-74fn+2+11fn+3)yn+1=yn+hyn′+h21440×(-21fn-1+472fn+330fn+1-72fn+2+11fn+3)yn+2′=yn′+h90(-fn-1+34fn+114fn+1+34fn+2-fn+3)yn+2=yn+2hyn′+h290(-3fn-1+68fn+114fn+1+fn+3)yn+3′=yn′+3h80(-fn-1+14fn+24fn+1+34fn+2+9fn+3)yn+3=yn+3hyn′+3h2160(-3fn-1+64fn+126fn+1+48fn+2+5fn+3).Figure 1
Three-step Adam’s method.


## 3. Analysis of the Method

In this section we will discuss the order, consistency, stability, and convergence of the three-step Adam’s method. The three-step Adam’s method belongs to the class of linear multistep method (LMM):(10)∑j=0kαjyn+j-h∑j=0kβjyn+j′-h2∑j=0kγjyn+j′′=Cp+2hp+2yp+2.


### 3.1. Order of the Method

Definition 1 (Fatunla [14] and Lambert [15]).
The linear multistep method is said to be of orderp if
(11)C0=C1=C2=⋯=Cp+1=0,Cp+2≠0,
where C0=∑j=0kαj,
(12)C1=∑j=0kjαj-∑j=0kβjC1⋮Cp=1p!∑j=0kjpαj-1(p-1)!∑j=0kjp-1βj-1(p-2)!∑j=0kjp-2γj,ghhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhgpip=2,3,…,
where Cp+2 is the error constant.
Rewrite (9) as (10) with k being 4 in matrix form:
(13)[0000000-110000000000-101000000000-10010][ynyn+1yn+2yn+3yn+4yn+5]=h[01-1000010000010-1000200000100-10030000][yn′yn+1′yn+2′yn+3′yn+4′yn+5′]+h2[-19720346720456720-74720117200-21144047214403301440-7214401114400-1903490114903490-1900-39068901149001900-380428072801028027800-9160192160378160144160151600][fnfn+1fn+2fn+3fn+4fn+5].
From Definition 1, we obtain
(14)C0=C1=⋯=C6=0,C7=[-111440-1100800-1126-3160-9560]T≠0.
Therefore, the order of three-step Adam’s method is five; (p=5) with error constant [-11/1440-1/100800-1/126-3/160-9/560]T.


### 3.2. Consistency of the Method

Definition 2 (Lambert [15]).
The linear multistep method is said to be consistent if it has orderp≥1.
Since the order of three-step Adam’s method isp=5>1, therefore, the method is consistent according to Definition 2.


### 3.3. Stability of the Method

Definition 3 (Lambert [15]).
A linear multistep method is zero-stable provided the rootsξj, j=0(1)k of the first characteristics polynomial ρ(ξ) specified as ρ(ξ)=det⁡⁡|∑j=0kA(j)ξ(k-j)|=0 satisfies |ξj|≤1 and for those roots with |ξj|=1 the multiplicity must not exceed two.
Rewrite (9) in matrix form as follows:
(15)[100000010000001000000100000010000001][yn+1′yn+1yn+2′yn+2yn+3′yn+3]=[000010000001000010000001000010000001][yn-2′yn-2yn-1′yn-1yn′yn]+h[0-19720346720456720-74720117200000000-1903490114903490-1900000000-38042807280102802780000000][fn-2fn-1fnfn+1fn+2fn+3]+h2[0000000-21144047214403301440-7214401114400000000-39068901149001900000000-916019216037816014416015160][fnfn+1fn+2fn+3fn+4fn+5].
From (15), the first characteristic polynomial, ρ(ξ)=det⁡⁡|ξA0-A1|=0, where
(16)A0=[100000010000001000000100000010000001],A1=[000010000001000010000001000010000001],ρ(ξ)=det⁡|[ξ000000ξ000000ξ000000ξ000000ξ000000ξ]-[000010000001000010000001000010000001]|hhhhhhhhhhhh0=ξ4(ξ-1)2ξ=0,0,0,0,1,1.
According to Definition 3, three-step Adam’s method is zero-stable.


### 3.4. Convergence of the Method

Definition 4 (Lambert [15]).
The linear multistep method is convergent if and only if it is consistent and zero-stable.Since the consistency and zero-stable of the method have been established, then the three-step Adam’s method is convergent.


## 3.1. Order of the Method

Definition 1 (Fatunla [14] and Lambert [15]).
The linear multistep method is said to be of orderp if
(11)C0=C1=C2=⋯=Cp+1=0,Cp+2≠0,
where C0=∑j=0kαj,
(12)C1=∑j=0kjαj-∑j=0kβjC1⋮Cp=1p!∑j=0kjpαj-1(p-1)!∑j=0kjp-1βj-1(p-2)!∑j=0kjp-2γj,ghhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhgpip=2,3,…,
where Cp+2 is the error constant.
Rewrite (9) as (10) with k being 4 in matrix form:
(13)[0000000-110000000000-101000000000-10010][ynyn+1yn+2yn+3yn+4yn+5]=h[01-1000010000010-1000200000100-10030000][yn′yn+1′yn+2′yn+3′yn+4′yn+5′]+h2[-19720346720456720-74720117200-21144047214403301440-7214401114400-1903490114903490-1900-39068901149001900-380428072801028027800-9160192160378160144160151600][fnfn+1fn+2fn+3fn+4fn+5].
From Definition 1, we obtain
(14)C0=C1=⋯=C6=0,C7=[-111440-1100800-1126-3160-9560]T≠0.
Therefore, the order of three-step Adam’s method is five; (p=5) with error constant [-11/1440-1/100800-1/126-3/160-9/560]T.


## 3.2. Consistency of the Method

Definition 2 (Lambert [15]).
The linear multistep method is said to be consistent if it has orderp≥1.
Since the order of three-step Adam’s method isp=5>1, therefore, the method is consistent according to Definition 2.


## 3.3. Stability of the Method

Definition 3 (Lambert [15]).
A linear multistep method is zero-stable provided the rootsξj, j=0(1)k of the first characteristics polynomial ρ(ξ) specified as ρ(ξ)=det⁡⁡|∑j=0kA(j)ξ(k-j)|=0 satisfies |ξj|≤1 and for those roots with |ξj|=1 the multiplicity must not exceed two.
Rewrite (9) in matrix form as follows:
(15)[100000010000001000000100000010000001][yn+1′yn+1yn+2′yn+2yn+3′yn+3]=[000010000001000010000001000010000001][yn-2′yn-2yn-1′yn-1yn′yn]+h[0-19720346720456720-74720117200000000-1903490114903490-1900000000-38042807280102802780000000][fn-2fn-1fnfn+1fn+2fn+3]+h2[0000000-21144047214403301440-7214401114400000000-39068901149001900000000-916019216037816014416015160][fnfn+1fn+2fn+3fn+4fn+5].
From (15), the first characteristic polynomial, ρ(ξ)=det⁡⁡|ξA0-A1|=0, where
(16)A0=[100000010000001000000100000010000001],A1=[000010000001000010000001000010000001],ρ(ξ)=det⁡|[ξ000000ξ000000ξ000000ξ000000ξ000000ξ]-[000010000001000010000001000010000001]|hhhhhhhhhhhh0=ξ4(ξ-1)2ξ=0,0,0,0,1,1.
According to Definition 3, three-step Adam’s method is zero-stable.


## 3.4. Convergence of the Method

Definition 4 (Lambert [15]).
The linear multistep method is convergent if and only if it is consistent and zero-stable.Since the consistency and zero-stable of the method have been established, then the three-step Adam’s method is convergent.


## 4. Implementation of the Method

The three-step block method of Adam’s type (3SAM) will be implemented for solving the boundary value problems via multiple shooting techniques. The idea for shooting technique is to form the initial condition from the boundary condition with the guessing value. Multiple shooting techniques are indeed a combination of several shooting techniques by dividing the given intervala≤x≤b into jth subinterval. For the Dirichlet boundary condition, the missing initial condition is y′(a). Equations (1) and (2) can be written as
(17)d2yjdx2=fj(x,yj,yj′),
with initial conditions
(18)y1(a)=α,y1′(a)=svy2(x1)=y(x1,sv-1),y2′(x1)=y′(x1,sv-1)⋮yj(xj-1)=y′(xj-1,sv-1),yj′(xj-1)=y′(xj-1,sv-1).
Therefore, we obtain the jth stop conditions
(19)|y(a,sv)-y(x1,sv-1)|≤ε.|y(x2,sv)-y(x1,sv-1)|≤ε.hhhhhhhhhhhhhhhhhh⋮hhhhhhh|y(xj,sv)-β|≤ε.
The iteration is repeated until we reach the stop conditions, the value of sv will be generated by three-step iterative method as follows:
(20)Tv=sv-y(b,sv)-βz(b,sv)Uv=-y(b,Tv)-βz(b,sv)sv+1=Tv-y(b,sv)-βz(b,sv)-y(b,Tv+Uv)-βz(b,sv).
Detail for the three-step iterative method can be refered to in Yun [16].The missing initial condition for the Neumann boundary condition isy(a), therefore the first initial condition is y1(a)=sv, y1′(a)=α. The stop condition and the three-step iterative method implementation depend to the value of y′(b,sv).


## 5. Numerical Result

In this section, four problems are tested to study the accuracy and the efficiency of the developed codes. The results obtained by the proposed method are compared to the existing method. The following notations are used in the tables:3SAM: Three-step Adam's method variable step size via multiple shooting techniques adapted with three-step iterative method;2P1BVS: two-point block method with variable step size proposed by Phang et al. [17];bvp4c: MATLAB solver proposed by Kierzenka and Shampine [18];MLAM: multilevel augmentation method proposed by Chen [19];COLHW: collocation method with Haar wavelets proposed by Siraj-ul-Islam et al. [20];SCM: sinc-collocation method proposed by Mohamed [21];LRBFM: local radial basis function method proposed by Mehdi and Ahmad [22];TOL: tolerance;TS: total number of steps;MAXE: maximum error;TFC: total function call;TIME: execution time in seconds;RELTOL: the tolerance use to measure the error relative use by bvp4c;ABSTOL: the absolute error tolerances use by bvp4c;MP: total mesh point use by bvp4c;—: no data in the references;1e-2: 1×10-2.Problem 1 (nonlinear Dirichlet boundary value problem).
Consider(21)y′′=-y2+sin2⁡(πx)-π2sin(πx),0≤x≤1.Dirichlet boundary condition:y(0)=0, y(1)=0;
exact solution:y=sin(πx);
source: Chen [19].Problem 2 (nonlinear Neumann boundary value problem).
Consider(22)y′′=2y3,0≤x≤1.Neumann boundary condition:(23)y′(0)=-1,y′(1)=-14.
Exact solution:y=1/(1+x);
source: Siraj-ul-Islam et al. [20].Problem 3 (nonlinear system of boundary value problem).
Consider(24)y1′′=-xy1′-cos⁡(πx)y2′+sin⁡(x)+(x2-x+2)cos⁡⁡(x)+(1-2x)cos⁡⁡(πx)y2′′=-xy′-x(y2′)2-2+xsin⁡(x)+(x2-x)cos⁡⁡(x)+x(1-2x)2.Two Dirichlet boundary condition:(25)y1(0)=0,y1(1)=0,y2(0)=0,y2(1)=0.
Exact solution:y1=(x-1)sin(x), y2=x-x2;
source: Mehdi and Ahmad [22].Problem 4 (nonlinear system of boundary value problem).
The equations governing the free convective boundary-layer flow above a heated impermeable horizontal surface are(26)f′′+mh+(m-23)ηh′=0h′′+(m+13)fh′-mf′h=0.Mixed boundary condition:(27)f(0)=0,f′⟶0asη⟶∞h(0)=1,h⟶0asη⟶∞.Source: Merkin and Zhang [23].Problems1–3 are solved by 3SAM with tolerances 1e-2, 1e-4, 1e-6, 1e-8, and 1e-10. In Problem 1, we solved the boundary value problem subject to the Dirichlet type boundary conditions by 3SAM and compared our result with 2P1BVS, MLAM, and bvp4c. In Problem 2, 3SAM solved the boundary value problem subject to the Neumann type boundary conditions and compared our result to 2P1BVS, COLHW, and bvp4c. In Problem 3, 3SAM solved the system of boundary value problem subject to two boundary conditions and compared our result to SCM and LRBFM. 3SAM, 2P1BVS, and bvp4c are implemented by variable step size strategy and controlled by the tolerances while MLAM and COLHW used constant step size.Tables1–3 showed the comparison of the numerical result for solving Problems 1–3. We have observed that the accuracy for 3SAM is better as the total number of steps was increased for all problems tested. Firstly, we will discuss the block method (3SAM and 2P1BVS) implemented in variable step size strategy. The total number of steps and the total function call taken by 3SAM is less than 2P1BVS in all problems tested; this is expected because 3SAM can obtain the solution at three points simultaneously while 2P1BVS obtains two points simultaneously per step. We also noticed that the maximum error for 2P1BVS is comparable to or better than 3SAM but the accuracy for 3SAM is still within the tolerance. This is because the total step size of 3SAM is less than 2P1BVS therefore the step size used by 3SAM is larger than 2P1BVS. For example, in Table 1 the maximum error for 3SAM is 3.40e-4 with 4 steps and the maximum error for 2P1BVS is 7.66e-5 with 11 steps.Table 1
Comparison of the numerical result for solving Problem1.
3SAM
TOL
1e −2
1e −4
1e −6
1e −8
1e −10
TS
4
6
11
22
45
MAXE
3.40e −4
5.26e −6
8.17e −8
2.36e −10
6.06e −12
TFC
127
184
319
616
1237
Time
0.000139
0.000151
0.000216
0.000307
0.000470
2P1BVS
TOL
1e −2
1e −4
1e −6
1e −8
1e −10
TS
11
16
26
32
57
MAXE
7.66e −5
2.47e −6
3.55e −8
1.16e −8
1.07e −10
TFC
265
377
593
737
1321
Time
0.000141
0.000162
0.000227
0.000275
0.000444
bpv4c
ABSTOL
1e −2
1e −4
1e −6
1e −8
1e −10
RELTOL
1.00
1e −2
1e −4
1e −6
1e −8
MP
15
15
19
69
233
MAXE
6.70e −4
3.35e −5
3.14e −6
1.27e −8
1.25e −10
TFC
157
187
374
1568
5043
Time
0.0313
0.0469
0.0313
0.0938
0.1563
MLAM
h
1/15
1/31
1/63
1/127
1/255
TS
15
31
63
127
255
MAXE
2.46e −4
3.94e −5
4.94e −6
6.17e −7
5.00e −8
Time
0.109
0.235
0.328
0.687
1.953Table 2
Comparison of the numerical result for solving Problem2.
3SAM
TOL
1e −2
1e −4
1e −6
1e −8
1e −10
TS
4
6
10
18
37
MAXE
8.55e −5
8.37e −6
9.68e −7
1.05e −7
1.09e −8
TFC
41
61
98
173
341
Time
0.000082
0.000101
0.000131
0.000229
0.000388
2P1BVS
TOL
1e −2
1e −4
1e −6
1e −8
1e −10
TS
12
16
23
36
59
MAXE
1.72e −5
1.33e −5
2.55e −7
6.20e −9
1.72e −10
TFC
65
89
129
205
341
Time
0.000090
0.000109
0.000141
0.000218
0.000315
bpv4c
ABSTOL
1e −2
1e −4
1e −6
1e −8
1e −10
RELTOL
1.00
1e −2
1e −4
1e −6
1e −8
MP
15
15
29
43
160
MAXE
3.60e −3
4.16e −4
3.17e −7
6.14e −9
2.50e −11
TFC
228
254
491
689
2745
Time
0.0156
0.0313
0.0313
0.0469
0.1406
COLHW
h
1/16
1/32
1/64
1/128
1/256
TS
16
32
64
128
256
MAXE
5.95e −4
1.55e −4
3.98e −5
1.01e −5
2.53e −6Table 3
Comparison of the numerical result for solving Problem3.
3SAM
TOL
1e − 2
1e − 4
1e − 6
1e − 8
1e − 10
TS
4
5
7
13
24
MAXE
3.79e − 04
7.14e − 06
3.57e − 08
8.38e − 10
2.01e − 11
TFC
114
153
231
366
672
Time
0.000203
0.000231
0.000286
0.000396
0.000412
SCM
TS
10
20
30
40
—
MAXE
4.06e − 1
1.70e − 6
8.14e − 8
6.40e − 9
—
LRBFM
TS
21
41
61
—
—
MAXE
2.10e − 4
2.37e − 4
2.01e − 5
—
—Next we discuss the comparison between 3SAM and bvp4c. The bvp4c is a MATLAB solver which uses the collocation formula and a mesh of points to divide the interval of integration into subintervals. If the solution does not satisfy the tolerance, the solver adapts the mesh and repeats the process. In Problems1-2, we choose the initial mesh point as 15. The total function call taken by 3SAM is less than bvp4c in all problems tested. We also noticed that the maximum error for 3SAM is comparable to or better than bvp4c when the tolerance is larger. For the tolerance getting smaller, the accuracy for bvp4c is comparable to or better than 3SAM but the accuracy for 3SAM is still within the tolerance. For example, in Table 2 the maximum error for 3SAM is 1.09e-8 with 341 total functions call and the maximum error for bvp4c is 2.50e-11 with 2745 total functions call. The accuracy of 3SAM is still within the tolerance and compared to the bvp4c and 3SAM is cheaper in terms of total function call and execution time.Finally we discuss the comparison between 3SAM with the method implemented using constant step size. The total number of steps taken by 3SAM is less than MLAM, COLHW, SCM, and LRBFM. This is expected because the 3SAM is using the variable step size strategy while MLAM, COLHW, SCM, and LRBFM are using the constant step size. We noticed that as the tolerance was getting smaller, 3SAM has obtained better accuracy compared to MLAM, COLHW, SCM, and LRBFM. For example, in Table3 the maximum error for 3SAM is 3.57e-8 with 7 steps, SCM has obtained the maximum error 8.14e-8 with 30 steps and LRBFM has obtained the maximum error 2.01e-5 with 61 steps. We also noticed that 3SAM has superiority in terms of execution time compared to MLAM.In Problem4, 3SAM solved the system of boundary value problem subject to the Neumann-Dirichlet type boundary conditions which is a free convective boundary-layer flow in a porous medium above a heated horizontal impermeable surface or below a cooled horizontal impermeable surface where wall temperature is a power function of distance from the origin. Figures 2-3 display the approximate solution of the value h and f′ for selected values m in solving Problem 4.Figure 2
Approximate solution of the valueh against η for selected values of m.Figure 3
Approximate solution of the valuef′ against η for selected values of m.


## 6. Conclusion

In this paper, we have shown that the proposed three-step block method of Adam’s type using variable step size with multiple shooting technique is suitable for solving nonlinear second-order two-point boundary value problems of Dirichlet type, Neumann type, and mixed type of boundary conditions. The numerical results showed that the proposed block method has superiority in terms of accuracy, total function call, total steps, and execution time.



---

*Source: 379829-2014-06-05.xml*