# New Algorithm of Two-Point Block Method for Solving Boundary Value Problem with Dirichlet and Neumann Boundary Conditions

Authors: Pei See Phang; Zanariah Abdul Majid; Fudziah Ismail; Khairil Iskandar Othman; Mohamed Suleiman
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Two-point block method with variable step-size strategy is presented to obtain the solutions for boundary value problems directly. Dirichlet type and Neumann type of boundary conditions are studied in this paper. Multiple shooting techniques adapted with the three-step iterative method are employed for generating the guessing value. Six boundary value problems are solved using the proposed method, and the numerical results are compared to the existing methods. The results suggest a significant improvement in the efficiency of the proposed methods in terms of the number of steps, execution time, and accuracy.


---

## Body

## 1. Introduction

Numerical solutions for boundary value problem have great importance in scientific computation, as they were widely used to model the real-world problems. There are several methods that can be used to solve the two-point boundary value problems such as the Adomian decomposition method modified by Duan and Rach [1] and Ebaid [2] to solve boundary value problem. Based on the Adomian decomposition method, a new analytical and numerical treatment is introduced to investigate linear and nonlinear two-point boundary value problems. Lang and Xu [3] studied a new quintic B-spline collocation method for linear and nonlinear second-order boundary value problems. Islam and Shirin [4] used Bernoulli polynomials to find the numerical solutions of the second-order linear and nonlinear boundary value problems. Besides that, Liu [5] had solved the boundary value problem with Neumann type using polynomial spline approach. The numerical solutions of second-order boundary value problems by collocation method with the Haar wavelets were presented by Siraj-ul-Islam et al. [6]. A new kind of finite difference scheme presented for special second-order nonlinear two-point boundary value problem has been proposed by Erdogan and Ozis [7], while Prentice [8] considered the error control in a finite difference solution of a two-point boundary value problem.In this paper, we are concerned for solving two-point second-order boundary value problems (BVPs) with two types of boundary condition, that is, Dirichlet and Neumann type. Two-point second-order boundary value problem is as follows:(1)y′′=f(x,y,y′),a≤x≤b.
Dirichlet boundary condition is as follows:
(2)y(a)=α,y(b)=β.
Neumann boundary condition is as follows:
(3)y′(a)=α,y′(b)=β.
The BVPs as (1) will be solved using two-point block method with variable step size. The block methods are commonly used to solve the ODEs such as in Majid et al. [9]. The authors have used the two-point block one-step method of Runge-Kutta type to solve the general second-order ODEs with variable step-size strategy.In Majid et al. [10], the authors have solved the BVPs of Dirichlet type only by two-point block method with constant step size. Hence, in this paper, we aim to extend the work in Majid et al. [10] for solving BVPs of Dirichlet and Neumann type using variable step size. This block method has the advantage to solve the second-order differential equation directly and obtain two approximate solutions simultaneously in block. We adapted the multiple shooting techniques to obtain the missing initial value, and the three-step iterative method proposed by Yun [11] was employed to generate the missing guessing value. The numerical computations have been performed using the C language.


## 2. Formulation of the Two-Point Block Method

In Figure1, we have divided the interval [a,b] into a series of blocks with each block containing two points. Both approximate solutions are simultaneously found using the same back values. The approximate solution of yn+1 and yn+2 at the point, xn+1 and xn+2, respectively, with step size h will be computed simultaneously using three back values at the points, xn, xn-1 and xn-2 with step size rh. The value, yn+1 and yn+2 will be obtained by integrate (1) once and twice over the intervals [xn,xn+1] and [xn,xn+2], respectively, as followsFigure 1
Two-point block method.First point is(4)y′(xn+1)-y′(xn)=∫xnxn+1f(x,y,y′)dx,y(xn+1)-y(xn)-hy′(xn)=∫xnxn+1(xn+1-x)f(x,y,y′)dx.Second point is(5)y′(xn+2)-y′(xn)=∫xnxn+2f(x,y,y′)dx,y(xn+2)-y(xn)-2hy′(xn)=∫xnxn+2(xn+2-x)f(x,y,y′)dx.
The method is derived using Lagrange interpolation polynomial and the five interpolating points; that is, (xn-2,fn-2), (xn-1,fn-1), (xn,fn), (xn+1,fn+1), and (xn+2,fn+2) are involved in the corrector formulae of the two-point block method. The function f(x,y,y′) in (4) and (5) will be replaced by Lagrange interpolating polynomial. Let s=(x-xn+2)/h, replacing dx=hds and changing the limit of integration from −2 to −1 for the first point and changing the limit of integration to −2 and 0 for the second point.Evaluate these integrals using MAPLE, the corrector formulae of the two-point block method will be obtained as follows.First point is(6)yn+1′=yn′+h240(1+r)(2+r)(1+2r)r2×[(14+37r+15r2)fn-2+(-28-176r-240r2)fn-1+(14+139r+564r2+1029r3+790r4+200r5)fn+(144r2+672r3+940r4+320r5)fn+1+(-3r2-21r3-50r4-40r5)fn+2],yn+1=yn+hyn′+h2240(2+r)(1+r)(1+2r)r2×[(6+19r+8r2)fn-2+(-12-88r-128r2)fn-1+(6+69r+329r2+644r3+538r4+140r5)fn+(32r2+184r3+324r4+120r5)fn+1+(-r2-8r3-22r4-20r5)fn+2].Second point is(7)yn+2′=yn′+h15(2+r)(1+r)(1+2r)r2×[10r5(-2-r)fn-2+(4-8r)fn-1+(-2-7r+3r2+33r3+35r4+10r5)fn+(48r2+144r3+140r4+40r5)fn+1+(9r2+33r3+35r4+10r5)fn+2],yn+2=yn+2hyn′+h215r(2+r)(1+r)(1+2r)×[6r3(2+r)fn-2+(-8-16r)fn-1+(6+41r+91r2+76r3+20r4)fn+(32r+112r2+128r3+40r4)fn+1+(2r+7r2+6r3)fn+2].
The method is the combination of predictor which is one order less than the corrector. The two-point block method with variable step-size strategy will be implemented for solving the boundary value problems via multiple shooting techniques.


## 3. Analysis of the Method

In this section, stability analysis, stability region, order of the method, and error constant of the two-point block method are discussed.


### 3.1. Stability Analysis

Definition 1.
The method is zero stable provided the rootsRj of the first characteristic polynomial ρ(R) specified as ρ(R)=det[∑i=0kA(i)Rk-i]=0 satisfy |Rj|≤1.
Substituter=1 into (7), and rewrite in matrix from:
(8)[1000010000100001][yn+2yn+2′yn+1yn+1′]=[1000010010000100][ynyn′yn-1yn-1′]+h[0200000001000000][ynyn′yn-1yn-1′]+h[000029901249024904900000-19720346720456720-74720][fn+2fn+1fnfn-1]+h[0000-190000000011720000][fn-2fn-3fn-4fn-5]+h2[590104907890-8900000-17144022014405821440-7614400000][fn-2fn-1fnfn-1]+h2[19000000001114400000000][fn-2fn-3fn-4fn-5].
The first characteristic polynomial of the two-point block method is given as follows:
(9)ρ(R)=det[RA0-A1]=0,whereA0=[1000010000100001],A1=[1000010010000100],ρ(R)=det[R-10000R-100-10R00-10R],0=(R-1)2R2,R=0,0,1,1,|R|≤1.
From Definition 1 and (9) the two-point one block method is zero stable.


### 3.2. Stability Region

The stability polynomial of the two-point block method applied to the test equation is as follows:(10)y′′=f=θy′+λy.The stability polynomials of two-point block method are as follows.Forr=0.5:
(11)(1-147200H1+8475400H12-83450H2+1271350H1H2+21113500H22)t8+(-2+19100H1-43072700H12)t7×(-839225H2-16331350H1H2-7491125H22)t7+(1-124H1+2846327000H12+149450H2)t6+(-127196750H1H2+561113500H22)t6+(4475H1+99250H12-34225H2-7192250H1H2+19150H22)t5+(-8675H12+8675H1H2-83375H22)t4=0,
where H1=h2λ and H2=h2θ.Forr=1.0:
(12)(1-289360H1+4132160H12-524H2+53432H1H2+23910800H22)t8+(-2-31120H1-209216H12)t7+(-1307360H2-161270H1H2-473900H22)t7+(1+3740H1+2536H12-41360H2-283360H1H2+9375400H22)t6+9375400H22t6+(49360H1+891080H12-17360H2-227H1H2+17900H22)t5-12160H12t4+(12160H1H2-110800H22)t4=0,
where H1=h2λ and H2=h2θ.Forr=2.0:
(13)(1-87100H1+6232700H12-103450H2+8475400H1H2+1374500H22)t8+(-2-223960H1-517675H12)t7+(-5121714400H2-51160H1H2-102029216000H22)t7+(1+25792400H1+1500128800H12-14637200H2)t6+(-4000186400H1H2+48139432000H22)t6+(1334800H1+61943200H12-16114400H2-917200H1H2)t5+16154000H22t5+(-186400H12+186400H1H2-1432000H22)t4=0,
where H1=h2λ and H2=h2θ.Figure2 shows the regions of absolute stability for the two-point block method when the step-size ratio r=0.5,1.0,and2.0. The stability region is the bounded shaded region and the region is larger as the step-size ratio increases. This is expected since smaller step-size will give larger step size ratio.Figure 2
Stability region of the two-point block method.


### 3.3. Order and Error Constant

Definition 2.
The block method is(14)∑j=0kαjyn+j-h∑j=0kβjyn+j′-h2∑j=0kγjyn+j′′=Cp+1yp+1+O(hp+2),
where
(15)Cp=1p!∑j=0kjpαj-1(p-1)!∑j=0kjp-1βj-1(p-2)!∑j=0kjp-2γj,C0=C1=C2=⋯=Cp+1=0,Cp+2≠0p is the order of the LMM method, and O(hp+2) is the local truncation error.
Substitutingr=1 into (7) and apply the formulae, we obtaining
(16)C0=C1=⋯=C6=0,C7=[1114403710080-1901315]T.
From Definition 2, the order of two-point block method is five with error constant
(17)C7=[1114403710080-1901315]T.


## 3.1. Stability Analysis

Definition 1.
The method is zero stable provided the rootsRj of the first characteristic polynomial ρ(R) specified as ρ(R)=det[∑i=0kA(i)Rk-i]=0 satisfy |Rj|≤1.
Substituter=1 into (7), and rewrite in matrix from:
(8)[1000010000100001][yn+2yn+2′yn+1yn+1′]=[1000010010000100][ynyn′yn-1yn-1′]+h[0200000001000000][ynyn′yn-1yn-1′]+h[000029901249024904900000-19720346720456720-74720][fn+2fn+1fnfn-1]+h[0000-190000000011720000][fn-2fn-3fn-4fn-5]+h2[590104907890-8900000-17144022014405821440-7614400000][fn-2fn-1fnfn-1]+h2[19000000001114400000000][fn-2fn-3fn-4fn-5].
The first characteristic polynomial of the two-point block method is given as follows:
(9)ρ(R)=det[RA0-A1]=0,whereA0=[1000010000100001],A1=[1000010010000100],ρ(R)=det[R-10000R-100-10R00-10R],0=(R-1)2R2,R=0,0,1,1,|R|≤1.
From Definition 1 and (9) the two-point one block method is zero stable.


## 3.2. Stability Region

The stability polynomial of the two-point block method applied to the test equation is as follows:(10)y′′=f=θy′+λy.The stability polynomials of two-point block method are as follows.Forr=0.5:
(11)(1-147200H1+8475400H12-83450H2+1271350H1H2+21113500H22)t8+(-2+19100H1-43072700H12)t7×(-839225H2-16331350H1H2-7491125H22)t7+(1-124H1+2846327000H12+149450H2)t6+(-127196750H1H2+561113500H22)t6+(4475H1+99250H12-34225H2-7192250H1H2+19150H22)t5+(-8675H12+8675H1H2-83375H22)t4=0,
where H1=h2λ and H2=h2θ.Forr=1.0:
(12)(1-289360H1+4132160H12-524H2+53432H1H2+23910800H22)t8+(-2-31120H1-209216H12)t7+(-1307360H2-161270H1H2-473900H22)t7+(1+3740H1+2536H12-41360H2-283360H1H2+9375400H22)t6+9375400H22t6+(49360H1+891080H12-17360H2-227H1H2+17900H22)t5-12160H12t4+(12160H1H2-110800H22)t4=0,
where H1=h2λ and H2=h2θ.Forr=2.0:
(13)(1-87100H1+6232700H12-103450H2+8475400H1H2+1374500H22)t8+(-2-223960H1-517675H12)t7+(-5121714400H2-51160H1H2-102029216000H22)t7+(1+25792400H1+1500128800H12-14637200H2)t6+(-4000186400H1H2+48139432000H22)t6+(1334800H1+61943200H12-16114400H2-917200H1H2)t5+16154000H22t5+(-186400H12+186400H1H2-1432000H22)t4=0,
where H1=h2λ and H2=h2θ.Figure2 shows the regions of absolute stability for the two-point block method when the step-size ratio r=0.5,1.0,and2.0. The stability region is the bounded shaded region and the region is larger as the step-size ratio increases. This is expected since smaller step-size will give larger step size ratio.Figure 2
Stability region of the two-point block method.


## 3.3. Order and Error Constant

Definition 2.
The block method is(14)∑j=0kαjyn+j-h∑j=0kβjyn+j′-h2∑j=0kγjyn+j′′=Cp+1yp+1+O(hp+2),
where
(15)Cp=1p!∑j=0kjpαj-1(p-1)!∑j=0kjp-1βj-1(p-2)!∑j=0kjp-2γj,C0=C1=C2=⋯=Cp+1=0,Cp+2≠0p is the order of the LMM method, and O(hp+2) is the local truncation error.
Substitutingr=1 into (7) and apply the formulae, we obtaining
(16)C0=C1=⋯=C6=0,C7=[1114403710080-1901315]T.
From Definition 2, the order of two-point block method is five with error constant
(17)C7=[1114403710080-1901315]T.


## 4. Implementation of the Method

### 4.1. Multiple Shooting Technique

The initial conditions being imposed at the same point in the independent variablex, but the boundary conditions are imposed at different values of the independent variable. The idea of 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 interval a≤x≤b into jth subinterval.


#### 4.1.1. Dirichlet-Type Boundary Condition

The missing initial condition isy′(a). Equation (1) can be written as
(18)d2yjdx2=fj(x,yj,yj′),
with conditions
(19)1:y1(a)=α,y1′(a)=sv2:y2(x1)=y(x1,sv-1),y2′(x1)=y′(x1,sv-1)⋮j:yj(xj-1)=y(xj-1,sv-1),yj′(xj-1)=y′(xj-1,sv-1).
Therefore, we obtain the jth stopping conditions as follows:
(20)|y(x1,sv)-y(x1,sv-1)|≤ε|y(x2,sv)-y(x2,sv-1)|≤ε⋮|y(xj,sv)-β|≤ε.
The iteration is repeated until we reached the stopping conditions, and the value of sv will be generated by the three-step iterative method as follows:
(21)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),
where z(b,sv) is the solution of
(22)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=0,z′(a)=1.


#### 4.1.2. Neumann-Type Boundary Condition

The missing initial condition isy(a). Equation (1) can be written as
(23)d2yjdx2=fj(x,yj,yj′),
with conditions given as in (19), but condition 1 will be replaced by
(24)y1(a)=sv,y1′(a)=α.
The stopping conditions are
(25)|y′(x1,sv)-y′(x1,sv-1)|≤ε|y′(x2,sv)-y′(x2,sv-1)|≤ε⋮|y′(xj,sv)-β|≤ε.
The three-step iterative method is
(26)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),
where the value of z′(b,sv) can be obtained from solving the z′′ as follows:
(27)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=1,z′(a)=0.


### 4.2. Variable Step-Size Strategy

The choices of the next step size will be restricted to half, doubled or the same as the current step size. The adjustment is based on the local truncation error (LTE). If the local truncation error is less or equal to the tolerance (TOL), the choice for the next step will be double or remain the same. In the code developed, when the next step size is doubled, the ratior is 0.5, while the step size remains constant, r is 1. If this condition fails, the current step size will be halved from the previous step size; that is, r is 2, and the approximate solutions in the block will be recalculate.Case 1.
LTE≤ TOL (successful step).
Substitutingr=1 in (6) and (7) will produce the following corrector formulae:
(28)yn+1′=yn′+h720×(11fn-2-74fn-1+456fn+346fn+1-19fn+2),yn+1=yn+hyn′+h21440×(11fn-2-76fn-1+582fn+220fn+1-17fn+2),yn+2′=yn′+h90×(-fn-2+4fn-1+24fn+124fn+1+29fn+2),yn+2=yn+2hyn′+h290×(fn-2-8fn-1+78fn+104fn+1+5fn+2).
Substituting r=0.5 in (6) and (7) will produce the following corrector formulae:
(29)yn+1′=yn′+h1800×(145fn-2-704fn-1+1635fn+755fn+1-31fn+2),yn+1=yn+hyn′+h21800×(70fn-2-352fn-1+975fn+220fn+1-13fn+2),yn+2′=yn′+h225×(-20fn-2+64fn-1+15fn+320fn+1+71fn+2),yn+2=yn+2hyn′+2h2225×(5fn-2-32fn-1+120fn+125fn+1+7fn+2).Case 2.
LTE> TOL (failure step).
Substitutingr=2 in (6) and (7) will produce the following corrector formulae:
(30)yn+1′=yn′+h14400×(37fn-2-335fn-1+7456fn+7808fn+1-565fn+237fn-2),yn+1=yn+hyn′+h214400×(19fn-2-175fn-1+4965fn+2656fn+1-265fn+2),yn+2′=yn′+h900×(295fn-2+1216fn-1+285fn+5fn+1-fn+2),yn+2=yn+2hyn′+h2450×(fn-2-10fn-1+345fn+544fn+1+20fn+2).


### 4.3.  Algorithm  2P1BVS

Step 1.
SetTOL and s0.Step 2.
Forn=1,2,3. setxn=a+nh, and evaluate the approximate values yn and zn with direct Adams-Bashforth method.Compute functionsfn and zn′′.Step 3.
Whilexn<b, do Step 4.Step 4.
Fori=1,2. Setxn+i=xn+ih.Evaluate the approximate valuesyn+i and zn+i with two-point block method.Compute functionsfn+i and zn+i′′.Step 5.
If fulfill stopping condition, go to Step8.Step 6.
Generate the new guessing valuessv by three-step iterative method. Set v=v+1.Step 7.
Ifv≥10, set b=(a+b)/2, and go to Step 2.Step 8.
Complete.This algorithm was developed in C language.


## 4.1. Multiple Shooting Technique

The initial conditions being imposed at the same point in the independent variablex, but the boundary conditions are imposed at different values of the independent variable. The idea of 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 interval a≤x≤b into jth subinterval.


### 4.1.1. Dirichlet-Type Boundary Condition

The missing initial condition isy′(a). Equation (1) can be written as
(18)d2yjdx2=fj(x,yj,yj′),
with conditions
(19)1:y1(a)=α,y1′(a)=sv2:y2(x1)=y(x1,sv-1),y2′(x1)=y′(x1,sv-1)⋮j:yj(xj-1)=y(xj-1,sv-1),yj′(xj-1)=y′(xj-1,sv-1).
Therefore, we obtain the jth stopping conditions as follows:
(20)|y(x1,sv)-y(x1,sv-1)|≤ε|y(x2,sv)-y(x2,sv-1)|≤ε⋮|y(xj,sv)-β|≤ε.
The iteration is repeated until we reached the stopping conditions, and the value of sv will be generated by the three-step iterative method as follows:
(21)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),
where z(b,sv) is the solution of
(22)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=0,z′(a)=1.


### 4.1.2. Neumann-Type Boundary Condition

The missing initial condition isy(a). Equation (1) can be written as
(23)d2yjdx2=fj(x,yj,yj′),
with conditions given as in (19), but condition 1 will be replaced by
(24)y1(a)=sv,y1′(a)=α.
The stopping conditions are
(25)|y′(x1,sv)-y′(x1,sv-1)|≤ε|y′(x2,sv)-y′(x2,sv-1)|≤ε⋮|y′(xj,sv)-β|≤ε.
The three-step iterative method is
(26)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),
where the value of z′(b,sv) can be obtained from solving the z′′ as follows:
(27)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=1,z′(a)=0.


## 4.1.1. Dirichlet-Type Boundary Condition

The missing initial condition isy′(a). Equation (1) can be written as
(18)d2yjdx2=fj(x,yj,yj′),
with conditions
(19)1:y1(a)=α,y1′(a)=sv2:y2(x1)=y(x1,sv-1),y2′(x1)=y′(x1,sv-1)⋮j:yj(xj-1)=y(xj-1,sv-1),yj′(xj-1)=y′(xj-1,sv-1).
Therefore, we obtain the jth stopping conditions as follows:
(20)|y(x1,sv)-y(x1,sv-1)|≤ε|y(x2,sv)-y(x2,sv-1)|≤ε⋮|y(xj,sv)-β|≤ε.
The iteration is repeated until we reached the stopping conditions, and the value of sv will be generated by the three-step iterative method as follows:
(21)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),
where z(b,sv) is the solution of
(22)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=0,z′(a)=1.


## 4.1.2. Neumann-Type Boundary Condition

The missing initial condition isy(a). Equation (1) can be written as
(23)d2yjdx2=fj(x,yj,yj′),
with conditions given as in (19), but condition 1 will be replaced by
(24)y1(a)=sv,y1′(a)=α.
The stopping conditions are
(25)|y′(x1,sv)-y′(x1,sv-1)|≤ε|y′(x2,sv)-y′(x2,sv-1)|≤ε⋮|y′(xj,sv)-β|≤ε.
The three-step iterative method is
(26)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),
where the value of z′(b,sv) can be obtained from solving the z′′ as follows:
(27)z′′=∂f(x,y,y′)∂yz+∂f(x,y,y′)∂y′z′,z(a)=1,z′(a)=0.


## 4.2. Variable Step-Size Strategy

The choices of the next step size will be restricted to half, doubled or the same as the current step size. The adjustment is based on the local truncation error (LTE). If the local truncation error is less or equal to the tolerance (TOL), the choice for the next step will be double or remain the same. In the code developed, when the next step size is doubled, the ratior is 0.5, while the step size remains constant, r is 1. If this condition fails, the current step size will be halved from the previous step size; that is, r is 2, and the approximate solutions in the block will be recalculate.Case 1.
LTE≤ TOL (successful step).
Substitutingr=1 in (6) and (7) will produce the following corrector formulae:
(28)yn+1′=yn′+h720×(11fn-2-74fn-1+456fn+346fn+1-19fn+2),yn+1=yn+hyn′+h21440×(11fn-2-76fn-1+582fn+220fn+1-17fn+2),yn+2′=yn′+h90×(-fn-2+4fn-1+24fn+124fn+1+29fn+2),yn+2=yn+2hyn′+h290×(fn-2-8fn-1+78fn+104fn+1+5fn+2).
Substituting r=0.5 in (6) and (7) will produce the following corrector formulae:
(29)yn+1′=yn′+h1800×(145fn-2-704fn-1+1635fn+755fn+1-31fn+2),yn+1=yn+hyn′+h21800×(70fn-2-352fn-1+975fn+220fn+1-13fn+2),yn+2′=yn′+h225×(-20fn-2+64fn-1+15fn+320fn+1+71fn+2),yn+2=yn+2hyn′+2h2225×(5fn-2-32fn-1+120fn+125fn+1+7fn+2).Case 2.
LTE> TOL (failure step).
Substitutingr=2 in (6) and (7) will produce the following corrector formulae:
(30)yn+1′=yn′+h14400×(37fn-2-335fn-1+7456fn+7808fn+1-565fn+237fn-2),yn+1=yn+hyn′+h214400×(19fn-2-175fn-1+4965fn+2656fn+1-265fn+2),yn+2′=yn′+h900×(295fn-2+1216fn-1+285fn+5fn+1-fn+2),yn+2=yn+2hyn′+h2450×(fn-2-10fn-1+345fn+544fn+1+20fn+2).


## 4.3.  Algorithm  2P1BVS

Step 1.
SetTOL and s0.Step 2.
Forn=1,2,3. setxn=a+nh, and evaluate the approximate values yn and zn with direct Adams-Bashforth method.Compute functionsfn and zn′′.Step 3.
Whilexn<b, do Step 4.Step 4.
Fori=1,2. Setxn+i=xn+ih.Evaluate the approximate valuesyn+i and zn+i with two-point block method.Compute functionsfn+i and zn+i′′.Step 5.
If fulfill stopping condition, go to Step8.Step 6.
Generate the new guessing valuessv by three-step iterative method. Set v=v+1.Step 7.
Ifv≥10, set b=(a+b)/2, and go to Step 2.Step 8.
Complete.This algorithm was developed in C language.


## 5. Problem Tested

In this section, we have tested the algorithm 2P1BVS to six numerical examples to illustrate its accuracy and efficiency.Problem 1.
We have(31)y′′=y′-siny-y2+x4-2x3+3+sin(x2-x),0≤x≤1.
Dirichlet boundary condition is as follows:y(0)=0, y(1)=0.
Exact solution is as follows:y=x2-x.Problem 2.
We have(32)y′′=-y2+sin2(πx)-π2sin(πx),0≤x≤1.
Dirichlet boundary condition is as follows:y(0)=0, y(1)=0.
Exact solution is as follows:y=sin(πx).Problem 3.
We have(33)y′′=-(2-4x2)y,0≤x≤1.
Neumann boundary condition is as follows:y′(0)=0, y′(1)=-(-2/exp(1)).
Exact solution is as follows:y=exp(-x2).Problem 4.
We have(34)y′′=2y3,0≤x≤1.
Neumann boundary condition is as follows:y′(0)=-1, y′(1)=-1/4.
Exact solution is as follows:y=1/(1+x).Problem 5.
We have(35)y′′=-y-1,0≤x≤b.
Neumann boundary condition:y′(0)=(1-cos(1))/sin(1), y′(b)=cos(b)+((1-cos(1))/sin(1))sin(b)-1.
Exact solution is as follows:y=cos(x)+((1-cos(1))/sin(1))sin(x)-1.Problem 6 (Troesch’s Problem).
We have(36)y′′=λsinh(λy),0≤x≤b.
Dirichlet boundary condition is as follows:y(0)=0,y(1)=1.
Troesch’s problem comes from the investigation of the confinement of a plasma column under radiation pressure. The closed form solution to this problem has been given in terms of Jacobian elliptic function by Lin et al. [12].


## 6. Numerical Result and Discussion

The following notations are used in the tables:TOL: ToleranceTS: Total number of stepsMAXE: Maximum errorTime (sec): Execution time in secondTFC: Total function calls2P1BVS: Two-point block method with variable step size via multiple shooting technique adapted with three-step iterative methodMLAM: Multilevel augmentation method proposed by Chen [13]COLHW: Collocation method with Haar wavelets proposed by Siraj-ul-Islam [6]FDM: Finite different method proposed by Erdogan and Ozis [7].In Problems1 and 2, we solved the boundary value problem with Dirichlet-type boundary conditions by 2P1BVS and compare our result with MLAM. In Table 1, the 2P1BVS has superiority in terms of accuracy and execution time when compared to MLAM. For example, the maximum error for the 2P1BVS and MLAM with 15 steps is 1.94e-16 and 2.42e-4 respectively. This result is expected since 2P1BVS solves the second-order BVPs directly and obtains two approximate solutions simultaneously. In Table 2, we observed that the 2P1BVS can obtain the maximum error 1.07e-10 with 57 steps, but MLAM can obtain the maximum error 5.00e-8 with 255 steps.Table 1
Comparison of the numerical result for solving Problem1.
2P1BVS
MLAM
TOL
TS
MAXE
Time (sec)
TS
MAXE
Time (sec)
1
.
00
e
-
2
12
5.34
e
-
5
0.000118
15
7.43
e
-
4
1.719
1.00
e
-
4
15
1.94
e
-
16
0.000130
31
2.42
e
-
4
3.750
1.00
e
-
6
19
1.30
e
-
16
0.000162
63
6.06
e
-
5
8.141
1.00
e
-
8
22
1.80
e
-
16
0.000169
127
1.51
e
-
5
16.875
1.00
e
-
10
25
6.71
e
-
17
0.000203
255
3.78
e
-
6
37.328
511
9.45
e
-
7
85.016
1023
2.36
e
-
7
190.12
2047
6.76
e
-
8
435.23Table 2
Comparison of the numerical result for solving Problem2.
2P1BVS
MLAM
TOL
TS
MAXE
Time (sec)
TS
MAXE
Time (sec)
1
.
00
e
-
2
11
7.66
e
-
5
0.000141
15
2.46
e
-
4
0.109
1.00
e
-
4
16
2.47
e
-
6
0.000162
31
3.94
e
-
5
0.235
1.00
e
-
6
26
3.55
e
-
8
0.000227
63
4.94
e
-
6
0.328
1.00
e
-
8
32
1.16
e
-
8
0.000275
127
6.17
e
-
7
0.687
1.00
e
-
10
57
1.07
e
-
10
0.000444
255
5.00
e
-
8
1.953In Problems3 and 4, 2P1BVS solved the boundary value problem with Neumann-type boundary conditions when was compared to COLHW. The maximum errors for 2P1BVS in Tables 3 and 4 are comparable to COLHW for the larger tolerance. For example, in Table 3, the maximum error for 2P1BVS with 12 steps and COLHW with 16 steps is 3.64e-4 and 2.91e-4, respectively. As the tolerance is getting smaller, the maximum error and total number of steps for 2P1BVS are setting better than COLHW. In Table 4, the maximum error for 2P1BVS with 108 steps and COLHW with 512 steps is 3.47e-12 and 6.34e-7, respectively.Table 3
Comparison of the numerical result for solving Problem3.
2P1BVS
COLHW
TOL
TS
MAXE
Time (sec)
TS
MAXE
1.00
e
-
2
12
3.64
e
-
4
0.000092
16
2.91
e
-
4
1.00
e
-
4
16
3.83
e
-
5
0.000111
32
7.48
e
-
5
1.00
e
-
6
23
5.36
e
-
7
0.000149
64
1.90
e
-
5
1.00
e
-
8
32
6.31
e
-
8
0.000193
128
4.77
e
-
6
1.00
e
-
10
54
1.09
e
-
9
0.000306
256
1.20
e
-
6Table 4
Comparison of the numerical result for solving Problem4.
2P1BVS
COLHW
TOL
TS
MAXE
Time (sec)
TS
MAXE
1.00
e
-
2
12
1.72
e
-
5
0.000090
16
5.95
e
-
4
1.00
e
-
4
16
1.33
e
-
5
0.000109
32
1.55
e
-
4
1.00
e
-
6
23
2.55
e
-
7
0.000141
64
3.98
e
-
5
1.00
e
-
8
36
6.20
e
-
9
0.000218
128
1.01
e
-
5
1.00
e
-
10
59
1.72
e
-
10
0.000315
256
2.53
e
-
6Table5 displays the numerical results for Problem 5 with larger interval at b=5 and b=10 for solving Problem 5. The 2P1BVS managed to obtain the accuracy within the given tolerances. We also observed that the 2P1BVS manages to converge rapidly for the larger interval, for example, at TOL = 10-10, the execution time is 0.000551 and 0.001083 seconds when b=5 and b=10, respectively.Table 5
Numerical result for solving Problem5 with different endpoint, b.
TOL
b
=
5
b
=
10
TS
TFC
MAXE
Time (sec)
TS
TFC
MAXE
Time (sec)
1
.
00
e
-
2
14
81
2.40
e
-
3
0.000116
17
111
4.35
e
-
2
0.000139
1.00
e
-
4
22
129
3.49
e
-
5
0.000177
28
177
3.64
e
-
4
0.000210
1.00
e
-
6
31
179
1.54
e
-
6
0.000245
46
273
3.10
e
-
6
0.000322
1.00
e
-
8
55
323
8.71
e
-
8
0.000376
92
553
2.99
e
-
7
0.000577
1.00
e
-
10
87
509
3.38
e
-
9
0.000551
175
1045
5.57
e
-
9
0.001083We consider the cases whenλ=0.5 and λ=1.0 to solve the Troesch problem at TOL = 10-5. The numerical results in Tables 6 and 7 show that the 2P1BVS manage to solve the problem as accurate as or better than FDM. Figures 3, 4, 5, and 6 display the comparison of the maximum errors versus the total steps for the numerical results in Tables 1–7.Table 6
Comparison of the numerical result for solving Troesch’s Problem withλ=0.5.
x
Solution
Error
Exact
FDM
2P1BVS
FDM
2P1BVS
0.1
0.0959443493
0.0959443492
0.0959443493
5.35
e
-
11
0.00
0.2
0.1921287477
0.1921287476
0.1921287477
1.04
e
-
10
0.00
0.3
0.2887944009
0.2887944007
0.2887944009
1.47
e
-
10
0.00
0.4
0.3861848464
0.3861848462
0.3861848464
1.80
e
-
10
9.99
e
-
16
0.5
0.4845471647
0.4845471645
0.4845471647
1.99
e
-
10
3.60
e
-
14
0.6
0.5841332484
0.5841332482
0.5841332484
2.02
e
-
10
5.61
e
-
13
0.7
0.6852011483
0.6852011481
0.6852011483
1.86
e
-
10
5.79
e
-
12
0.8
0.7880165227
0.7880165225
0.7880165227
1.49
e
-
10
4.41
e
-
11
0.9
0.8928542161
0.8928542161
0.8928542158
9.12
e
-
11
2.66
e
-
10Table 7
Comparison of the numerical result for solving Troesch’s Problem withλ=1.0.
x
Solution
Error
Exact
FDM
2P1BVS
FDM
2P1BVS
0.1
0.0846612565
0.0846612556
0.0846612565
9.30
e
-
10
0.00
0.2
0.1701713582
0.1701713565
0.1701713582
1.70
e
-
9
0.00
0.3
0.2573939080
0.2573939059
0.2573939080
2.18
e
-
9
1.90
e
-
13
0.4
0.3472228551
0.3472228528
0.3472228551
2.23
e
-
9
1.48
e
-
11
0.5
0.4405998351
0.4405998333
0.4405998347
1.80
e
-
9
4.46
e
-
10
0.6
0.5385343980
0.5385343971
0.5385343906
9.10
e
-
10
7.37
e
-
09
0.7
0.6421286091
0.6421286094
0.6421285278
3.02
e
-
10
8.13
e
-
08
0.8
0.7526080939
0.7526080954
0.7526080937
1.47
e
-
9
2.43
e
-
10
0.9
0.8713625196
0.8713625215
0.8713625152
1.87
e
-
9
4.42
e
-
09Figure 3
Comparison of maximum error and total step for Problem1.Figure 4
Comparison of maximum error and total step for problem 2.Figure 5
Comparison of maximum error and total step for problem 3.Figure 6
Comparison of maximum error and total step for problem 4.


## 7. Conclusion

In this paper, we have shown that the proposed two-point block method using variable step size (2P1BVS) is suitable for solving directly two-point second-order boundary value problems in Dirichlet- and Neumann-type-boundary conditions.



---

*Source: 917589-2013-04-03.xml*