# A Subdivision Based Iterative Collocation Algorithm for Nonlinear Third-Order Boundary Value Problems

Authors: Syeda Tehmina Ejaz; Ghulam Mustafa
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

We construct an algorithm for the numerical solution of nonlinear third-order boundary value problems. This algorithm is based on eight-point binary subdivision scheme. Proposed algorithm is stable and convergent and gives more accurate results than fourth-degree B-spline algorithm.


---

## Body

## 1. Introduction

Many problems in physics, chemistry, and engineering science are demonstrated mathematically by third-order boundary value problems. These boundary value problems can be found in different areas of applied mathematics and physics as, in the deflection of a curved beam having a constant or varying cross section, a three-layer beam, electromagnetic waves, or gravity driven flows. Third-order boundary value problems were discussed in many papers in recent years. Different types of techniques have been used to study such problems: finite difference method, spline method, Haar wavelets method, Adomian decomposition, modified Adomian decomposition method (MADM) shooting technique via Newton method, He’s homotopy perturbation method, and so forth. Haq et al. [1] studied the Haar wavelets method to solve the third-order boundary value problem. A modified form of Adomian decomposition method is applied to construct the numerical solution for solving third-order boundary value problems proposed by Duan and Rach [2] and Hasan [3]. Abdullah et al. [4] solved the third-order boundary value problem directly using the fourth-order two-point block method adapted with the shooting technique via Newton method. Caglar et al. [5] solved third-order linear and nonlinear boundary value problems by using fourth-degree B-splines. Their algorithm has first-order accuracy. A new algorithm for solving the general nonlinear third-order differential equation is developed by [6] using shifted Jacobi-Gauss collocation spectral method.For finite differences methods, only discrete approximate values of the unknowny
(
x
) can be obtained. We need further data processing techniques to get accurate fitted curve to data. For spline interpolation or approximation methods the unknown function y
(
x
) is assumed to be piecewise polynomial which in turn requires at least piecewise higher order differentiability of the function f
(
x
,
y
,
u
,
v
).To overcome the above disadvantages, Qu and Agarwal [7, 8] introduced the subdivision based algorithm for the solution of two-point second-order boundary value problems. Mustafa and Ejaz [9] solved linear third-order boundary value problems by using subdivision technique. Ejaz et al. [10] solved two-point fourth-order linear boundary value problem by subdivision based method. Higher order linear and nonlinear problems are not solved by subdivision techniques until now. This motivates us to solve nonlinear third-order boundary value problems by subdivision schemes based collocation iterative algorithms.An outline of the paper is as follows. In Section2, some results about existence and uniqueness of the solution of third-order boundary value problem are given. In Section 3, subdivision algorithm, basis function, and their derivatives are briefed. In Section 4, subdivision based iterative algorithm for the solution of nonlinear third-order BVPs using the derivatives of basis functions is formulated. Convergence of the proposed algorithm is also discussed in this section. Error analysis is given in Section 5. Numerical examples illustrating the usefulness of our proposed algorithm are given in Section 6.


## 2. Existence and Uniqueness of the Solution

In this section, we present some results about the existence and uniqueness of the solution of third-order nonlinear boundary value problems. The details of these results can be found in [11].The general third-order nonlinear boundary value problem can be prescribed as(1)
y
′
′
′
=
f
x
,
y
,
y
′
,
y
′
′with boundary conditions defined as(2)
y
0
=
α
1
,
y
′
0
=
α
2
,
y
1
=
α
3
,where α
i
,
i
=
1,2
,
3 are constants.Proposition 1.
If the functionf
x
,
y
,
y
′
,
y
′
′ is continuous and satisfies the following uniform Lipschitz condition:(3)
f
x
,
y
,
y
′
,
y
′
′
-
f
x
,
y
~
,
y
′
~
,
y
′
′
~
≤
M
0
y
-
y
~
+
M
1
y
′
-
y
′
~
+
M
2
y
′
′
-
y
′
′
~
,
∀
x
,
y
,
y
′
,
y
′
′
,
x
,
y
~
,
y
′
~
,
y
′
′
~
∈
a
,
b
×
R
3
,where the constants M
0, M
1, and M
2 satisfy(4)
2
81
M
0
b
-
a
3
+
1
6
M
1
b
-
a
2
+
2
3
M
2
b
-
a
<
1or(5)
3
160
M
0
b
-
a
3
+
17
150
M
1
b
-
a
2
+
2
3
M
2
b
-
a
<
1then the boundary value problem has one and only one solution:(6)
y
′
′
′
=
f
x
,
y
,
y
′
,
y
′
′
,
a
≤
x
≤
b
,
y
a
=
y
1
,
y
′
a
=
y
2
,
y
b
=
y
3
.The existence and uniqueness of the differential equation (1) with boundary conditions at two or three points are presented in [11].Remark 2.
Throughout this paper the functionf
x
,
y
,
y
′
,
y
′
′ satisfies Lipschitz conditions (3) along with condition (4)–(6). So the existence and uniqueness of the solutions of (1) are guaranteed.


## 3. Subdivision Scheme and Basis Function

In this section, we define binary subdivision algorithm and their basis function that are used to construct the approximate solutions of (1).


### 3.1. Interpolating Subdivision Scheme

We consider the following 8-point binary interpolating subdivision scheme introduced in [12, 13](7)
p
2
i
k
+
1
=
p
i
k
,
p
2
i
+
1
k
+
1
=
1225
2048
p
i
k
+
p
i
+
1
k
-
245
2048
p
i
-
1
k
+
p
i
+
2
k
+
49
2048
p
i
-
2
k
+
p
i
+
3
k
-
5
2048
p
i
-
3
k
+
p
i
+
4
k
,where p
i
k and p
i
k
+
1 are points at kth and (
k
+
1
)th iterative level. Scheme (7) is C
3-continuous and support length (
-
7,7
) with 8th-order approximation.


### 3.2. Basis Function and Their Derivatives

The basis function is the limit function resulting from cardinal data, where all vertices of the polygon have value zero except for one. Letg
(
x
)
,
x
∈
R be the fundamental solution of (7) that satisfies the two-scale equation: (8)
g
x
=
g
2
x
+
1
2048
1225
g
2
x
-
1
+
g
2
x
+
1
-
245
g
2
x
-
3
+
g
2
x
+
3
+
49
g
2
x
-
5
+
g
2
x
+
5
-
5
g
2
x
-
7
+
g
2
x
+
7
,
x
∈
R
,
g
x
∈
C
3
,
g
x
=
0
,
i
f
x
∉
-
6,6
,
g
i
=
δ
0
,
i
∈
Z
.Proposition 3.
The fundamental solution is three times continuously differentiable over the interval[
-
6,6
]. Its derivatives at integers are given by (9)
g
′
i
=
2
sign
⁡
i
E
i
T
η
1
,
g
′
′
i
=
2
2
E
i
T
η
2
,
g
′
′
′
i
=
2
3
sign
⁡
i
E
i
T
η
3
, where (10)
sgn
⁡
t
=
-
1
,
t
<
0
,
0
,
t
=
0
,
1
,
t
>
0
, and E
t’s for 0
≤
t
≤
6 are defined below (11)
E
t
=
e
6
t
,
e
5
t
,
e
4
t
,
e
3
t
,
e
2
t
,
e
1
t
,
e
0
t
,
e
-
1
t
,
e
-
2
t
,
e
-
3
t
,
e
-
4
t
,
e
-
5
t
,
e
-
6
t
T
, where (12)
e
i
t
=
1
,
i
=
t
,
0
,
i
≠
t
, and η
j, 1
≤
j
≤
3, are defined in [9].Furthermore, the numeric values of first, second, and third derivatives ofg
(
i
) at i
∈
[
-
6,6
] are(13)
g
′
0
=
0
,
g
′
±
1
=
∓
78592
49553
,
g
′
±
2
=
±
76113
198212
,
g
′
±
3
=
∓
3328
49553
,
g
′
±
4
=
±
2645
594636
,
g
′
±
5
=
±
256
743295
,
g
′
±
6
=
∓
1
594636
,
(14)
g
′
′
0
=
-
342643
41124
,
g
′
′
±
1
=
5704256
1079505
,
g
′
′
±
2
=
-
12053651
8636040
,
g
′
′
±
3
=
325952
1079505
,
g
′
′
±
4
=
-
60871
2878680
,
g
′
′
±
5
=
-
704
215901
,
g
′
′
±
6
=
55
1727208
,
(15)
g
′
′
′
0
=
0
,
g
′
′
′
±
1
=
∓
292352
117495
,
g
′
′
′
±
2
=
±
3047987
1879920
,
g
′
′
′
±
3
=
∓
3312
13055
,
g
′
′
′
±
4
=
∓
1369
234990
,
g
′
′
′
±
5
=
±
16
2611
,
g
′
′
′
±
6
=
∓
5
41776
.The above derivative values are found by using the left eigenvectors of the subdivision process (7). The detailed description about these left eigenvectors η
j
,
1
≤
j
≤
3, and derivatives can be found in [9].


## 3.1. Interpolating Subdivision Scheme

We consider the following 8-point binary interpolating subdivision scheme introduced in [12, 13](7)
p
2
i
k
+
1
=
p
i
k
,
p
2
i
+
1
k
+
1
=
1225
2048
p
i
k
+
p
i
+
1
k
-
245
2048
p
i
-
1
k
+
p
i
+
2
k
+
49
2048
p
i
-
2
k
+
p
i
+
3
k
-
5
2048
p
i
-
3
k
+
p
i
+
4
k
,where p
i
k and p
i
k
+
1 are points at kth and (
k
+
1
)th iterative level. Scheme (7) is C
3-continuous and support length (
-
7,7
) with 8th-order approximation.


## 3.2. Basis Function and Their Derivatives

The basis function is the limit function resulting from cardinal data, where all vertices of the polygon have value zero except for one. Letg
(
x
)
,
x
∈
R be the fundamental solution of (7) that satisfies the two-scale equation: (8)
g
x
=
g
2
x
+
1
2048
1225
g
2
x
-
1
+
g
2
x
+
1
-
245
g
2
x
-
3
+
g
2
x
+
3
+
49
g
2
x
-
5
+
g
2
x
+
5
-
5
g
2
x
-
7
+
g
2
x
+
7
,
x
∈
R
,
g
x
∈
C
3
,
g
x
=
0
,
i
f
x
∉
-
6,6
,
g
i
=
δ
0
,
i
∈
Z
.Proposition 3.
The fundamental solution is three times continuously differentiable over the interval[
-
6,6
]. Its derivatives at integers are given by (9)
g
′
i
=
2
sign
⁡
i
E
i
T
η
1
,
g
′
′
i
=
2
2
E
i
T
η
2
,
g
′
′
′
i
=
2
3
sign
⁡
i
E
i
T
η
3
, where (10)
sgn
⁡
t
=
-
1
,
t
<
0
,
0
,
t
=
0
,
1
,
t
>
0
, and E
t’s for 0
≤
t
≤
6 are defined below (11)
E
t
=
e
6
t
,
e
5
t
,
e
4
t
,
e
3
t
,
e
2
t
,
e
1
t
,
e
0
t
,
e
-
1
t
,
e
-
2
t
,
e
-
3
t
,
e
-
4
t
,
e
-
5
t
,
e
-
6
t
T
, where (12)
e
i
t
=
1
,
i
=
t
,
0
,
i
≠
t
, and η
j, 1
≤
j
≤
3, are defined in [9].Furthermore, the numeric values of first, second, and third derivatives ofg
(
i
) at i
∈
[
-
6,6
] are(13)
g
′
0
=
0
,
g
′
±
1
=
∓
78592
49553
,
g
′
±
2
=
±
76113
198212
,
g
′
±
3
=
∓
3328
49553
,
g
′
±
4
=
±
2645
594636
,
g
′
±
5
=
±
256
743295
,
g
′
±
6
=
∓
1
594636
,
(14)
g
′
′
0
=
-
342643
41124
,
g
′
′
±
1
=
5704256
1079505
,
g
′
′
±
2
=
-
12053651
8636040
,
g
′
′
±
3
=
325952
1079505
,
g
′
′
±
4
=
-
60871
2878680
,
g
′
′
±
5
=
-
704
215901
,
g
′
′
±
6
=
55
1727208
,
(15)
g
′
′
′
0
=
0
,
g
′
′
′
±
1
=
∓
292352
117495
,
g
′
′
′
±
2
=
±
3047987
1879920
,
g
′
′
′
±
3
=
∓
3312
13055
,
g
′
′
′
±
4
=
∓
1369
234990
,
g
′
′
′
±
5
=
±
16
2611
,
g
′
′
′
±
6
=
∓
5
41776
.The above derivative values are found by using the left eigenvectors of the subdivision process (7). The detailed description about these left eigenvectors η
j
,
1
≤
j
≤
3, and derivatives can be found in [9].


## 4. Subdivision Based Iterative Algorithm

In this section, we describe the algorithm for the numerical solution of nonlinear boundary value problem (1).


### 4.1. The Collocation Algorithm

In this subsection, we construct the collocation method based on the interpolating subdivision scheme (7). Let U
(
x
) be the assumed solution of (1):(16)
U
x
=
∑
i
=
-
6
N
+
6
u
i
g
x
-
x
i
h
,
0
⩽
x
⩽
1
,where N
(
⩾
6
)
∈
Z
+, h
=
1
/
N, x
i
=
i
/
N
=
i
h, and {
u
i
} are unknown to be determined for the solution of (1). The collocation algorithm together with the boundary conditions is defined as follows:(17)
U
′
′
′
x
j
=
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
,with the following type of boundary conditions:(18)
U
0
=
α
1
,
U
′
0
=
α
2
,
U
N
=
α
3
,by taking the third derivative of (16) we get(19)
U
′
′
′
x
=
1
h
3
∑
i
=
-
6
N
+
6
u
i
g
′
′
′
x
-
x
i
h
,
0
⩽
x
⩽
1
,using (19) into (17), we get(20)
∑
i
=
-
6
N
+
6
u
i
g
′
′
′
x
j
-
x
i
h
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.This can be written as(21)
∑
i
=
-
6
N
+
6
u
i
g
j
-
i
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.As we know that g
i
′
′
′
=
-
g
-
i
′
′
′, then above system of equations becomes(22)
∑
i
=
-
6
N
+
6
-
1
u
i
g
i
-
j
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.Now we simplify the nonlinear system of (22) in the following theorems.Theorem 4.
The nonlinear system of (22) for j
=
0 becomes(23)
∑
i
=
-
6
6
-
1
u
i
g
i
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.Proof.
By expanding (22) for j
=
0, we get(24)
-
1
u
-
6
g
-
6
′
′
′
+
u
-
5
g
-
5
′
′
′
+
⋯
+
u
N
+
5
g
N
+
5
′
′
′
+
u
N
+
6
g
N
+
6
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.As we know g
′
′
′
(
i
) exists only for the interval for i
∈
[
-
6,6
] and outside the interval it will be zero. Then above equation can be written as (25)
-
1
u
-
6
g
-
6
′
′
′
+
u
-
6
g
-
6
′
′
′
+
⋯
+
u
5
g
5
′
′
′
+
u
6
g
6
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.This implies (23).Theorem 5.
Forj
=
1,2
,
…
,
N, the nonlinear system of (22) becomes(26)
∑
i
=
-
6
+
j
6
+
j
-
1
u
i
g
i
-
j
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
.Proof.
Substitutingj
=
1 in (22), it becomes (27)
-
1
u
-
6
g
-
7
′
′
′
+
u
-
5
g
-
6
′
′
′
+
⋯
+
u
7
g
6
′
′
′
+
⋯
+
u
N
+
6
g
N
+
5
′
′
′
=
h
3
f
x
1
,
U
x
1
,
U
′
x
1
.
g
′
′
′
(
i
) is nonzero only for the interval for i
∈
[
-
6,6
] and outside the interval it will be zero. Then the above equation becomes(28)
-
1
u
-
5
g
-
6
′
′
′
+
u
-
4
g
-
5
′
′
′
+
⋯
+
u
4
g
5
′
′
′
+
u
5
h
6
′
′
′
=
h
3
f
x
1
,
U
x
1
,
U
′
x
1
.For j
=
2, (22) becomes (29)
-
1
u
-
6
g
-
8
′
′
′
+
u
-
5
g
-
7
′
′
′
+
u
-
4
g
-
6
′
′
′
+
⋯
+
u
4
g
2
′
′
′
+
u
5
g
3
′
′
′
+
u
6
g
4
′
′
′
+
u
7
g
5
′
′
′
u
8
g
6
′
′
′
+
⋯
+
u
N
+
5
g
N
+
3
′
′
′
+
u
N
+
6
g
N
+
4
′
′
′
=
h
3
f
x
2
,
U
x
2
,
U
′
x
2
.This implies(30)
-
1
u
-
4
g
-
6
′
′
′
+
u
-
3
g
-
5
′
′
′
+
⋯
+
u
7
g
5
′
′
′
+
u
8
g
6
′
′
′
=
h
4
f
x
2
,
U
x
2
,
U
′
x
2
.Similarly, we can find the expression for j
=
3,4
,
…
,
N; that is,(31)
-
1
u
-
3
g
-
6
′
′
′
+
u
-
4
g
-
5
′
′
′
+
⋯
+
u
8
g
5
′
′
′
+
u
9
g
6
′
′
′
=
h
3
f
x
3
,
U
x
3
,
U
′
x
3
,
-
1
u
-
2
g
-
6
′
′
′
+
u
-
1
g
-
5
′
′
′
+
⋯
+
u
9
g
5
′
′
′
+
u
10
g
6
′
′
′
=
h
3
f
x
4
,
U
x
4
,
U
′
x
4
,
⋮
-
1
u
N
-
6
g
-
6
′
′
′
+
u
N
-
5
g
-
5
′
′
′
+
⋯
+
u
N
+
5
g
5
′
′
′
+
u
N
+
6
g
6
′
′
′
=
h
3
f
x
N
,
U
x
N
,
U
′
x
N
.Hence, by combining (28), (30), and (31) we get (26).


### 4.2. Unstable Nonlinear System

The nonlinear system of (22) is equivalent to the following nonlinear system of N
+
1 equations with (N
+
13) unknowns {
u
i
}:(32)
A
1
U
=
F
u
,where A
1 is banded matrix of order (
N
+
1
)
×
(
N
+
13
), U is the unknown vector of order N
+
13, and F
(
u
) is the vector of order N
+
1 that depends on u. The matrix A and vectors U and F
(
u
) are given explicitly by(33)
A
1
=
-
1
g
p
q
′
′
′
q
-
p
-
6
N
+
1
×
N
+
13
,where p and q represent the row and column, respectively; that is, p
=
1,2
,
3
,
…
,
N
+
1 and q
=
1,2
,
3
,
…
,
N
+
13,(34)
F
u
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
,
…
,
h
3
f
x
N
,
U
x
N
,
U
′
x
N
T
,
(35)
U
=
u
-
6
,
u
-
5
,
z
-
4
,
…
,
u
N
+
4
,
u
N
+
5
,
u
N
+
6
T
,
U
′
x
j
=
∑
i
=
-
6
N
+
6
u
j
g
′
x
j
-
x
i
h
.
g
′
i is defined in (13) and g
(
i
)
=
g
i. System (32) is unstable and we need to make it stable to get unique solution. The detail for the stable system of nonlinear equations is given in next section.


### 4.3. Stable Nonlinear System

For unique solution of nonlinear systems (32), we need twelve more conditions. Three conditions can be attained from given boundary conditions for nonlinear systems of equations and remaining nine conditions are attained by setting some extrapolation method. The details of the given boundary conditions and extrapolation method are given below.


#### 4.3.1. Approximated Boundary Condition

The given boundary conditions are(36)
U
0
=
α
1
,
U
′
0
=
α
2
,
U
N
=
α
3
.We see that first derivative is involved in the given boundary conditions, since approximation order of interpolating scheme (7) is eight, so we approximate derivative boundary conditions at end point with approximation order eight. The approximation of derivative conditions at end point is defined as(37)
U
′
0
=
N
840
-
2283
u
0
+
6720
u
1
-
11760
u
2
+
15680
u
3
-
14700
u
4
+
9408
u
5
-
3920
u
6
+
960
u
7
-
105
u
8
+
O
h
8
.


#### 4.3.2. Imposed Boundary Conditions

Remaining nine conditions for the nonlinear systems (32) to get stable systems for the solution of (1) are obtained by setting the following extrapolation method.We define five conditions at left end points and four conditions at the right end points. Since subdivision scheme reproduces seven-degree polynomials, so we define boundary conditions of order eight for solution of (32). For simplicity only the left end points u
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 are discussed and the values of right end points u
N
+
2
,
u
N
+
3
,
u
N
+
4
,
u
N
+
5 can be treated similarly.The valuesu
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 can be determined by the septic polynomial R
(
x
) interpolating (
x
i
,
u
i
)
,
0
≤
i
≤
5. Precisely, we have (38)
u
-
i
=
R
-
x
i
,
i
=
1,2
,
3,4
,
5
, where (39)
R
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
U
x
i
-
r
. Since by (16), U
(
x
i
)
=
u
i for i
=
1,2
,
3,4
,
5 then by replacing x
i by -
x
i, we have (40)
R
-
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
u
-
i
+
r
. Hence, the following boundary conditions can be employed at the left end:(41)
∑
r
=
0
8
8
r
-
1
r
u
-
i
+
r
=
0
,
i
=
5,4
,
3,2
,
1
.Similarly for the right end, we can define u
i
=
R
(
-
x
i
), i
=
N
+
2, N
+
3, N
+
4, N
+
5, and (42)
R
x
i
=
∑
r
=
0
8
8
r
-
1
r
+
1
u
i
-
r
. So we have the following boundary conditions at the right end:(43)
∑
r
=
0
8
8
r
-
1
r
u
i
-
r
=
0
,
i
=
N
+
2
,
N
+
3
,
N
+
4
,
N
+
5
.Finally, we get the following new system of (
N
+
13
) linear equations with (
N
+
13
) unknowns {
u
i
}, in which N
+
1 equations are obtained from (22), three equations from boundary conditions (36), and nine from boundary conditions (41) and (43).Hence, the stable nonlinear system of equations is defined as(44)
A
U
=
G
u
,where the matrix A is given by(45)
A
=
A
0
T
,
A
1
T
,
A
2
T
T
,and the matrix A
1 is defined in (33). The matrix (
A
0
)
7
×
(
N
+
13
) is constructed by taking first five rows from (41), sixth row from (37), and the last row U
(
0
) from (36). Hence,(46)
A
0
=
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
⋯
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
⋯
0
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
⋯
0
0
0
0
0
0
0
0
-
2283
N
840
6720
N
840
-
11760
N
840
15680
N
840
-
14700
N
840
9408
N
840
-
3920
N
840
960
N
840
-
105
N
840
0
⋯
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
⋯
0
0
. The matrix (
A
2
)
6
×
(
N
+
13
) is constructed by taking first row U
(
N
) from (36) and last five rows are obtained from (43). Hence,(47)
A
2
=
0
0
⋯
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
⋯
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
0
⋯
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
. The column vector U is defined in (35) and G is defined as(48)
G
u
=
0,0
,
0,0
,
0
,
U
′
0
,
U
0
,
F
T
u
,
U
1
,
0,0
,
0,0
T
,where F
(
u
) is defined by (34).


#### 4.3.3. Nonsingularity of a Matrix

The nonsingularity of coefficient matrixA, which is defined in (45), can be checked by different methods.Since the determinant of matrixA is nonzero for N
≤
2000 then the nonlinear system of equations give solution for N
≤
2000.It is also observed that all the eigenvalues are nonzero forN
≤
2000. Hence, by [14] matrix A is nonsingular. For large N
>
2000 the matrix may or may not be nonsingular.The coefficient matrixA is neither symmetric nor diagonally dominant. Though it can be proved that A is nonsingular/invertible, the matrix is almost symmetric except the first and last few rows and columns due to its boundary treatment for large value of N. Therefore, we first consider the symmetric part of it, that is, the square band matrix B of order N
+
1 defined as (49)
-
1
×
B
=
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
g
4
′
′
′
⋯
⋯
0
0
0
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
⋯
⋯
0
0
0
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
⋯
⋯
0
0
0
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
⋯
⋯
0
0
0
g
-
3
′
′
′
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
⋯
⋯
0
0
0
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
⋯
⋯
g
2
′
′
′
g
2
′
′
′
g
3
′
′
′
0
0
0
0
⋯
⋯
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
0
0
0
0
⋯
⋯
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
.Kılıç and Stanica [15] presented a method to find the inverse of banded matrix of order n by using the LU factorization of the banded matrix. As B is a band matrix of order (
N
+
1
) so by LU factorization method its inverse exists.


### 4.4. Iterative Algorithm and Its Convergence

In this section, we propose an iterative algorithm and discuss its convergence.


#### 4.4.1. Iterative Algorithm Based on Basis Function

The iterative algorithm is based on basis function of the subdivision scheme (7) as defined in the following three steps.First Step (initial approximation). Initial approximation is important because numerical solution depends on the initial approximation. We define the process for finding the initial approximation as follows.Let initial approximate solutionU
0 be the solution of the following linear system:(50)
A
U
0
=
G
0
,where (51)
G
0
=
0,0
,
0,0
,
0
,
y
′
a
,
y
a
,
p
0
,
p
1
,
p
2
,
…
,
p
N
,
y
b
,
0,0
,
0,0
T
,
p
i
=
h
3
p
x
i
,
m
i
,
S
,
i
=
0,1
,
2
,
…
,
N
,
a
≤
x
≤
b
,
m
i
=
y
a
+
i
h
y
b
-
y
a
b
-
a
,
S
=
y
b
-
y
a
.
G
0 is the initial linear approximation of the nonlinear vector G
(
u
). By solving the linear system of (50) we get initial approximate solution.Second Step (numerical solution). The numerical solutions U
∗ of the nonlinear system are obtained by using simple iterative scheme:(52)
A
U
m
+
1
=
G
U
m
,
m
=
0,1
,
2,3
,
…
.Third Step (stopping condition). The above iterative processes will terminate when the following condition is satisfied(53)
U
m
-
U
m
-
1
≤
t
o
l
,where tol is supposed value; that is, t
o
l
=
1
0
-
6. The convergence of the above iterative algorithm is guaranteed by the following propositions. The solutions of linear system of (50) and (52) are obtained by Gaussian eliminations method.Proposition 6.
The successive solutions{
U
m
} for the nonlinear system (44) generated by the iterative algorithm (52) linearly converge to the solution U
∗ provided that M
0 and M
1 are Lipschitz constants and step size h is small. That is,(54)
A
-
1
M
0
h
3
+
867307
212370
M
1
h
2
≤
1
.Proof.
LetU
∗ and U
(
m
) be the solutions of the nonlinear system (44) then by definition, for small h, we have (55)
A
U
∗
=
G
U
∗
,
A
U
m
+
1
=
G
U
m
. Let the error vector be defined as e
(
m
)
=
U
m
-
U
∗ at mth iteration which satisfies (56)
A
U
m
+
1
-
A
U
∗
=
G
U
m
-
G
U
∗
,
G
U
m
+
1
-
U
∗
=
G
U
m
-
G
U
∗
,
A
e
m
+
1
=
G
U
m
-
G
U
∗
. For i
=
0,1
,
2
,
…
,
N, (57)
D
3
e
i
m
+
1
=
F
U
m
-
F
U
∗
i
,by using the mean value theorem, which is stated as “if a function f
(
x
,
y
,
z
) is continuously differentiable in an open set of R
3 containing points (
x
1
,
y
1
,
z
1
) and (
x
2
,
y
2
,
z
2
) and the line segment connecting them, then an equation (58)
f
x
2
,
y
2
,
z
2
-
f
x
1
,
y
1
,
z
1
=
f
x
′
r
,
s
,
t
x
2
-
x
1
+
f
y
′
r
,
s
,
t
y
2
-
y
1
+
f
z
′
r
,
s
,
t
z
2
-
z
1 is valid for the interior point (
a
,
b
,
c
) of the segment.” We have (59)
D
3
e
i
m
+
1
=
f
x
i
,
U
i
m
,
U
′
m
-
f
x
i
,
U
i
∗
,
U
′
∗
. The above equation can be written as (by using mean value theorem) (60)
D
3
e
i
m
+
1
=
f
x
∗
x
i
-
x
i
+
f
y
∗
U
i
m
-
U
i
∗
+
f
y
′
∗
U
′
m
-
U
′
∗
, and by using the definition of error vector, we have (61)
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
e
i
′
m
,
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
D
1
e
i
m
, where D
3 and D
1 are the derivative difference operators defined as (62)
D
1
f
i
=
1
2973180
h
-
5
f
i
+
6
-
f
i
-
6
+
1024
f
i
+
5
-
f
i
-
5
+
13225
f
i
+
4
-
f
i
-
4
-
199680
f
i
+
3
-
f
i
-
3
+
1141695
f
i
+
2
-
f
i
-
2
-
4715520
f
i
+
1
-
f
i
-
1
,
D
3
f
i
=
1
1879920
h
3
-
225
f
i
+
6
-
f
i
-
6
+
11520
f
i
+
5
-
f
i
-
5
-
10952
f
i
+
4
-
f
i
-
4
-
476928
f
i
+
3
-
f
i
-
3
+
3047987
f
i
+
2
-
f
i
-
2
-
4677632
f
i
+
1
-
f
i
-
1
. This implies (63)
D
3
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. Since e
i
=
e
N
-
i
=
0
,
i
=
0
,
-
1
,
-
2
,
…
,
-
6, therefore, we have (64)
A
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. This can be written as (65)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By taking norm on both sides, we get (66)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By using the definition of Lipschitz condition, we get(67)
e
i
m
+
1
≤
A
-
1
h
3
M
0
e
i
m
+
h
2
M
1
D
1
e
i
m
.This implies (68)
e
i
m
+
1
e
i
m
≤
A
-
1
h
3
M
0
+
h
2
M
1
D
1
, which is equivalent to (69)
e
i
m
+
1
e
i
m
≈
h
2
M
1
A
-
1
D
1
. The results follow immediately from this inequality and the following fact: (70)
D
1
=
867307
212370
.A simple approximation of condition by omitting the cubic term is (71)
h
≤
212370
867307
M
1
-
1
A
-
1
-
1
1
/
2
.This completes the proof.Proposition 7.
Iff satisfies the Lipschitz condition (3) and Lipschitz constants M
0, M
1, and M
2 and mesh size h are small enough then nonlinear system (32) has a unique solution U
∗. A sufficient condition for the existence of a solution is given by (54).Proof.
From the proof of previous proposition, we observe that if (54) holds, then an inequality similar to that of (67) holds, which implies that the sequence {
U
m
} is contracting and hence converges. The limit U
∗ of (52) also satisfies (44) due to the continuity of the right hand side function f
x
,
y
,
y
′
,
y
′
′.Remark 8.
The numerical complexity for the solution of the linear system (52), where the matrix A is almost band matrix with half band 7, using Gaussian elimination method is about 49
(
N
+
13
) multiplications. The number of complexities depends upon the efficient boundary treatment. If more efficient boundary treatment is constructed then number of complexities will be reduced.


## 4.1. The Collocation Algorithm

In this subsection, we construct the collocation method based on the interpolating subdivision scheme (7). Let U
(
x
) be the assumed solution of (1):(16)
U
x
=
∑
i
=
-
6
N
+
6
u
i
g
x
-
x
i
h
,
0
⩽
x
⩽
1
,where N
(
⩾
6
)
∈
Z
+, h
=
1
/
N, x
i
=
i
/
N
=
i
h, and {
u
i
} are unknown to be determined for the solution of (1). The collocation algorithm together with the boundary conditions is defined as follows:(17)
U
′
′
′
x
j
=
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
,with the following type of boundary conditions:(18)
U
0
=
α
1
,
U
′
0
=
α
2
,
U
N
=
α
3
,by taking the third derivative of (16) we get(19)
U
′
′
′
x
=
1
h
3
∑
i
=
-
6
N
+
6
u
i
g
′
′
′
x
-
x
i
h
,
0
⩽
x
⩽
1
,using (19) into (17), we get(20)
∑
i
=
-
6
N
+
6
u
i
g
′
′
′
x
j
-
x
i
h
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.This can be written as(21)
∑
i
=
-
6
N
+
6
u
i
g
j
-
i
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.As we know that g
i
′
′
′
=
-
g
-
i
′
′
′, then above system of equations becomes(22)
∑
i
=
-
6
N
+
6
-
1
u
i
g
i
-
j
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
,
j
=
0,1
,
2
,
…
,
N
.Now we simplify the nonlinear system of (22) in the following theorems.Theorem 4.
The nonlinear system of (22) for j
=
0 becomes(23)
∑
i
=
-
6
6
-
1
u
i
g
i
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.Proof.
By expanding (22) for j
=
0, we get(24)
-
1
u
-
6
g
-
6
′
′
′
+
u
-
5
g
-
5
′
′
′
+
⋯
+
u
N
+
5
g
N
+
5
′
′
′
+
u
N
+
6
g
N
+
6
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.As we know g
′
′
′
(
i
) exists only for the interval for i
∈
[
-
6,6
] and outside the interval it will be zero. Then above equation can be written as (25)
-
1
u
-
6
g
-
6
′
′
′
+
u
-
6
g
-
6
′
′
′
+
⋯
+
u
5
g
5
′
′
′
+
u
6
g
6
′
′
′
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
.This implies (23).Theorem 5.
Forj
=
1,2
,
…
,
N, the nonlinear system of (22) becomes(26)
∑
i
=
-
6
+
j
6
+
j
-
1
u
i
g
i
-
j
′
′
′
=
h
3
f
x
j
,
U
x
j
,
U
′
x
j
.Proof.
Substitutingj
=
1 in (22), it becomes (27)
-
1
u
-
6
g
-
7
′
′
′
+
u
-
5
g
-
6
′
′
′
+
⋯
+
u
7
g
6
′
′
′
+
⋯
+
u
N
+
6
g
N
+
5
′
′
′
=
h
3
f
x
1
,
U
x
1
,
U
′
x
1
.
g
′
′
′
(
i
) is nonzero only for the interval for i
∈
[
-
6,6
] and outside the interval it will be zero. Then the above equation becomes(28)
-
1
u
-
5
g
-
6
′
′
′
+
u
-
4
g
-
5
′
′
′
+
⋯
+
u
4
g
5
′
′
′
+
u
5
h
6
′
′
′
=
h
3
f
x
1
,
U
x
1
,
U
′
x
1
.For j
=
2, (22) becomes (29)
-
1
u
-
6
g
-
8
′
′
′
+
u
-
5
g
-
7
′
′
′
+
u
-
4
g
-
6
′
′
′
+
⋯
+
u
4
g
2
′
′
′
+
u
5
g
3
′
′
′
+
u
6
g
4
′
′
′
+
u
7
g
5
′
′
′
u
8
g
6
′
′
′
+
⋯
+
u
N
+
5
g
N
+
3
′
′
′
+
u
N
+
6
g
N
+
4
′
′
′
=
h
3
f
x
2
,
U
x
2
,
U
′
x
2
.This implies(30)
-
1
u
-
4
g
-
6
′
′
′
+
u
-
3
g
-
5
′
′
′
+
⋯
+
u
7
g
5
′
′
′
+
u
8
g
6
′
′
′
=
h
4
f
x
2
,
U
x
2
,
U
′
x
2
.Similarly, we can find the expression for j
=
3,4
,
…
,
N; that is,(31)
-
1
u
-
3
g
-
6
′
′
′
+
u
-
4
g
-
5
′
′
′
+
⋯
+
u
8
g
5
′
′
′
+
u
9
g
6
′
′
′
=
h
3
f
x
3
,
U
x
3
,
U
′
x
3
,
-
1
u
-
2
g
-
6
′
′
′
+
u
-
1
g
-
5
′
′
′
+
⋯
+
u
9
g
5
′
′
′
+
u
10
g
6
′
′
′
=
h
3
f
x
4
,
U
x
4
,
U
′
x
4
,
⋮
-
1
u
N
-
6
g
-
6
′
′
′
+
u
N
-
5
g
-
5
′
′
′
+
⋯
+
u
N
+
5
g
5
′
′
′
+
u
N
+
6
g
6
′
′
′
=
h
3
f
x
N
,
U
x
N
,
U
′
x
N
.Hence, by combining (28), (30), and (31) we get (26).


## 4.2. Unstable Nonlinear System

The nonlinear system of (22) is equivalent to the following nonlinear system of N
+
1 equations with (N
+
13) unknowns {
u
i
}:(32)
A
1
U
=
F
u
,where A
1 is banded matrix of order (
N
+
1
)
×
(
N
+
13
), U is the unknown vector of order N
+
13, and F
(
u
) is the vector of order N
+
1 that depends on u. The matrix A and vectors U and F
(
u
) are given explicitly by(33)
A
1
=
-
1
g
p
q
′
′
′
q
-
p
-
6
N
+
1
×
N
+
13
,where p and q represent the row and column, respectively; that is, p
=
1,2
,
3
,
…
,
N
+
1 and q
=
1,2
,
3
,
…
,
N
+
13,(34)
F
u
=
h
3
f
x
0
,
U
x
0
,
U
′
x
0
,
…
,
h
3
f
x
N
,
U
x
N
,
U
′
x
N
T
,
(35)
U
=
u
-
6
,
u
-
5
,
z
-
4
,
…
,
u
N
+
4
,
u
N
+
5
,
u
N
+
6
T
,
U
′
x
j
=
∑
i
=
-
6
N
+
6
u
j
g
′
x
j
-
x
i
h
.
g
′
i is defined in (13) and g
(
i
)
=
g
i. System (32) is unstable and we need to make it stable to get unique solution. The detail for the stable system of nonlinear equations is given in next section.


## 4.3. Stable Nonlinear System

For unique solution of nonlinear systems (32), we need twelve more conditions. Three conditions can be attained from given boundary conditions for nonlinear systems of equations and remaining nine conditions are attained by setting some extrapolation method. The details of the given boundary conditions and extrapolation method are given below.


### 4.3.1. Approximated Boundary Condition

The given boundary conditions are(36)
U
0
=
α
1
,
U
′
0
=
α
2
,
U
N
=
α
3
.We see that first derivative is involved in the given boundary conditions, since approximation order of interpolating scheme (7) is eight, so we approximate derivative boundary conditions at end point with approximation order eight. The approximation of derivative conditions at end point is defined as(37)
U
′
0
=
N
840
-
2283
u
0
+
6720
u
1
-
11760
u
2
+
15680
u
3
-
14700
u
4
+
9408
u
5
-
3920
u
6
+
960
u
7
-
105
u
8
+
O
h
8
.


### 4.3.2. Imposed Boundary Conditions

Remaining nine conditions for the nonlinear systems (32) to get stable systems for the solution of (1) are obtained by setting the following extrapolation method.We define five conditions at left end points and four conditions at the right end points. Since subdivision scheme reproduces seven-degree polynomials, so we define boundary conditions of order eight for solution of (32). For simplicity only the left end points u
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 are discussed and the values of right end points u
N
+
2
,
u
N
+
3
,
u
N
+
4
,
u
N
+
5 can be treated similarly.The valuesu
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 can be determined by the septic polynomial R
(
x
) interpolating (
x
i
,
u
i
)
,
0
≤
i
≤
5. Precisely, we have (38)
u
-
i
=
R
-
x
i
,
i
=
1,2
,
3,4
,
5
, where (39)
R
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
U
x
i
-
r
. Since by (16), U
(
x
i
)
=
u
i for i
=
1,2
,
3,4
,
5 then by replacing x
i by -
x
i, we have (40)
R
-
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
u
-
i
+
r
. Hence, the following boundary conditions can be employed at the left end:(41)
∑
r
=
0
8
8
r
-
1
r
u
-
i
+
r
=
0
,
i
=
5,4
,
3,2
,
1
.Similarly for the right end, we can define u
i
=
R
(
-
x
i
), i
=
N
+
2, N
+
3, N
+
4, N
+
5, and (42)
R
x
i
=
∑
r
=
0
8
8
r
-
1
r
+
1
u
i
-
r
. So we have the following boundary conditions at the right end:(43)
∑
r
=
0
8
8
r
-
1
r
u
i
-
r
=
0
,
i
=
N
+
2
,
N
+
3
,
N
+
4
,
N
+
5
.Finally, we get the following new system of (
N
+
13
) linear equations with (
N
+
13
) unknowns {
u
i
}, in which N
+
1 equations are obtained from (22), three equations from boundary conditions (36), and nine from boundary conditions (41) and (43).Hence, the stable nonlinear system of equations is defined as(44)
A
U
=
G
u
,where the matrix A is given by(45)
A
=
A
0
T
,
A
1
T
,
A
2
T
T
,and the matrix A
1 is defined in (33). The matrix (
A
0
)
7
×
(
N
+
13
) is constructed by taking first five rows from (41), sixth row from (37), and the last row U
(
0
) from (36). Hence,(46)
A
0
=
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
⋯
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
⋯
0
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
⋯
0
0
0
0
0
0
0
0
-
2283
N
840
6720
N
840
-
11760
N
840
15680
N
840
-
14700
N
840
9408
N
840
-
3920
N
840
960
N
840
-
105
N
840
0
⋯
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
⋯
0
0
. The matrix (
A
2
)
6
×
(
N
+
13
) is constructed by taking first row U
(
N
) from (36) and last five rows are obtained from (43). Hence,(47)
A
2
=
0
0
⋯
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
⋯
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
0
⋯
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
. The column vector U is defined in (35) and G is defined as(48)
G
u
=
0,0
,
0,0
,
0
,
U
′
0
,
U
0
,
F
T
u
,
U
1
,
0,0
,
0,0
T
,where F
(
u
) is defined by (34).


### 4.3.3. Nonsingularity of a Matrix

The nonsingularity of coefficient matrixA, which is defined in (45), can be checked by different methods.Since the determinant of matrixA is nonzero for N
≤
2000 then the nonlinear system of equations give solution for N
≤
2000.It is also observed that all the eigenvalues are nonzero forN
≤
2000. Hence, by [14] matrix A is nonsingular. For large N
>
2000 the matrix may or may not be nonsingular.The coefficient matrixA is neither symmetric nor diagonally dominant. Though it can be proved that A is nonsingular/invertible, the matrix is almost symmetric except the first and last few rows and columns due to its boundary treatment for large value of N. Therefore, we first consider the symmetric part of it, that is, the square band matrix B of order N
+
1 defined as (49)
-
1
×
B
=
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
g
4
′
′
′
⋯
⋯
0
0
0
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
⋯
⋯
0
0
0
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
⋯
⋯
0
0
0
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
⋯
⋯
0
0
0
g
-
3
′
′
′
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
⋯
⋯
0
0
0
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
⋯
⋯
g
2
′
′
′
g
2
′
′
′
g
3
′
′
′
0
0
0
0
⋯
⋯
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
0
0
0
0
⋯
⋯
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
.Kılıç and Stanica [15] presented a method to find the inverse of banded matrix of order n by using the LU factorization of the banded matrix. As B is a band matrix of order (
N
+
1
) so by LU factorization method its inverse exists.


## 4.3.1. Approximated Boundary Condition

The given boundary conditions are(36)
U
0
=
α
1
,
U
′
0
=
α
2
,
U
N
=
α
3
.We see that first derivative is involved in the given boundary conditions, since approximation order of interpolating scheme (7) is eight, so we approximate derivative boundary conditions at end point with approximation order eight. The approximation of derivative conditions at end point is defined as(37)
U
′
0
=
N
840
-
2283
u
0
+
6720
u
1
-
11760
u
2
+
15680
u
3
-
14700
u
4
+
9408
u
5
-
3920
u
6
+
960
u
7
-
105
u
8
+
O
h
8
.


## 4.3.2. Imposed Boundary Conditions

Remaining nine conditions for the nonlinear systems (32) to get stable systems for the solution of (1) are obtained by setting the following extrapolation method.We define five conditions at left end points and four conditions at the right end points. Since subdivision scheme reproduces seven-degree polynomials, so we define boundary conditions of order eight for solution of (32). For simplicity only the left end points u
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 are discussed and the values of right end points u
N
+
2
,
u
N
+
3
,
u
N
+
4
,
u
N
+
5 can be treated similarly.The valuesu
-
5
,
u
-
4
,
u
-
3
,
u
-
2
,
u
-
1 can be determined by the septic polynomial R
(
x
) interpolating (
x
i
,
u
i
)
,
0
≤
i
≤
5. Precisely, we have (38)
u
-
i
=
R
-
x
i
,
i
=
1,2
,
3,4
,
5
, where (39)
R
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
U
x
i
-
r
. Since by (16), U
(
x
i
)
=
u
i for i
=
1,2
,
3,4
,
5 then by replacing x
i by -
x
i, we have (40)
R
-
x
i
=
∑
r
=
1
8
8
r
-
1
r
+
1
u
-
i
+
r
. Hence, the following boundary conditions can be employed at the left end:(41)
∑
r
=
0
8
8
r
-
1
r
u
-
i
+
r
=
0
,
i
=
5,4
,
3,2
,
1
.Similarly for the right end, we can define u
i
=
R
(
-
x
i
), i
=
N
+
2, N
+
3, N
+
4, N
+
5, and (42)
R
x
i
=
∑
r
=
0
8
8
r
-
1
r
+
1
u
i
-
r
. So we have the following boundary conditions at the right end:(43)
∑
r
=
0
8
8
r
-
1
r
u
i
-
r
=
0
,
i
=
N
+
2
,
N
+
3
,
N
+
4
,
N
+
5
.Finally, we get the following new system of (
N
+
13
) linear equations with (
N
+
13
) unknowns {
u
i
}, in which N
+
1 equations are obtained from (22), three equations from boundary conditions (36), and nine from boundary conditions (41) and (43).Hence, the stable nonlinear system of equations is defined as(44)
A
U
=
G
u
,where the matrix A is given by(45)
A
=
A
0
T
,
A
1
T
,
A
2
T
T
,and the matrix A
1 is defined in (33). The matrix (
A
0
)
7
×
(
N
+
13
) is constructed by taking first five rows from (41), sixth row from (37), and the last row U
(
0
) from (36). Hence,(46)
A
0
=
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
⋯
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
⋯
0
0
0
0
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
⋯
0
0
0
0
0
0
0
0
-
2283
N
840
6720
N
840
-
11760
N
840
15680
N
840
-
14700
N
840
9408
N
840
-
3920
N
840
960
N
840
-
105
N
840
0
⋯
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
⋯
0
0
. The matrix (
A
2
)
6
×
(
N
+
13
) is constructed by taking first row U
(
N
) from (36) and last five rows are obtained from (43). Hence,(47)
A
2
=
0
0
⋯
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
⋯
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
0
⋯
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
0
⋯
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
0
0
0
⋯
0
0
0
1
-
8
28
-
56
70
-
56
28
-
8
1
0
0
. The column vector U is defined in (35) and G is defined as(48)
G
u
=
0,0
,
0,0
,
0
,
U
′
0
,
U
0
,
F
T
u
,
U
1
,
0,0
,
0,0
T
,where F
(
u
) is defined by (34).


## 4.3.3. Nonsingularity of a Matrix

The nonsingularity of coefficient matrixA, which is defined in (45), can be checked by different methods.Since the determinant of matrixA is nonzero for N
≤
2000 then the nonlinear system of equations give solution for N
≤
2000.It is also observed that all the eigenvalues are nonzero forN
≤
2000. Hence, by [14] matrix A is nonsingular. For large N
>
2000 the matrix may or may not be nonsingular.The coefficient matrixA is neither symmetric nor diagonally dominant. Though it can be proved that A is nonsingular/invertible, the matrix is almost symmetric except the first and last few rows and columns due to its boundary treatment for large value of N. Therefore, we first consider the symmetric part of it, that is, the square band matrix B of order N
+
1 defined as (49)
-
1
×
B
=
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
g
4
′
′
′
⋯
⋯
0
0
0
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
g
3
′
′
′
⋯
⋯
0
0
0
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
⋯
⋯
0
0
0
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
⋯
⋯
0
0
0
g
-
3
′
′
′
g
-
2
′
′
′
g
-
1
′
′
′
g
0
′
′
′
⋯
⋯
0
0
0
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
⋯
0
0
0
0
⋯
⋯
g
2
′
′
′
g
2
′
′
′
g
3
′
′
′
0
0
0
0
⋯
⋯
g
0
′
′
′
g
1
′
′
′
g
2
′
′
′
0
0
0
0
⋯
⋯
g
-
1
′
′
′
g
0
′
′
′
g
1
′
′
′
.Kılıç and Stanica [15] presented a method to find the inverse of banded matrix of order n by using the LU factorization of the banded matrix. As B is a band matrix of order (
N
+
1
) so by LU factorization method its inverse exists.


## 4.4. Iterative Algorithm and Its Convergence

In this section, we propose an iterative algorithm and discuss its convergence.


### 4.4.1. Iterative Algorithm Based on Basis Function

The iterative algorithm is based on basis function of the subdivision scheme (7) as defined in the following three steps.First Step (initial approximation). Initial approximation is important because numerical solution depends on the initial approximation. We define the process for finding the initial approximation as follows.Let initial approximate solutionU
0 be the solution of the following linear system:(50)
A
U
0
=
G
0
,where (51)
G
0
=
0,0
,
0,0
,
0
,
y
′
a
,
y
a
,
p
0
,
p
1
,
p
2
,
…
,
p
N
,
y
b
,
0,0
,
0,0
T
,
p
i
=
h
3
p
x
i
,
m
i
,
S
,
i
=
0,1
,
2
,
…
,
N
,
a
≤
x
≤
b
,
m
i
=
y
a
+
i
h
y
b
-
y
a
b
-
a
,
S
=
y
b
-
y
a
.
G
0 is the initial linear approximation of the nonlinear vector G
(
u
). By solving the linear system of (50) we get initial approximate solution.Second Step (numerical solution). The numerical solutions U
∗ of the nonlinear system are obtained by using simple iterative scheme:(52)
A
U
m
+
1
=
G
U
m
,
m
=
0,1
,
2,3
,
…
.Third Step (stopping condition). The above iterative processes will terminate when the following condition is satisfied(53)
U
m
-
U
m
-
1
≤
t
o
l
,where tol is supposed value; that is, t
o
l
=
1
0
-
6. The convergence of the above iterative algorithm is guaranteed by the following propositions. The solutions of linear system of (50) and (52) are obtained by Gaussian eliminations method.Proposition 6.
The successive solutions{
U
m
} for the nonlinear system (44) generated by the iterative algorithm (52) linearly converge to the solution U
∗ provided that M
0 and M
1 are Lipschitz constants and step size h is small. That is,(54)
A
-
1
M
0
h
3
+
867307
212370
M
1
h
2
≤
1
.Proof.
LetU
∗ and U
(
m
) be the solutions of the nonlinear system (44) then by definition, for small h, we have (55)
A
U
∗
=
G
U
∗
,
A
U
m
+
1
=
G
U
m
. Let the error vector be defined as e
(
m
)
=
U
m
-
U
∗ at mth iteration which satisfies (56)
A
U
m
+
1
-
A
U
∗
=
G
U
m
-
G
U
∗
,
G
U
m
+
1
-
U
∗
=
G
U
m
-
G
U
∗
,
A
e
m
+
1
=
G
U
m
-
G
U
∗
. For i
=
0,1
,
2
,
…
,
N, (57)
D
3
e
i
m
+
1
=
F
U
m
-
F
U
∗
i
,by using the mean value theorem, which is stated as “if a function f
(
x
,
y
,
z
) is continuously differentiable in an open set of R
3 containing points (
x
1
,
y
1
,
z
1
) and (
x
2
,
y
2
,
z
2
) and the line segment connecting them, then an equation (58)
f
x
2
,
y
2
,
z
2
-
f
x
1
,
y
1
,
z
1
=
f
x
′
r
,
s
,
t
x
2
-
x
1
+
f
y
′
r
,
s
,
t
y
2
-
y
1
+
f
z
′
r
,
s
,
t
z
2
-
z
1 is valid for the interior point (
a
,
b
,
c
) of the segment.” We have (59)
D
3
e
i
m
+
1
=
f
x
i
,
U
i
m
,
U
′
m
-
f
x
i
,
U
i
∗
,
U
′
∗
. The above equation can be written as (by using mean value theorem) (60)
D
3
e
i
m
+
1
=
f
x
∗
x
i
-
x
i
+
f
y
∗
U
i
m
-
U
i
∗
+
f
y
′
∗
U
′
m
-
U
′
∗
, and by using the definition of error vector, we have (61)
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
e
i
′
m
,
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
D
1
e
i
m
, where D
3 and D
1 are the derivative difference operators defined as (62)
D
1
f
i
=
1
2973180
h
-
5
f
i
+
6
-
f
i
-
6
+
1024
f
i
+
5
-
f
i
-
5
+
13225
f
i
+
4
-
f
i
-
4
-
199680
f
i
+
3
-
f
i
-
3
+
1141695
f
i
+
2
-
f
i
-
2
-
4715520
f
i
+
1
-
f
i
-
1
,
D
3
f
i
=
1
1879920
h
3
-
225
f
i
+
6
-
f
i
-
6
+
11520
f
i
+
5
-
f
i
-
5
-
10952
f
i
+
4
-
f
i
-
4
-
476928
f
i
+
3
-
f
i
-
3
+
3047987
f
i
+
2
-
f
i
-
2
-
4677632
f
i
+
1
-
f
i
-
1
. This implies (63)
D
3
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. Since e
i
=
e
N
-
i
=
0
,
i
=
0
,
-
1
,
-
2
,
…
,
-
6, therefore, we have (64)
A
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. This can be written as (65)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By taking norm on both sides, we get (66)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By using the definition of Lipschitz condition, we get(67)
e
i
m
+
1
≤
A
-
1
h
3
M
0
e
i
m
+
h
2
M
1
D
1
e
i
m
.This implies (68)
e
i
m
+
1
e
i
m
≤
A
-
1
h
3
M
0
+
h
2
M
1
D
1
, which is equivalent to (69)
e
i
m
+
1
e
i
m
≈
h
2
M
1
A
-
1
D
1
. The results follow immediately from this inequality and the following fact: (70)
D
1
=
867307
212370
.A simple approximation of condition by omitting the cubic term is (71)
h
≤
212370
867307
M
1
-
1
A
-
1
-
1
1
/
2
.This completes the proof.Proposition 7.
Iff satisfies the Lipschitz condition (3) and Lipschitz constants M
0, M
1, and M
2 and mesh size h are small enough then nonlinear system (32) has a unique solution U
∗. A sufficient condition for the existence of a solution is given by (54).Proof.
From the proof of previous proposition, we observe that if (54) holds, then an inequality similar to that of (67) holds, which implies that the sequence {
U
m
} is contracting and hence converges. The limit U
∗ of (52) also satisfies (44) due to the continuity of the right hand side function f
x
,
y
,
y
′
,
y
′
′.Remark 8.
The numerical complexity for the solution of the linear system (52), where the matrix A is almost band matrix with half band 7, using Gaussian elimination method is about 49
(
N
+
13
) multiplications. The number of complexities depends upon the efficient boundary treatment. If more efficient boundary treatment is constructed then number of complexities will be reduced.


## 4.4.1. Iterative Algorithm Based on Basis Function

The iterative algorithm is based on basis function of the subdivision scheme (7) as defined in the following three steps.First Step (initial approximation). Initial approximation is important because numerical solution depends on the initial approximation. We define the process for finding the initial approximation as follows.Let initial approximate solutionU
0 be the solution of the following linear system:(50)
A
U
0
=
G
0
,where (51)
G
0
=
0,0
,
0,0
,
0
,
y
′
a
,
y
a
,
p
0
,
p
1
,
p
2
,
…
,
p
N
,
y
b
,
0,0
,
0,0
T
,
p
i
=
h
3
p
x
i
,
m
i
,
S
,
i
=
0,1
,
2
,
…
,
N
,
a
≤
x
≤
b
,
m
i
=
y
a
+
i
h
y
b
-
y
a
b
-
a
,
S
=
y
b
-
y
a
.
G
0 is the initial linear approximation of the nonlinear vector G
(
u
). By solving the linear system of (50) we get initial approximate solution.Second Step (numerical solution). The numerical solutions U
∗ of the nonlinear system are obtained by using simple iterative scheme:(52)
A
U
m
+
1
=
G
U
m
,
m
=
0,1
,
2,3
,
…
.Third Step (stopping condition). The above iterative processes will terminate when the following condition is satisfied(53)
U
m
-
U
m
-
1
≤
t
o
l
,where tol is supposed value; that is, t
o
l
=
1
0
-
6. The convergence of the above iterative algorithm is guaranteed by the following propositions. The solutions of linear system of (50) and (52) are obtained by Gaussian eliminations method.Proposition 6.
The successive solutions{
U
m
} for the nonlinear system (44) generated by the iterative algorithm (52) linearly converge to the solution U
∗ provided that M
0 and M
1 are Lipschitz constants and step size h is small. That is,(54)
A
-
1
M
0
h
3
+
867307
212370
M
1
h
2
≤
1
.Proof.
LetU
∗ and U
(
m
) be the solutions of the nonlinear system (44) then by definition, for small h, we have (55)
A
U
∗
=
G
U
∗
,
A
U
m
+
1
=
G
U
m
. Let the error vector be defined as e
(
m
)
=
U
m
-
U
∗ at mth iteration which satisfies (56)
A
U
m
+
1
-
A
U
∗
=
G
U
m
-
G
U
∗
,
G
U
m
+
1
-
U
∗
=
G
U
m
-
G
U
∗
,
A
e
m
+
1
=
G
U
m
-
G
U
∗
. For i
=
0,1
,
2
,
…
,
N, (57)
D
3
e
i
m
+
1
=
F
U
m
-
F
U
∗
i
,by using the mean value theorem, which is stated as “if a function f
(
x
,
y
,
z
) is continuously differentiable in an open set of R
3 containing points (
x
1
,
y
1
,
z
1
) and (
x
2
,
y
2
,
z
2
) and the line segment connecting them, then an equation (58)
f
x
2
,
y
2
,
z
2
-
f
x
1
,
y
1
,
z
1
=
f
x
′
r
,
s
,
t
x
2
-
x
1
+
f
y
′
r
,
s
,
t
y
2
-
y
1
+
f
z
′
r
,
s
,
t
z
2
-
z
1 is valid for the interior point (
a
,
b
,
c
) of the segment.” We have (59)
D
3
e
i
m
+
1
=
f
x
i
,
U
i
m
,
U
′
m
-
f
x
i
,
U
i
∗
,
U
′
∗
. The above equation can be written as (by using mean value theorem) (60)
D
3
e
i
m
+
1
=
f
x
∗
x
i
-
x
i
+
f
y
∗
U
i
m
-
U
i
∗
+
f
y
′
∗
U
′
m
-
U
′
∗
, and by using the definition of error vector, we have (61)
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
e
i
′
m
,
D
3
e
i
m
+
1
=
f
y
∗
e
i
m
+
f
y
′
∗
D
1
e
i
m
, where D
3 and D
1 are the derivative difference operators defined as (62)
D
1
f
i
=
1
2973180
h
-
5
f
i
+
6
-
f
i
-
6
+
1024
f
i
+
5
-
f
i
-
5
+
13225
f
i
+
4
-
f
i
-
4
-
199680
f
i
+
3
-
f
i
-
3
+
1141695
f
i
+
2
-
f
i
-
2
-
4715520
f
i
+
1
-
f
i
-
1
,
D
3
f
i
=
1
1879920
h
3
-
225
f
i
+
6
-
f
i
-
6
+
11520
f
i
+
5
-
f
i
-
5
-
10952
f
i
+
4
-
f
i
-
4
-
476928
f
i
+
3
-
f
i
-
3
+
3047987
f
i
+
2
-
f
i
-
2
-
4677632
f
i
+
1
-
f
i
-
1
. This implies (63)
D
3
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. Since e
i
=
e
N
-
i
=
0
,
i
=
0
,
-
1
,
-
2
,
…
,
-
6, therefore, we have (64)
A
e
i
m
+
1
=
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. This can be written as (65)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By taking norm on both sides, we get (66)
e
i
m
+
1
=
A
-
1
h
3
f
y
∗
e
i
m
+
h
2
f
y
′
∗
D
1
e
i
m
. By using the definition of Lipschitz condition, we get(67)
e
i
m
+
1
≤
A
-
1
h
3
M
0
e
i
m
+
h
2
M
1
D
1
e
i
m
.This implies (68)
e
i
m
+
1
e
i
m
≤
A
-
1
h
3
M
0
+
h
2
M
1
D
1
, which is equivalent to (69)
e
i
m
+
1
e
i
m
≈
h
2
M
1
A
-
1
D
1
. The results follow immediately from this inequality and the following fact: (70)
D
1
=
867307
212370
.A simple approximation of condition by omitting the cubic term is (71)
h
≤
212370
867307
M
1
-
1
A
-
1
-
1
1
/
2
.This completes the proof.Proposition 7.
Iff satisfies the Lipschitz condition (3) and Lipschitz constants M
0, M
1, and M
2 and mesh size h are small enough then nonlinear system (32) has a unique solution U
∗. A sufficient condition for the existence of a solution is given by (54).Proof.
From the proof of previous proposition, we observe that if (54) holds, then an inequality similar to that of (67) holds, which implies that the sequence {
U
m
} is contracting and hence converges. The limit U
∗ of (52) also satisfies (44) due to the continuity of the right hand side function f
x
,
y
,
y
′
,
y
′
′.Remark 8.
The numerical complexity for the solution of the linear system (52), where the matrix A is almost band matrix with half band 7, using Gaussian elimination method is about 49
(
N
+
13
) multiplications. The number of complexities depends upon the efficient boundary treatment. If more efficient boundary treatment is constructed then number of complexities will be reduced.


## 5. Error Analysis

From the approximation properties of the basis functiong
(
x
), it is shown that the collocation method (16) with nonic precision treatments at end points has at least the power of approximation O
(
h
3
). Here we present our main results for error estimation. Proof of these results is similar to the proof of Proposition in [9].Proposition 9.
Suppose the exact solutiony
(
x
)
∈
C
3
[
0,1
] and {
u
i
} are obtained by (67) then absolute error by interpolating collocation algorithm is (72)
e
r
r
x
∞
=
U
l
x
-
y
l
x
∞
=
O
h
3
-
l
,
l
=
0,1
,
2,3
,where l denotes the order of derivative.Proof.
Since the order of approximation of subdivision scheme (2) is ten so by direct calculation (third left eigenvector), we can find derivative of smooth function y
(
x
) as (73)
y
′
′
′
x
j
=
1
1879920
h
3
-
225
y
x
j
+
6
h
-
y
x
j
-
6
h
+
11520
y
x
j
+
5
h
-
y
x
j
-
5
h
-
10952
y
x
j
+
4
h
+
y
x
j
-
4
h
-
476928
y
x
j
+
3
h
+
y
x
j
-
3
h
+
3047987
y
x
j
+
2
h
-
y
x
j
-
2
h
-
4677632
y
x
j
+
h
-
y
x
j
-
h
+
O
h
8
. This can be written as(74)
y
′
′
′
x
j
=
1
1879920
h
3
-
225
y
j
+
6
-
y
j
-
6
+
11520
y
j
+
5
-
y
j
-
5
-
10952
y
j
+
4
-
y
j
-
4
-
476928
y
j
+
3
-
y
j
-
3
+
3047987
y
j
+
2
-
y
j
-
2
-
4677632
y
j
+
1
-
y
j
-
1
+
O
h
8
.Similarly, we have(75)
U
′
′
′
x
j
=
1
1879920
h
3
-
225
u
j
+
6
-
u
j
-
6
+
11520
u
j
+
5
-
u
j
-
5
-
10952
u
j
+
4
-
u
j
-
4
-
476928
u
j
+
3
-
u
j
-
3
+
3047987
u
j
+
2
-
u
j
-
2
-
4677632
u
j
+
1
-
u
j
-
1
+
O
h
8
.If we define error function e
(
x
)
=
U
(
x
)
-
y
(
x
) and error vectors at the nodes by (76)
e
x
j
=
U
x
j
-
y
x
j
+
j
h
,
-
6
≤
j
≤
N
+
6
,or equivalently e
j
=
Z
j
-
y
j
,
-
6
≤
j
≤
N
+
6, then this implies (77)
e
j
′
=
U
j
′
-
y
j
′
,
e
j
′
′
=
U
j
′
′
-
y
j
′
′
,
e
j
′
′
′
=
U
j
′
′
′
-
y
j
′
′
′
. By subtracting (75) from (74), we get(78)
y
j
′
′
′
-
U
j
′
′
′
=
1
1879920
h
3
-
225
y
j
+
6
-
y
j
-
6
-
u
j
+
6
-
u
j
-
6
+
11520
y
j
+
5
-
y
j
-
5
-
u
j
+
5
-
u
j
-
5
-
10952
y
j
+
4
-
y
j
-
4
-
u
j
+
4
-
u
j
-
4
-
476928
y
j
+
3
-
y
j
-
3
-
u
j
+
3
-
u
j
-
3
+
3047987
y
j
+
2
-
y
j
-
2
-
u
j
+
2
-
u
j
-
2
-
4677632
y
j
+
1
-
y
j
-
1
-
u
j
+
1
-
u
j
-
1
+
O
h
8
.This implies(79)
e
j
′
′
′
=
1
1879920
h
3
-
225
e
j
+
6
-
e
j
-
6
+
11520
e
j
+
5
-
e
j
-
5
-
10952
e
j
+
4
-
e
j
-
4
-
476928
e
j
+
3
-
e
j
-
3
+
3047987
e
j
+
2
-
e
j
-
2
-
4677632
e
j
+
1
-
e
j
-
1
+
O
h
8
.From (1), (16), and (79) and by assuming the eighth order boundary treatments at the end points, we have(80)
e
j
′
′
′
=
a
j
e
j
+
b
j
e
j
′
,
0
≤
i
≤
N
,
e
j
=
max
0
≤
k
≤
6
⁡
e
k
O
h
8
,
-
6
≤
j
≤
0
max
N
-
6
≤
k
≤
N
⁡
e
k
O
h
8
,
N
≤
j
≤
N
+
6
,where j
=
0,1
,
…
,
N, (81)
a
j
=
f
y
t
j
,
y
j
∗
,
y
j
′
∗
,
b
j
=
f
y
′
t
j
,
y
j
∗
,
y
j
′
∗
,
y
j
∗
=
y
j
+
θ
j
e
j
,
y
j
′
∗
=
y
j
′
+
θ
j
e
j
′
,
0
≤
θ
j
≤
1
.Using the results (79) and(82)
D
1
U
i
=
1
2973180
h
-
5
u
i
+
6
-
u
i
-
6
+
1024
u
i
+
5
-
u
i
-
5
+
13225
u
i
+
4
-
u
i
-
4
-
199680
u
i
+
3
-
u
i
-
3
+
1141695
u
i
+
2
-
u
i
-
2
-
4715520
u
i
+
1
-
u
i
-
1
, it can be concluded that relation (80) is equivalent to (83)
A
+
O
h
6
-
O
h
3
-
D
1
O
h
E
=
O
h
8
E
,where E
=
(
e
-
6
,
e
-
5
,
…
,
e
5
,
e
6
).
Hence, for smallh, the coefficient matrix A
+
O
(
h
2
) will be invertible; thus using the standard result from algebra and effect of A
-
1, we have the following estimate: (84)
E
≤
A
-
1
1
-
O
h
2
O
h
8
=
O
h
3
. This completes the result.


## 6. Examples, Comparison, and Conclusion

In this section, we use subdivision based collocation algorithm to find the solution of some nonlinear third-order boundary value problems. We present numerical results in table format along with their graphical representations. We also give comparison of the results obtained by our algorithm and the results computed by existing algorithms. We end this section with precise conclusion.


### 6.1. Numerical Examples

We find the approximate solutions of the following nonlinear problems to check the accuracy and convergence of subdivision based iterative collocation algorithm.Example 1.
The nonlinear boundary value problem is(85)
y
′
′
′
=
-
2
e
-
3
y
+
4
1
+
x
-
3
,with boundary conditions(86)
y
0
=
0
,
y
′
0
=
1
,
y
1
=
ln
⁡
2
.Example 2.
The nonlinear boundary value problem is(87)
y
′
′
′
x
=
e
-
x
y
2
xsubject to the boundary conditions:(88)
y
0
=
1
,
y
′
0
=
1
,
y
1
=
e
.Example 3.
The nonlinear boundary value problem is(89)
y
′
′
′
x
=
-
e
x
y
2
xsubject to the boundary conditions: (90)
y
0
=
1
,
y
′
0
=
-
1
,
y
1
=
1
e
.The exact solutions of problems (85), (87), and (89) are y
=
ln
⁡
(
1
+
x
), y
=
e
x, and y
=
e
-
x, respectively.Example 4.
We consider the third-order ordinary differential equation:(91)
y
′
′
′
=
y
-
k
,
y
⩾
0
,where k is constant. The initial conditions imposed by Tanner [16] are(92)
y
0
=
1
,
y
′
0
=
0
.The problem is closed by the boundary condition:(93)
y
r
=
0
,and the problem becomes singular at y
=
0. The boundary condition is imposed [17] as(94)
y
r
=
ϵ
,where r is a constant satisfying r
>
0. The analytical solution of (91) by [17](95)
y
x
=
1
6
6
-
7
x
2
+
6
ϵ
x
2
+
x
3
.Case 1 (see [17]).
Fork
=
0 problem (91) takes the form(96)
y
′
′
′
=
y
0
=
1
.
We solve (96) along with conditions (92) and (94) by letting r
=
1; then ϵ
=
11
/
12 and we obtain the results which also support our algorithm. That is, the numerical results have the order of approximation O
(
h
3
). The numerical results are tabulated in Table 4 and graphical representation of these results is shown in Figure 4. These results are obtained after first iteration level. The maximum absolute error is 6.1250
×
1
0
-
3.Case 2 (see [17]).
Fork
=
1
/
2 problem (91) takes the form(97)
y
′
′
′
=
y
-
1
/
2
,with the boundary conditions y
(
0
)
=
1, y
′
0
=
0, and y
(
r
)
=
ϵ. The numerical solution of (97) is obtained by using the proposed numerical algorithm. The solution after third iteration level is presented in Table 5 and their graphical representation is presented in Figure 5. The maximum absolute error for this problem is 6.4004
×
1
0
-
3.Case 3 (see [17, 18]).
In this case, we numerically solve problem (91) for k
=
2; that is,(98)
y
′
′
′
=
y
-
2together with the boundary conditions y
(
0
)
=
1, y
′
0
=
0, and y
(
r
)
=
ϵ. Its exact solution is given in (95). The numerical solution of (98) with given and imposed condition at node points is tabulated in Table 6 and graphical representation is given in Figure 6. The maximum absolute error for this problem is 4.93276
×
1
0
-
3.


### 6.2. Comparison and Discussion

For Examples1, 2, and 3, we use the iterative collocation algorithm described in Section 4, for h
=
1
0
-
1
,
2
0
-
1
,
5
0
-
1 (i.e., for N
=
10,20,50) and t
o
l
=
1
0
-
6 along with eighth-order boundary treatment at end points, to get solutions of nonlinear boundary value problems. The numerical results are obtained after third iteration with condition (53).(i)
The numerical solutions of problems (85), (87), and (89) are presented in Tables 1, 2, and 3, respectively.
(ii)
Caglar et al. [5] solved problem (85) by fourth-degree B-spline algorithm. The maximum absolute errors obtained by the proposed algorithm and by [5] are 4.77
×
1
0
-
3 and 5.80
×
1
0
-
2, respectively. The graphical comparison between exact and approximate solutions is shown in Figure 1. We observe that numerical results obtained by proposed algorithm are better than the results of [5].
(iii)
Hasan [3] solved problem (87) by modified Adomian decomposition method (MADM). We also solve this problem by subdivision based collocation algorithm. Here we observe that the order of approximation by proposed and MADM algorithms is the same (i.e., O
(
h
3
)). The graphical comparison between exact and approximate solutions is shown in Figure 2.
(iv)
The comparison between exact and approximate solutions of problem (89) is given in Figure 3.In Example 4, we consider the problem related to thin film flows. We solve (91) by assuming different values of k and we observer from the numerical results tabulated in Tables 4, 5, and 6 that accuracy of the approximate solution is O
(
h
3
). The numerical results are obtained after first and third iterations for k
=
0 and k
=
1
/
2
,
1 with condition (53), respectively.Table 1
Numerical results of Example1: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
Error by method [5]
0.0
0
0
0
0
0.1
0.0953101798
0.0957100706
0.0003998908
0.009
0.2
0.1823215568
0.1836519852
0.0013304284
0.013
0.3
0.2623642645
0.2648256813
0.0024614168
0.031
0.4
0.3364722366
0.3400108156
0.0035385790
0.045
0.5
0.4054651081
0.4098275640
0.0043624559
0.054
0.6
0.4700036292
0.4747771017
0.0047734725
0.058
0.7
0.5306282511
0.5352702927
0.0046420416
0.058
0.8
0.5877866649
0.5916474874
0.0038608225
0.053
0.9
0.6418538862
0.6441934825
0.0023395963
0.044
1.0
0.6931471806
0.6931471806
0
0.000Table 2
Numerical results of Example2: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
1.105170918
1.1056664532
0.000495535
0.2
1.221402758
1.2232054699
0.001802712
0.3
1.349858808
1.3534912400
0.003632432
0.4
1.491824698
1.4974920669
0.005667369
0.5
1.648721271
1.6562805667
0.007559296
0.6
1.822118800
1.8310443828
0.008925583
0.7
2.013752707
2.0230974700
0.009344763
0.8
2.225540928
2.2338920118
0.008351084
0.9
2.459603111
2.4650310392
0.005427928
1.0
2.718281828
2.7182818285
0Table 3
Numerical results of Example3: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9048374180
0.9045439989
0.0002934191
0.2
0.8187307531
0.8177120092
0.0010187439
0.3
0.7408182207
0.7388586762
0.0019595445
0.4
0.6703200460
0.6674013997
0.0029186463
0.5
0.6065306597
0.6028146540
0.0037160057
0.6
0.5488116363
0.5446245817
0.0041870544
0.7
0.4965853038
0.4924038346
0.0041814692
0.8
0.4493289641
0.4457666384
0.0035623257
0.9
0.4065696597
0.4043640577
0.0022056020
1.0
0.3678794412
0.3678794412
0Table 4
Numerical results of Example4: k
=
0 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980416667
0.0003750000
0.2
0.9913333333
0.9926666667
0.0013333334
0.3
0.9820000000
0.9846250000
0.0026250000
0.4
0.9706666667
0.9746666667
0.0040000000
0.5
0.9583333333
0.9635416667
0.0052083334
0.6
0.9460000000
0.9520000000
0.0060000000
0.7
0.9346666667
0.9407916667
0.0061250000
0.8
0.9253333333
0.9306666667
0.0053333334
0.9
0.9190000000
0.9223750000
0.0033750000
1.0
0.9166666667
0.9166666667
0Table 5
Numerical results of Example4: k
=
1
/
2 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980508018
0.0003841351
0.2
0.9913333333
0.9927028673
0.0013695340
0.3
0.9820000000
0.9847044600
0.0013695340
0.4
0.9706666667
0.9748013302
0.0041346635
0.5
0.9583333333
0.9637358441
0.0054025108
0.6
0.9460000000
0.9522463844
0.0062463844
0.7
0.9346666667
0.9410670255
0.0064003588
0.8
0.9253333333
0.9309274865
0.0055941532
0.9
0.9190000000
0.9225533666
0.0035533666
1.0
0.9166666667
0.9166666692
0Table 6
Numerical results of Example4: k
=
2 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980024045
0.0003357378
0.2
0.9913333333
0.9925110306
0.0011776973
0.3
0.9820000000
0.9842830603
0.0022830603
0.4
0.9706666667
0.9740862436
0.0034195769
0.5
0.9583333333
0.9627028681
0.0043695348
0.6
0.9460000000
0.9509327630
0.0049327630
0.7
0.9346666667
0.9395952385
0.0049285718
0.8
0.9253333333
0.9295297470
0.0041964137
0.9
0.9190000000
0.9215950588
0.0025950588
1.0
0.9166666667
0.9166666667
0Figure 1
Comparison of the analytic and approximate solution of Example1 by proposed algorithm and [5]. In this figure solid line shows exact solution, dotted lines show approximate solution by proposed algorithm, and dashed lines show the solution obtained by [5].Figure 2
Comparison of the analytic and approximate solution of Example2 by proposed algorithm and [3]. In this figure solid line shows exact solution, dotted lines show approximate solution by proposed algorithm, and dashed lines show the solution obtained by [3].Figure 3
Comparison of the analytic and approximate solution of Example3 by proposed algorithm. In this figure solid line shows exact solution and dashed lines show approximate solution by proposed algorithm.Figure 4
Comparison between exact and approximate solutions of Example4 for k
=
0. Solid line represents exact solution and dash line represents approximate solution.Figure 5
Comparison between exact and approximate solutions of Example4 for k
=
1
/
2. Solid line represents exact solution and dash line represents approximate solution.Figure 6
Comparison between exact and approximate solutions of Example4 for k
=
2. Solid line represents exact solution and dash line represents approximate solution.The numerical solutions of Examples1, 2, and 3 at different step sizes are shown in Figure 7. From this figure, we see that step sizes have small effect on the numerical solutions of BVPs.Figure 7
Approximate solutions of Examples1, 2, and 3 at different step sizes.
(a)
(b)
(c)


### 6.3. Conclusion

In this paper, we have presented subdivision based iterative collocation algorithm for the solution of nonlinear third-order boundary value problems. The proposed algorithm has been applied on different nonlinear third-order boundary value problems. Numerical results show that the accuracy of approximate solution isO
(
h
3
). We have also observed that the accuracy of the solution can be improved by choosing different subdivision schemes with the proper adjustment of boundary conditions. Our proposed algorithm gives better results comparative to the solution obtained by fourth-degree B-spline [5]. The order of approximation by the proposed algorithm and modified Adomian decomposition method [3] is the same.


## 6.1. Numerical Examples

We find the approximate solutions of the following nonlinear problems to check the accuracy and convergence of subdivision based iterative collocation algorithm.Example 1.
The nonlinear boundary value problem is(85)
y
′
′
′
=
-
2
e
-
3
y
+
4
1
+
x
-
3
,with boundary conditions(86)
y
0
=
0
,
y
′
0
=
1
,
y
1
=
ln
⁡
2
.Example 2.
The nonlinear boundary value problem is(87)
y
′
′
′
x
=
e
-
x
y
2
xsubject to the boundary conditions:(88)
y
0
=
1
,
y
′
0
=
1
,
y
1
=
e
.Example 3.
The nonlinear boundary value problem is(89)
y
′
′
′
x
=
-
e
x
y
2
xsubject to the boundary conditions: (90)
y
0
=
1
,
y
′
0
=
-
1
,
y
1
=
1
e
.The exact solutions of problems (85), (87), and (89) are y
=
ln
⁡
(
1
+
x
), y
=
e
x, and y
=
e
-
x, respectively.Example 4.
We consider the third-order ordinary differential equation:(91)
y
′
′
′
=
y
-
k
,
y
⩾
0
,where k is constant. The initial conditions imposed by Tanner [16] are(92)
y
0
=
1
,
y
′
0
=
0
.The problem is closed by the boundary condition:(93)
y
r
=
0
,and the problem becomes singular at y
=
0. The boundary condition is imposed [17] as(94)
y
r
=
ϵ
,where r is a constant satisfying r
>
0. The analytical solution of (91) by [17](95)
y
x
=
1
6
6
-
7
x
2
+
6
ϵ
x
2
+
x
3
.Case 1 (see [17]).
Fork
=
0 problem (91) takes the form(96)
y
′
′
′
=
y
0
=
1
.
We solve (96) along with conditions (92) and (94) by letting r
=
1; then ϵ
=
11
/
12 and we obtain the results which also support our algorithm. That is, the numerical results have the order of approximation O
(
h
3
). The numerical results are tabulated in Table 4 and graphical representation of these results is shown in Figure 4. These results are obtained after first iteration level. The maximum absolute error is 6.1250
×
1
0
-
3.Case 2 (see [17]).
Fork
=
1
/
2 problem (91) takes the form(97)
y
′
′
′
=
y
-
1
/
2
,with the boundary conditions y
(
0
)
=
1, y
′
0
=
0, and y
(
r
)
=
ϵ. The numerical solution of (97) is obtained by using the proposed numerical algorithm. The solution after third iteration level is presented in Table 5 and their graphical representation is presented in Figure 5. The maximum absolute error for this problem is 6.4004
×
1
0
-
3.Case 3 (see [17, 18]).
In this case, we numerically solve problem (91) for k
=
2; that is,(98)
y
′
′
′
=
y
-
2together with the boundary conditions y
(
0
)
=
1, y
′
0
=
0, and y
(
r
)
=
ϵ. Its exact solution is given in (95). The numerical solution of (98) with given and imposed condition at node points is tabulated in Table 6 and graphical representation is given in Figure 6. The maximum absolute error for this problem is 4.93276
×
1
0
-
3.


## 6.2. Comparison and Discussion

For Examples1, 2, and 3, we use the iterative collocation algorithm described in Section 4, for h
=
1
0
-
1
,
2
0
-
1
,
5
0
-
1 (i.e., for N
=
10,20,50) and t
o
l
=
1
0
-
6 along with eighth-order boundary treatment at end points, to get solutions of nonlinear boundary value problems. The numerical results are obtained after third iteration with condition (53).(i)
The numerical solutions of problems (85), (87), and (89) are presented in Tables 1, 2, and 3, respectively.
(ii)
Caglar et al. [5] solved problem (85) by fourth-degree B-spline algorithm. The maximum absolute errors obtained by the proposed algorithm and by [5] are 4.77
×
1
0
-
3 and 5.80
×
1
0
-
2, respectively. The graphical comparison between exact and approximate solutions is shown in Figure 1. We observe that numerical results obtained by proposed algorithm are better than the results of [5].
(iii)
Hasan [3] solved problem (87) by modified Adomian decomposition method (MADM). We also solve this problem by subdivision based collocation algorithm. Here we observe that the order of approximation by proposed and MADM algorithms is the same (i.e., O
(
h
3
)). The graphical comparison between exact and approximate solutions is shown in Figure 2.
(iv)
The comparison between exact and approximate solutions of problem (89) is given in Figure 3.In Example 4, we consider the problem related to thin film flows. We solve (91) by assuming different values of k and we observer from the numerical results tabulated in Tables 4, 5, and 6 that accuracy of the approximate solution is O
(
h
3
). The numerical results are obtained after first and third iterations for k
=
0 and k
=
1
/
2
,
1 with condition (53), respectively.Table 1
Numerical results of Example1: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
Error by method [5]
0.0
0
0
0
0
0.1
0.0953101798
0.0957100706
0.0003998908
0.009
0.2
0.1823215568
0.1836519852
0.0013304284
0.013
0.3
0.2623642645
0.2648256813
0.0024614168
0.031
0.4
0.3364722366
0.3400108156
0.0035385790
0.045
0.5
0.4054651081
0.4098275640
0.0043624559
0.054
0.6
0.4700036292
0.4747771017
0.0047734725
0.058
0.7
0.5306282511
0.5352702927
0.0046420416
0.058
0.8
0.5877866649
0.5916474874
0.0038608225
0.053
0.9
0.6418538862
0.6441934825
0.0023395963
0.044
1.0
0.6931471806
0.6931471806
0
0.000Table 2
Numerical results of Example2: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
1.105170918
1.1056664532
0.000495535
0.2
1.221402758
1.2232054699
0.001802712
0.3
1.349858808
1.3534912400
0.003632432
0.4
1.491824698
1.4974920669
0.005667369
0.5
1.648721271
1.6562805667
0.007559296
0.6
1.822118800
1.8310443828
0.008925583
0.7
2.013752707
2.0230974700
0.009344763
0.8
2.225540928
2.2338920118
0.008351084
0.9
2.459603111
2.4650310392
0.005427928
1.0
2.718281828
2.7182818285
0Table 3
Numerical results of Example3: h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9048374180
0.9045439989
0.0002934191
0.2
0.8187307531
0.8177120092
0.0010187439
0.3
0.7408182207
0.7388586762
0.0019595445
0.4
0.6703200460
0.6674013997
0.0029186463
0.5
0.6065306597
0.6028146540
0.0037160057
0.6
0.5488116363
0.5446245817
0.0041870544
0.7
0.4965853038
0.4924038346
0.0041814692
0.8
0.4493289641
0.4457666384
0.0035623257
0.9
0.4065696597
0.4043640577
0.0022056020
1.0
0.3678794412
0.3678794412
0Table 4
Numerical results of Example4: k
=
0 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980416667
0.0003750000
0.2
0.9913333333
0.9926666667
0.0013333334
0.3
0.9820000000
0.9846250000
0.0026250000
0.4
0.9706666667
0.9746666667
0.0040000000
0.5
0.9583333333
0.9635416667
0.0052083334
0.6
0.9460000000
0.9520000000
0.0060000000
0.7
0.9346666667
0.9407916667
0.0061250000
0.8
0.9253333333
0.9306666667
0.0053333334
0.9
0.9190000000
0.9223750000
0.0033750000
1.0
0.9166666667
0.9166666667
0Table 5
Numerical results of Example4: k
=
1
/
2 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980508018
0.0003841351
0.2
0.9913333333
0.9927028673
0.0013695340
0.3
0.9820000000
0.9847044600
0.0013695340
0.4
0.9706666667
0.9748013302
0.0041346635
0.5
0.9583333333
0.9637358441
0.0054025108
0.6
0.9460000000
0.9522463844
0.0062463844
0.7
0.9346666667
0.9410670255
0.0064003588
0.8
0.9253333333
0.9309274865
0.0055941532
0.9
0.9190000000
0.9225533666
0.0035533666
1.0
0.9166666667
0.9166666692
0Table 6
Numerical results of Example4: k
=
2 and h
=
1
0
-
1.
x
i
Analytic solutionY
i
Approximate solutionZ
i
Error by proposed method
0.0
1
1
0
0.1
0.9976666667
0.9980024045
0.0003357378
0.2
0.9913333333
0.9925110306
0.0011776973
0.3
0.9820000000
0.9842830603
0.0022830603
0.4
0.9706666667
0.9740862436
0.0034195769
0.5
0.9583333333
0.9627028681
0.0043695348
0.6
0.9460000000
0.9509327630
0.0049327630
0.7
0.9346666667
0.9395952385
0.0049285718
0.8
0.9253333333
0.9295297470
0.0041964137
0.9
0.9190000000
0.9215950588
0.0025950588
1.0
0.9166666667
0.9166666667
0Figure 1
Comparison of the analytic and approximate solution of Example1 by proposed algorithm and [5]. In this figure solid line shows exact solution, dotted lines show approximate solution by proposed algorithm, and dashed lines show the solution obtained by [5].Figure 2
Comparison of the analytic and approximate solution of Example2 by proposed algorithm and [3]. In this figure solid line shows exact solution, dotted lines show approximate solution by proposed algorithm, and dashed lines show the solution obtained by [3].Figure 3
Comparison of the analytic and approximate solution of Example3 by proposed algorithm. In this figure solid line shows exact solution and dashed lines show approximate solution by proposed algorithm.Figure 4
Comparison between exact and approximate solutions of Example4 for k
=
0. Solid line represents exact solution and dash line represents approximate solution.Figure 5
Comparison between exact and approximate solutions of Example4 for k
=
1
/
2. Solid line represents exact solution and dash line represents approximate solution.Figure 6
Comparison between exact and approximate solutions of Example4 for k
=
2. Solid line represents exact solution and dash line represents approximate solution.The numerical solutions of Examples1, 2, and 3 at different step sizes are shown in Figure 7. From this figure, we see that step sizes have small effect on the numerical solutions of BVPs.Figure 7
Approximate solutions of Examples1, 2, and 3 at different step sizes.
(a)
(b)
(c)


## 6.3. Conclusion

In this paper, we have presented subdivision based iterative collocation algorithm for the solution of nonlinear third-order boundary value problems. The proposed algorithm has been applied on different nonlinear third-order boundary value problems. Numerical results show that the accuracy of approximate solution isO
(
h
3
). We have also observed that the accuracy of the solution can be improved by choosing different subdivision schemes with the proper adjustment of boundary conditions. Our proposed algorithm gives better results comparative to the solution obtained by fourth-degree B-spline [5]. The order of approximation by the proposed algorithm and modified Adomian decomposition method [3] is the same.



---

*Source: 5026504-2016-06-09.xml*