# The Convergence of Geometric Mesh Cubic Spline Finite Difference Scheme for Nonlinear Higher Order Two-Point Boundary Value Problems

Authors: Navnit Jha; R. K. Mohanty; Vinod Chauhan
Journal: International Journal of Computational Mathematics
 (2014)

---

## Abstract

An efficient algorithm for the numerical solution of higher (even) orders two-point nonlinear boundary value problems has been developed. The method is third order accurate and applicable to both singular and nonsingular cases. We have used cubic spline polynomial basis and geometric mesh finite difference technique for the generation of this new scheme. The irreducibility and monotone property of the iteration matrix have been established and the convergence analysis of the proposed method has been discussed. Some numerical experiments have been carried out to demonstrate the computational efficiency in terms of convergence order, maximum absolute errors, and root mean square errors. The numerical results justify the reliability and efficiency of the method in terms of both order and accuracy.


---

## Body

## 1. Introduction

Consider the following nonlinear two-point boundary value problems of order2
m
,
m
∈
Z
+:
(1)
-
u
(
2
m
)
(
x
)
+
g
(
x
,
u
(
x
)
,
u
(
1
)
(
x
)
,
u
(
2
)
(
x
)
,
…
,
u
(
2
m
-
1
)
(
x
)
)
=
0
,
-
∞
<
a
≤
x
≤
b
<
∞
,
u
(
2
l
)
(
a
)
=
a
l
,
u
(
2
l
)
(
b
)
=
b
l
,
l
=
0
(
1
)
m
-
1
,
where g
∈
C
2
m
+
2
[
a
,
b
], g
u
>
0, and a
l, b
l are finite real constants and -
∞
<
a
≤
x
≤
b
<
∞.The higher order two-point boundary value problems play an important role in various areas of mathematical physics and engineering. The mathematical modeling of geological folding of rock layers [1], theory of plates and shell [2], waves on a suspension bridge [3], reaction diffusion equation [4], astrophysical narrow convection layers bounded by stable layers [5], viscoelastic and inelastic flows, deformation of beam and plate deflection theory [6–8], and so forth are some of the modeling problems in mathematical physics.The analytical solution of (1) for the arbitrary choice of g
(
•
) is difficult and thus we attempt to develop an economical computational method. The existence and uniqueness of the solutions of higher order boundary value problems have been discussed by Agarwal and Krishnamoorthy [9], O’Regan [10], Aftabizadeh [11], and Wei [12]. In the past, the approximate solution for the second, fourth, and/or sixth order two-point boundary value problems has been discussed using homotopy analysis by Liang and Jeffrey [13], reproducing kernel space by Yao and Lin [14], spline solution by Siddiqi and Twizell [15], and the Sinc-Galerkin and Sinc-Collocation methods by Rashidinia and Nabati [16]. The monotone iterative technique and quasilinearization method for the higher order ordinary differential equations have been analysed by Koleva and Vulkov [17]. The geometric mesh technique gains its importance from the theory of electrochemical reaction-convection-diffusion problems in one-dimensional space geometry [18]. The formulation of the geometric mesh finite difference approximations for the two-point singular perturbation problems was developed by Jain et al. [19] and Iyengar et al. [20]. Later, the applications of geometric mesh in the context of second order two-point boundary value problems were studied extensively by Mohanty et al. [21] and Kadalbajoo and Kumar [22].There are some higher order approximation methods available for the numerical solution of differential equations. The most important issues associated with higher order schemes are more computing time and more memory spaces. Therefore, the lower order numerical method with very low computing time is worthwhile. This paper aims to develop a third order accurate numerical method based on cubic polynomial spline basis and geometric mesh finite difference approximations for the numerical solution of second order two-point boundary value problems. The method can be easily extended to fourth, sixth, and even higher order problems. The simplicity of the proposed method lies in its three-point discretizations with evaluations at two neighbouring nodes and one central node. The scheme is inherently compact and minor modifications are required for the singular problems. The resulting systems of algebraic equations are solved using block gauss elimination method obtained from the discretizations of linear differential equations. Classical Newton’s method has been applied to the nonlinear coupled difference equations.Section2 discussed how the third order accurate method based on geometric mesh polynomial spline basis and finite difference discretizations for second order differential equations was derived. The convergence analysis has been discussed briefly in Section 3. Further in Section 4, the proposed method has been extended to higher (even) order two-point boundary value problems; however, we restrict our computations up to the sixth order problems. In Section 5, some examples are presented which show singular and nonsingular nature both in linear and nonlinear cases of second, fourth, and sixth order problems.


## 2. Geometric Mesh Cubic Spline Methods

We introduce a finite set of geometric gridsa
=
x
0
<
x
1
<
⋯
<
x
n
<
x
n
+
1
=
b for the solution region [
a
,
b
]. Let h
k
=
x
k
-
x
k
-
1, k
=
1
(
1
)
n
+
1 be the nonuniform step size and let τ
k
=
h
k
+
1
/
h
k
>
0 be the geometric mesh ratio. Let U
k
=
u
(
x
k
) be the exact solution values and let u
k be the approximate solution of u at the mesh x
k, respectively. Let P
k
(
x
) be the cubic interpolating polynomial which interpolates the solution u
(
x
) of the second order differential equation u
(
2
)
(
x
)
=
g
(
x
) at x
k defined as follows:
(2)
P
k
(
x
)
=
α
⌢
k
(
x
-
x
k
)
3
+
β
⌢
k
(
x
-
x
k
)
2
+
γ
⌢
k
(
x
-
x
k
)
+
δ
⌢
k
,
k
=
1
(
1
)
n
,
where the coefficients α
⌢
k, β
⌢
k, γ
⌢
k, and δ
⌢
k are obtained using the following relations:
(3)
P
k
(
2
)
(
x
k
)
=
G
k
,
P
k
(
2
)
(
x
k
±
1
)
=
G
k
±
1
,
P
k
(
x
k
)
=
U
k
.
We obtain via algebraic calculations the following expressions:
(4)
α
⌢
k
=
1
6
h
k
(
G
k
-
G
k
-
1
)
,
β
⌢
k
=
1
2
G
k
,
γ
⌢
k
=
h
k
6
(
2
G
k
+
G
k
-
1
)
-
1
h
k
(
U
k
+
U
k
-
1
)
,
δ
⌢
k
=
U
k
.
Using the continuity of the first order derivatives P
k
-
1
(
1
)
(
x
k
)
=
P
k
(
1
)
(
x
k
), k
=
1
(
1
)
n, we obtain
(5)
-
τ
k
U
k
-
1
+
(
1
+
τ
k
)
U
k
-
U
k
+
1
+
h
k
2
τ
k
6
×
(
τ
k
G
k
+
1
+
2
(
1
+
τ
k
)
G
k
+
G
k
-
1
)
=
0
,
k
=
1
(
1
)
n
.
Now, we define the following finite difference geometric mesh approximations so as to obtain the third order of accuracy for the second order differential equations with significant first order derivative:
(6)
U
~
k
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
(
U
k
+
1
-
(
1
-
τ
k
2
)
U
k
-
τ
k
2
U
k
-
1
)
,
U
~
k
+
1
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
×
(
(
1
+
2
τ
k
)
U
k
+
1
-
(
1
+
τ
k
)
2
U
k
+
τ
k
2
U
k
-
1
)
,
U
~
k
-
1
(
1
)
=
-
1
h
k
τ
k
(
1
+
τ
k
)
×
(
U
k
+
1
-
(
1
+
τ
k
)
2
U
k
+
τ
k
(
2
+
τ
k
)
U
k
-
1
)
.
With the help of series expansion, we obtain
(7)
U
~
k
(
1
)
=
U
k
(
1
)
+
τ
k
h
k
2
6
U
k
(
3
)
+
O
(
h
k
3
)
,
(8)
U
~
k
+
1
(
1
)
=
U
k
+
1
(
1
)
-
τ
k
h
k
2
6
(
1
+
τ
k
)
U
k
(
3
)
+
O
(
h
k
3
)
,
(9)
U
~
k
-
1
(
1
)
=
U
k
-
1
(
1
)
-
h
k
2
6
(
1
+
τ
k
)
U
k
(
3
)
+
O
(
h
k
3
)
.
Define
(10)
G
~
k
=
g
(
x
k
,
U
k
,
U
~
k
(
1
)
)
G
~
k
±
1
=
g
(
x
k
±
1
,
U
k
±
1
,
U
~
k
±
1
(
1
)
)
.
With the help of (7)–(10), we obtain
(11)
G
~
k
=
g
k
+
h
k
2
6
A
k
τ
k
U
k
(
3
)
+
O
(
h
k
3
)
,
(12)
G
~
k
+
1
=
g
k
+
1
-
h
k
2
6
A
k
τ
k
(
1
+
τ
k
)
U
k
(
3
)
+
O
(
h
k
3
)
,
G
~
k
-
1
=
g
k
-
1
-
h
k
2
6
A
k
(
1
+
τ
k
)
U
k
(
3
)
+
O
(
h
k
3
)
,
where A
k
=
(
∂
g
/
∂
u
(
1
)
)
x
=
x
k.LetU
~
~
k
(
1
)
=
U
~
k
(
1
)
+
ζ
h
k
(
G
~
k
+
1
-
G
~
k
-
1
), where ζ is a free parameter to be determined. With the help of (7) and (12), we find
(13)
U
~
~
k
(
1
)
=
U
k
(
1
)
+
h
k
2
6
(
6
ζ
(
1
+
τ
k
)
-
τ
k
)
U
k
(
3
)
+
O
(
h
k
3
)
.
Therefore, if we choose ζ
=
τ
k
/
6
(
1
+
τ
k
), we achieve U
~
~
k
(
1
)
=
U
k
(
1
)
+
O
(
h
k
3
).Now consider the spline relation (2), for (
0
<
θ
<
1
); we can easily obtain the following approximations:
(14)
U
~
~
k
+
θ
=
P
k
(
x
k
+
θ
)
=
θ
h
k
2
6
τ
k
(
θ
-
τ
k
)
×
(
(
θ
+
τ
k
)
G
~
k
+
1
-
(
θ
-
2
τ
k
)
G
~
k
)
-
1
τ
k
(
(
θ
-
τ
k
)
U
k
-
θ
U
k
+
1
)
,
U
~
~
k
+
θ
(
1
)
=
P
k
(
1
)
(
x
k
+
θ
)
=
h
k
6
τ
k
×
(
(
3
θ
2
-
τ
k
2
)
G
~
k
+
1
-
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
G
~
k
)
-
1
τ
k
h
k
(
U
k
-
U
k
+
1
)
,
U
~
~
k
-
θ
=
P
k
(
x
k
-
θ
)
=
θ
h
k
2
6
×
(
θ
-
1
)
(
(
θ
+
1
)
G
~
k
-
1
-
(
θ
-
2
)
G
~
k
)
+
(
1
-
θ
)
U
k
+
θ
U
k
-
1
,
U
~
~
k
-
θ
(
1
)
=
P
k
(
1
)
(
x
k
-
θ
)
=
h
k
6
×
(
(
3
θ
2
-
6
θ
+
2
)
G
~
k
-
(
3
θ
2
-
1
)
G
~
k
-
1
)
+
1
h
k
(
U
k
-
U
k
-
1
)
.
It is easy to verify that
(15)
U
~
~
k
±
θ
=
U
k
±
θ
+
O
(
h
k
4
)
,
U
~
~
k
±
θ
(
1
)
=
U
k
±
θ
(
1
)
+
O
(
h
k
3
)
.
We define
(16)
G
~
~
k
±
θ
=
g
(
x
k
±
θ
,
U
~
~
k
±
θ
,
U
~
~
k
±
θ
(
1
)
)
,
G
~
~
k
=
g
(
x
k
,
U
~
~
k
,
U
~
~
k
(
1
)
)
.
Thus, we obtain
(17)
G
~
~
k
±
θ
=
G
k
±
θ
+
O
(
h
k
3
)
.
G
~
~
k
=
G
k
+
O
(
h
k
3
)
Analogous to spline relation (5), with the help of (17), consider the following relationship:
(18)
-
τ
k
U
k
-
1
+
(
1
+
τ
k
)
U
k
-
U
k
+
1
+
h
k
2
(
π
k
+
1
G
~
~
k
+
1
+
π
k
G
~
~
k
+
π
k
-
1
G
~
~
k
-
1
)
=
h
k
2
2
U
k
(
2
)
(
τ
k
2
+
τ
k
-
2
(
π
k
+
1
+
π
k
+
π
k
-
1
)
)
+
h
k
3
6
U
k
(
3
)
(
τ
k
3
-
τ
k
-
6
θ
(
π
k
+
1
-
π
k
-
1
)
)
+
h
k
4
24
U
k
(
4
)
(
τ
k
4
+
τ
k
-
12
θ
2
(
π
k
+
1
+
π
k
-
1
)
)
+
O
(
h
k
5
)
.
If the coefficients of h
k
ς, ς
=
2,3
,
4 in the relation (18) vanish, then the right hand side of (18) has local truncation error of O
(
h
k
5
) and hence we obtain
(19)
π
k
+
1
=
τ
k
24
θ
2
(
τ
k
3
+
2
θ
(
τ
k
2
-
1
)
+
1
)
,
π
k
=
-
τ
k
12
θ
2
(
τ
k
3
-
6
θ
2
(
τ
k
+
1
)
+
1
)
,
π
k
-
1
=
τ
k
24
θ
2
(
τ
k
3
-
2
θ
(
τ
k
2
-
1
)
+
1
)
.
Finally, we conclude the derivation for the approximate solution of second order differential equation with significant first order derivative
(20)
u
(
2
)
(
x
)
=
g
(
x
,
u
(
x
)
,
u
(
1
)
(
x
)
)
and the corresponding geometric mesh spline approximation is given by
(21)
-
U
k
+
1
+
(
1
+
τ
k
)
U
k
-
τ
k
U
k
-
1
+
h
k
2
(
π
k
+
1
G
~
~
k
+
1
+
π
k
G
~
~
k
+
π
k
-
1
G
~
~
k
-
1
)
=
O
(
h
k
5
)
.
The difference equation (21) has local truncation error of O
(
h
k
5
) for τ
k
≠
1. However, for the uniform mesh (τ
k
=
1), the local truncation errors become O
(
h
k
6
). If g
(
•
) contains any singularity in the given domain [
a
,
b
], then we use finite Taylor’s expansion for the terms containing x and then neglecting O
(
h
k
5
) terms from (21) and incorporating the boundary values U
0
=
u
(
a
)
=
a
0, U
n
+
1
=
u
(
b
)
=
b
0, we obtain a system of recurrence relations with tridiagonal nature. If g
(
•
) is linear, then the resulting matrix system can be solved by gauss elimination method; otherwise we use Newton’s solver for computing the numerical solutions. For the convergence of the difference equations (21), τ
k must be positive [23].


## 3. Convergence Analysis

In this section, we discuss the convergence theory for the geometric mesh cubic spline method (21) discussed for the numerical solution of second order two-point boundary value problems (20). At x
=
x
k, (20) is written as
(22)
U
k
(
2
)
=
g
(
x
k
,
U
k
,
U
k
(
1
)
)
≡
G
k
(
Say
)
.
Then, the geometric mesh cubic spline scheme (21) is given by
(23)
ψ
k
(
U
k
-
1
,
U
k
,
U
k
+
1
)
+
T
k
(
h
k
)
=
0
,
k
=
1
(
1
)
n
,
where
(24)
ψ
k
(
U
k
-
1
,
U
k
,
U
k
+
1
)
=
-
τ
k
U
k
-
1
+
(
1
+
τ
k
)
U
k
-
U
k
+
1
+
h
k
2
τ
k
24
θ
2
×
(
(
τ
k
3
+
2
θ
(
τ
k
2
-
1
)
+
1
)
G
~
~
k
+
θ
h
h
l
+
(
τ
k
3
-
2
θ
(
τ
k
2
-
1
)
+
1
)
G
~
~
k
-
θ
h
h
l
-
2
(
τ
k
3
-
6
θ
2
(
τ
k
+
1
)
+
1
)
G
~
~
k
)
,
(25)
T
k
(
h
k
)
=
O
(
h
k
5
)
.
Scheme (23) in the matrix vector notation is written as
(26)
ψ
(
U
)
+
T
=
0
,
where U
=
[
U
1
,
…
,
U
n
]
T, T
=
[
T
1
,
…
,
T
n
]
T.Our aim is to find an approximationu for U by solving n
×
n system
(27)
ψ
(
u
)
=
0
.
From (26) and (27), we obtain
(28)
ψ
(
u
)
-
ψ
(
U
)
=
T
.
Let ε
k
=
u
k
-
U
k
,
k
=
1
(
1
)
n and ε
=
u
-
U  be the discretization error vector.Let us define(29)
g
~
k
+
η
=
g
(
x
k
+
η
,
u
k
+
η
,
u
~
k
+
η
(
1
)
)
≃
G
~
k
+
η
,
η
=
0
,
±
1
,
g
~
~
k
±
θ
=
g
(
x
k
±
θ
,
u
~
~
k
±
θ
,
u
~
~
k
+
θ
(
1
)
)
≃
G
~
~
k
±
θ
,
0
<
θ
<
1
,
g
~
~
k
=
g
(
x
k
,
u
k
,
u
~
~
k
(
1
)
)
≃
G
~
~
k
,
(30)
E
~
k
+
η
=
g
~
k
+
η
-
G
~
k
+
η
,
η
=
0
,
±
1
,
E
~
~
k
±
θ
=
g
~
~
k
±
θ
-
G
~
~
k
±
θ
,
0
<
θ
<
1
,
E
~
~
k
=
g
~
~
k
-
G
~
~
k
,
(31)
ε
~
k
(
1
)
=
u
~
k
(
1
)
-
U
~
k
(
1
)
,
ε
~
k
±
1
(
1
)
=
u
~
k
±
1
(
1
)
-
U
~
k
±
1
(
1
)
,
ε
~
~
k
±
θ
=
u
~
~
k
±
θ
-
U
~
~
k
±
θ
,
ε
~
~
k
±
θ
(
1
)
=
u
~
~
k
±
θ
(
1
)
-
U
~
~
k
±
θ
(
1
)
,
ε
~
~
k
(
1
)
=
u
~
~
k
(
1
)
-
U
~
~
k
(
1
)
,
where
(32)
ε
~
k
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
(
ε
k
+
1
-
(
1
-
τ
k
2
)
ε
k
-
τ
k
2
ε
k
-
1
)
,
ε
~
k
+
1
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
(
(
1
+
2
τ
k
)
ε
k
+
1
-
(
1
+
τ
k
)
2
ε
k
+
τ
k
2
ε
k
-
1
)
,
ε
~
k
-
1
(
1
)
=
-
1
h
k
τ
k
(
1
+
τ
k
)
(
ε
k
+
1
-
(
1
+
τ
k
)
2
ε
k
+
τ
k
(
τ
k
+
2
)
ε
k
-
1
)
,
ε
~
~
k
+
θ
=
θ
h
k
2
6
τ
k
(
θ
-
τ
k
)
(
(
θ
+
τ
k
)
E
~
k
+
1
-
(
θ
-
2
τ
k
)
E
~
k
)
+
θ
τ
k
(
ε
k
+
1
-
ε
k
)
+
ε
k
,
ε
~
~
k
-
θ
=
θ
h
k
2
6
(
θ
-
1
)
(
(
θ
+
1
)
E
~
k
-
1
-
(
θ
-
2
)
E
~
k
)
+
θ
(
ε
k
-
1
-
ε
k
)
+
ε
k
,
ε
~
~
k
+
θ
(
1
)
=
h
k
6
τ
k
(
(
3
θ
2
-
τ
k
2
)
E
~
k
+
1
-
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
E
~
k
)
+
1
h
k
τ
k
(
ε
k
+
1
-
ε
k
)
,
ε
~
~
k
-
θ
(
1
)
=
h
k
6
(
(
3
θ
2
-
6
θ
+
2
)
E
~
k
-
(
3
θ
2
-
1
)
E
~
k
-
1
)
+
1
h
k
(
ε
k
-
ε
k
-
1
)
,
ε
~
~
k
(
1
)
=
ε
~
k
(
1
)
-
h
k
τ
k
6
(
1
+
τ
k
)
(
E
~
k
+
1
-
E
~
k
-
1
)
.
With the help of mean value theorem, (30) can be expressed as
(33)
E
~
k
+
η
=
α
k
+
η
ε
~
k
+
η
(
1
)
+
β
k
+
η
ε
k
+
η
,
η
=
0
,
±
1
,
E
~
~
k
±
θ
=
α
k
±
θ
ε
~
~
k
±
θ
(
1
)
+
β
k
±
θ
ε
k
±
θ
,
0
<
θ
<
1
,
E
~
~
k
=
α
k
ε
~
~
k
(
1
)
+
β
k
ε
k
,
where α
l
=
(
∂
g
/
∂
u
(
1
)
)
x
l, β
l
=
(
∂
g
/
∂
u
)
x
l
>
0.Now, from (24), we find that
(34)
ψ
k
(
u
k
-
1
,
u
k
,
u
k
+
1
)
-
ψ
k
(
U
k
-
1
,
U
k
,
U
k
+
1
)
=
-
ε
k
+
1
+
(
1
+
τ
k
)
ε
k
-
τ
k
ε
k
-
1
+
h
k
2
τ
k
24
θ
2
×
(
(
τ
k
3
+
2
θ
(
τ
k
2
-
1
)
+
1
)
E
~
~
k
+
θ
hhl
+
(
τ
k
3
-
2
θ
(
τ
k
2
-
1
)
+
1
)
E
~
~
k
-
θ
hhl
-
2
(
τ
k
3
-
6
θ
2
(
τ
k
+
1
)
+
1
)
E
~
~
k
)
.
Equivalently, in the matrix notation, we have
(35)
ψ
(
u
)
-
ψ
(
U
)
=
M
ε
,
where M
=
(
p
k
,
q
k
,
r
k
) is the tridiagonal matrix with the following entries:
(36)
p
k
=
-
τ
k
+
h
k
2
τ
k
(
τ
k
2
+
(
2
θ
-
1
)
(
τ
k
-
1
)
)
144
θ
2
×
(
α
k
+
θ
(
(
2
τ
k
2
-
6
θ
τ
k
+
3
θ
2
)
α
k
+
(
3
θ
2
-
τ
k
2
)
α
k
+
1
)
hhl
(
(
2
τ
k
2
-
6
θ
τ
k
+
3
θ
2
)
α
k
+
(
3
θ
2
-
τ
k
2
)
α
k
+
1
)
+
β
k
+
θ
h
k
(
θ
-
τ
k
)
θ
(
(
θ
-
2
τ
k
)
α
k
+
(
θ
+
τ
k
)
α
k
+
1
)
)
-
h
k
τ
k
(
τ
k
2
+
(
2
θ
+
1
)
(
τ
k
-
1
)
)
144
θ
2
×
(
α
k
-
θ
(
(
3
θ
2
-
6
θ
+
2
)
α
k
τ
k
h
k
-
(
3
θ
2
-
1
)
(
3
θ
2
-
6
θ
+
2
)
(
θ
-
1
)
(
θ
-
2
)
+
6
(
1
+
τ
k
)
)
hhhhhhhl
×
h
k
(
α
k
-
1
(
2
+
τ
k
)
-
β
k
-
1
h
k
(
1
+
τ
k
)
)
hhhhhhhl
(
3
θ
2
-
6
θ
+
2
)
+
6
(
1
+
τ
k
)
)
-
β
k
-
θ
h
k
θ
hhh
×
(
h
k
(
θ
-
1
)
(
θ
-
2
)
α
k
τ
k
-
h
k
(
θ
2
-
1
)
hhhhhh
×
(
α
k
-
1
(
2
+
τ
k
)
-
β
k
-
1
h
k
(
1
+
τ
k
)
)
hhhhhh
(
3
θ
2
-
6
θ
+
2
)
(
θ
-
1
)
(
θ
-
2
)
+
6
(
1
+
τ
k
)
)
)
+
h
k
τ
k
2
(
τ
k
2
-
τ
k
+
1
-
6
θ
2
)
α
k
72
(
1
+
τ
k
)
θ
2
×
(
6
(
1
+
τ
k
)
+
h
k
(
α
k
+
1
τ
k
+
α
k
-
1
(
2
+
τ
k
)
hhhhhhhhhhhhhh
-
β
k
-
1
h
k
(
1
+
τ
k
)
)
)
,
q
k
=
1
+
τ
k
-
h
k
(
τ
k
3
+
2
θ
(
τ
k
2
-
1
)
+
1
)
144
θ
2
τ
k
×
(
α
k
+
θ
(
h
k
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
h
h
h
h
h
h
l
h
+
h
k
(
3
θ
2
-
τ
k
2
)
α
k
+
1
(
1
+
τ
k
)
+
6
τ
k
)
hhl
+
β
k
+
θ
h
k
(
θ
-
τ
k
)
hhl
×
(
h
k
θ
(
θ
-
2
τ
k
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
hhhhh
α
k
+
θ
(
h
k
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
+
h
k
θ
(
θ
+
τ
k
)
α
k
+
1
(
1
+
τ
k
)
+
6
τ
k
)
)
+
h
k
(
τ
k
3
-
2
θ
(
τ
k
2
-
1
)
+
1
)
144
θ
2
×
(
α
k
-
θ
(
h
k
(
3
θ
2
-
6
θ
+
2
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
h
h
h
h
h
h
h
h
-
h
k
(
3
θ
2
-
1
)
α
k
-
1
(
1
+
τ
k
)
+
6
τ
k
)
hhl
-
β
k
-
θ
h
k
(
θ
-
1
)
hhl
×
(
h
k
θ
(
θ
-
2
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
hhhhl
(
h
k
(
3
θ
2
-
6
θ
+
2
)
(
α
k
(
τ
k
-
1
)
+
β
k
h
k
τ
k
)
-
h
k
θ
(
θ
+
1
)
α
k
-
1
(
1
+
τ
k
)
+
6
τ
k
)
)
+
h
k
(
τ
k
3
-
6
θ
2
(
τ
k
+
1
)
+
1
)
72
θ
2
×
(
α
k
(
6
(
τ
k
-
1
)
+
(
α
k
-
1
+
α
k
+
1
)
h
k
τ
k
)
+
6
β
k
h
k
τ
k
)
,
r
k
=
-
1
+
h
k
(
τ
k
3
+
2
θ
(
τ
k
2
-
1
)
+
1
)
144
θ
2
τ
k
(
1
+
τ
k
)
×
(
α
k
+
θ
(
-
(
2
τ
k
2
-
6
θ
τ
k
+
3
θ
2
)
h
k
α
k
+
6
τ
k
(
1
+
τ
k
)
)
)
hhhhhhh
+
h
k
(
3
θ
2
-
τ
k
2
)
hhhhhhh
×
(
α
k
+
1
(
1
+
2
τ
k
)
+
β
k
+
1
h
k
τ
k
(
1
+
τ
k
)
)
hhhhhhh
(
2
τ
k
2
-
6
θ
τ
k
+
3
θ
2
)
+
6
τ
k
(
1
+
τ
k
)
)
+
β
k
+
θ
θ
h
k
hhh
×
(
-
h
k
(
θ
-
τ
k
)
(
θ
-
2
τ
k
)
α
k
+
h
k
(
θ
-
τ
k
)
(
θ
+
τ
k
)
hhhhhh
×
(
α
k
+
1
(
1
+
2
τ
k
)
+
β
k
+
1
h
k
τ
k
(
1
+
τ
k
)
)
hhhhhh
(
-
(
2
τ
k
2
-
6
θ
τ
k
+
3
θ
2
)
h
k
α
k
+
6
τ
k
(
1
+
τ
k
)
)
)
+
h
k
2
(
τ
k
3
-
2
θ
(
τ
k
2
-
1
)
+
1
)
144
θ
2
(
1
+
τ
k
)
×
(
α
k
-
θ
(
(
3
θ
2
-
6
θ
+
2
)
α
k
+
(
3
θ
2
-
1
)
α
k
-
1
)
hhl
(
(
3
θ
2
-
6
θ
+
2
)
α
k
+
(
3
θ
2
-
1
)
α
k
-
1
)
-
β
k
+
θ
h
k
θ
(
θ
-
1
)
(
(
θ
-
2
)
α
k
+
(
θ
+
1
)
α
k
-
1
)
)
-
h
k
(
τ
k
2
-
τ
k
+
1
+
6
θ
2
)
α
k
72
(
1
+
τ
k
)
θ
2
×
(
6
(
1
+
τ
k
)
-
h
k
τ
k
h
h
l
×
(
α
k
+
1
(
1
+
2
τ
k
)
+
h
k
τ
k
(
1
+
τ
k
)
β
k
+
1
+
α
k
-
1
)
)
.
From (28) and (35), we obtain
(37)
M
ε
=
Τ
.
It is easy to see that for sufficiently small h
k, that is, lim
⁡
h
k
→
0
p
k
=
-
τ
k
≠
0 and lim
⁡
h
k
→
0
⁡
r
k
=
-
1
≠
0.Hence the graphG
(
M
) of matrix M is strongly connected and thus M is irreducible [24].Now let(38)
S
k
=
{
q
k
+
r
k
,
k
=
1
p
k
+
q
k
+
r
k
,
k
=
2
(
1
)
n
-
1
p
k
+
q
k
,
k
=
n
.
Let α
=
min
⁡
k
{
α
k
,
α
k
±
1
,
α
k
±
θ
} and β
=
min
⁡
k
{
β
k
,
β
k
±
1
,
β
k
±
θ
}.Then, we obtain(39)
S
k
≥
τ
k
+
τ
k
h
k
24
θ
2
(
1
+
2
θ
-
τ
k
3
-
2
(
θ
-
1
)
2
(
6
θ
2
-
1
)
τ
k
h
h
h
h
h
h
h
h
l
×
τ
k
2
+
2
(
6
θ
2
-
1
)
τ
k
)
α
+
O
(
h
k
2
)
,
k
=
1
,
S
k
≥
h
k
2
τ
k
2
β
(
1
+
τ
k
)
-
h
k
3
τ
k
48
θ
2
(
τ
k
2
-
1
)
×
(
1
+
2
θ
-
8
θ
2
+
τ
k
2
-
(
1
-
2
θ
)
τ
k
)
α
β
+
O
(
h
k
4
)
,
k
=
2
(
1
)
n
-
1
,
S
k
≥
1
+
h
k
24
θ
2
(
1
+
2
θ
-
12
θ
2
-
τ
k
3
h
h
h
h
h
h
h
h
12
θ
2
-
2
(
θ
-
1
)
τ
k
2
-
2
τ
k
)
α
+
O
(
h
k
2
)
,
k
=
n
.
Thus, for sufficiently small values of h
k, we find that
(40)
S
k
≥
τ
k
>
0
,
k
=
1
;
S
k
≥
0
,
k
=
2
(
1
)
n
-
1
;
S
k
≥
1
>
0
,
k
=
n
.
Hence M is monotone [25]; consequently we find that M
-
1 exists and M
-
1
≥
0 [26]. If M
i
j
-
1 denotes the (
i
,
j
)th element of M
-
1 and we define the matrix-vector norm as(41)
∥
M
-
1
∥
=
max
⁡
1
≤
i
≤
n
∑
j
=
1
n
|
M
ij
-
1
|
,
∥
T
∥
=
max
⁡
1
≤
i
≤
n
|
T
i
|
.
Also, from the theory of matrix, we know that
(42)
∑
j
=
1
n
M
i
j
-
1
S
j
=
1
,
i
=
1
(
1
)
n
.
With the help of Taylor’s expansion, we obtain
(43)
M
i
k
-
1
≤
1
S
k
≤
1
τ
k
-
h
k
24
τ
k
θ
2
×
(
1
+
2
θ
-
τ
k
3
-
2
(
θ
-
1
)
τ
k
2
+
2
(
6
θ
2
-
1
)
τ
k
)
α
+
O
(
h
k
2
)
,
k
=
1
∑
j
=
2
n
-
1
M
i
j
-
1
≤
1
min
⁡
2
≤
j
≤
n
-
1
⁡
S
j
≤
2
h
k
-
2
τ
k
(
1
+
τ
k
)
β
+
(
τ
k
-
1
)
(
τ
k
2
-
(
1
-
2
θ
)
τ
k
+
1
+
2
θ
-
8
θ
2
)
α
12
θ
2
β
τ
k
(
1
+
τ
k
)
h
k
+
O
(
h
k
ξ
)
,
ξ
≥
0
,
M
i
k
-
1
≤
1
S
k
≤
1
+
h
k
24
θ
2
×
(
τ
k
3
+
2
(
θ
-
1
)
τ
k
2
+
2
τ
k
+
12
θ
2
-
2
θ
-
1
)
α
+
O
(
h
k
2
)
,
k
=
n
.
Therefore, from (37) and (41), we obtain
(44)
∥
ε
∥
=
∥
M
-
1
∥
·
∥
T
∥
=
max
⁡
1
≤
i
≤
n
(
|
M
i
1
-
1
|
+
∑
j
=
2
n
-
1
|
M
i
j
-
1
|
+
|
M
i
n
-
1
|
)
·
O
(
h
k
5
)
≤
(
1
S
1
+
1
min
⁡
2
≤
j
≤
n
S
j
+
1
S
n
)
·
O
(
h
k
5
)
=
2
h
k
3
τ
k
(
1
+
τ
k
)
β
+
O
(
h
k
4
)
⟹
∥
ε
∥
=
O
(
h
k
3
)
.
This proves the third order convergence of the geometric mesh spline method. The convergence theory may be now formally stated.Theorem 1.
The geometric mesh cubic spline finite difference method given by (21) for the numerical solution of second order two-point boundary value problems (20) with sufficiently small h
k and 0
<
τ
k
≠
1 gives a third order convergent solution.


## 4. Geometric Mesh Spline Algorithms for Higher Order Boundary Value Problem

The geometric mesh spline method, developed in Section2, can be easily extended to higher order (even) boundary value problem (1) and it can be written as a system of second order boundary value problem:
(45)
-
y
l
(
2
)
(
x
)
+
y
l
+
1
(
x
)
=
0
,
l
=
1
(
1
)
m
-
1
,
-
y
m
(
2
)
(
x
)
+
Υ
(
x
,
y
1
(
x
)
,
y
1
(
1
)
(
x
)
,
…
,
y
m
(
x
)
,
y
m
(
1
)
(
x
)
)
=
0
,
a
<
x
<
b
,
y
l
(
a
)
=
a
l
,
y
l
(
b
)
=
b
l
,
l
=
1
(
1
)
m
.
Then, the third order accurate algorithm based on the geometric mesh cubic spline method for the numerical solution of (45) may be written as follows:
(46)
Y
~
l
k
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
×
(
Y
l
k
+
1
-
(
1
-
τ
k
2
)
Y
l
k
-
τ
k
2
Y
l
k
-
1
)
,
l
=
1
(
1
)
m
,
Y
~
l
k
+
1
(
1
)
=
1
h
k
τ
k
(
1
+
τ
k
)
×
(
(
1
+
2
τ
k
)
Y
l
k
+
1
-
(
1
+
τ
k
)
2
Y
l
k
+
τ
k
2
Y
l
k
-
1
)
,
hhhhhhhhhhhhhhhhhhhhhh
l
=
1
(
1
)
m
,
Y
~
l
k
-
1
(
1
)
=
-
1
h
k
τ
k
(
1
+
τ
k
)
×
(
Y
l
k
+
1
-
(
1
+
τ
k
)
2
Y
l
k
+
τ
k
(
τ
k
+
2
)
Y
l
k
-
1
)
,
hhhhhhhhhhhhhhhhhhhhhhl
l
=
1
(
1
)
m
,
G
~
k
+
η
=
g
(
x
k
+
η
,
Y
1
k
+
η
,
Y
1
k
+
η
(
1
)
,
…
,
Y
m
k
+
η
,
Y
m
k
+
η
(
1
)
)
,
η
=
0
,
±
1
,
Y
~
~
l
k
+
θ
=
θ
h
k
2
6
τ
k
(
θ
-
τ
k
)
×
(
(
θ
+
τ
k
)
Y
l
+
1
k
+
1
-
(
θ
-
2
τ
k
)
Y
l
+
1
k
)
+
θ
τ
k
(
Y
l
k
+
1
-
Y
l
k
)
+
Y
l
k
,
l
=
1
(
1
)
m
-
1
,
Y
~
~
m
k
+
θ
=
θ
h
k
2
6
τ
k
(
θ
-
τ
k
)
(
(
θ
+
τ
k
)
G
~
k
+
1
-
(
θ
-
2
τ
k
)
G
~
k
)
+
θ
τ
k
(
Y
m
k
+
1
-
Y
m
k
)
+
Y
m
k
,
Y
~
~
l
k
-
θ
=
θ
h
k
2
6
(
θ
-
1
)
(
(
θ
+
1
)
Y
l
+
1
k
-
1
-
(
θ
-
2
)
Y
l
+
1
k
)
+
θ
(
Y
l
k
-
1
-
Y
l
k
)
+
Y
l
k
,
l
=
1
(
1
)
m
-
1
,
Y
~
~
m
k
-
θ
=
θ
h
k
2
6
(
θ
-
1
)
(
(
θ
+
1
)
G
~
k
-
1
-
(
θ
-
2
)
G
~
k
)
+
θ
(
Y
m
k
-
1
-
Y
m
k
)
+
Y
m
k
,
Y
~
~
l
k
+
θ
(
1
)
=
h
k
6
τ
k
×
(
(
3
θ
2
-
τ
k
2
)
Y
l
+
1
k
+
1
-
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
Y
l
+
1
k
)
-
1
τ
k
h
k
(
Y
l
k
-
Y
l
k
+
1
)
,
l
=
1
(
1
)
m
-
1
,
Y
~
~
m
k
+
θ
(
1
)
=
h
k
6
τ
k
×
(
(
3
θ
2
-
τ
k
2
)
G
~
k
+
1
-
(
3
θ
2
-
6
θ
τ
k
+
2
τ
k
2
)
G
~
k
)
-
1
τ
k
h
k
(
Y
m
k
-
Y
m
k
+
1
)
,
Y
~
~
l
k
-
θ
(
1
)
=
h
k
6
(
(
3
θ
2
-
6
θ
+
2
)
Y
l
+
1
k
-
(
3
θ
2
-
1
)
Y
l
+
1
k
-
1
)
+
1
h
k
(
Y
l
k
-
Y
l
k
-
1
)
,
l
=
1
(
1
)
m
-
1
,
Y
~
~
m
k
-
θ
(
1
)
=
h
k
6
(
(
3
θ
2
-
6
θ
+
2
)
G
~
k
-
(
3
θ
2
-
1
)
G
~
k
-
1
)
+
1
h
k
(
Y
m
k
-
Y
m
k
-
1
)
,
Y
~
~
l
k
(
1
)
=
Y
~
l
k
(
1
)
+
h
k
τ
k
6
(
1
+
τ
k
)
(
Y
l
+
1
k
+
1
-
Y
l
+
1
k
-
1
)
,
hhhhhhhhhlhhhhh
l
=
1
(
1
)
m
-
1
,
Y
~
~
m
k
(
1
)
=
Y
~
m
k
(
1
)
+
h
k
τ
k
6
(
1
+
τ
k
)
(
G
~
k
+
1
-
G
~
k
-
1
)
,
G
~
~
k
±
θ
=
g
(
x
k
±
θ
,
Y
~
~
1
k
±
θ
,
Y
~
~
1
k
±
θ
(
1
)
,
…
,
Y
~
~
m
k
±
θ
,
Y
~
~
m
k
±
θ
(
1
)
)
,
G
~
~
k
=
g
(
x
k
,
Y
1
k
,
Y
~
~
1
k
(
1
)
,
…
,
Y
m
k
,
Y
~
~
m
k
(
1
)
)
,
(47)
-
Y
l
k
+
1
+
(
1
+
τ
k
)
Y
l
k
-
τ
k
Y
l
k
-
1
+
h
k
2
×
(
π
k
+
1
Y
~
~
l
+
1
k
+
1
+
π
k
Y
~
~
l
+
1
k
+
π
k
-
1
Y
~
~
l
+
1
k
-
1
)
=
O
(
h
k
5
)
,
l
=
1
(
1
)
m
-
1
-
Y
m
k
+
1
+
(
1
+
τ
k
)
Y
m
k
-
τ
k
Y
m
k
-
1
+
h
k
2
×
(
π
k
+
1
G
~
~
k
+
θ
+
π
k
G
~
~
k
+
π
k
-
1
G
~
~
k
-
1
)
=
O
(
h
k
5
)
,
k
=
1
(
1
)
n
.
The boundary values are used to obtain values at k
±
1 for k
=
1 and n, respectively. The numerical scheme may be implemented by neglecting O
(
h
k
5
) terms from the system of (47). The resulting difference equations in case of linear boundary value problems give a m
n
×
m
n block tridiagonal system of equations for the unknowns Y
l
k, l
=
1
(
1
)
m, k
=
1
(
1
)
n and can be easily solved using block gauss elimination method.


## 5. Numerical Results

To illustrate the efficiency of the cubic spline geometric mesh finite difference method, we have computed both nonsingular and singular equations for linear and nonlinear two-point boundary value problems of second, fourth, and sixth order. The boundary conditions are obtained using the analytical solution as a test procedure. The numerical accuracy of results is presented using maximum absolute errors (ε
u
(
m
)
(
∞
)), root mean square errors (ε
u
(
m
)
(
2
)), and computational order of convergence (Θ
m) for interpolating mth order derivative of u
(
x
) with the error tolerance being ≤
1
0
-
15:
(48)
ε
u
(
m
)
(
∞
)
=
max
⁡
1
≤
k
≤
n
|
u
k
(
m
)
-
u
(
m
)
(
x
k
)
|
,
ε
u
(
m
)
(
2
)
=
(
1
n
∑
k
=
1
n
|
u
k
(
m
)
-
u
(
m
)
(
x
k
)
|
2
)
1
/
2
,
Θ
m
=
lo
g
2
(
ε
u
(
m
)
(
2
)
|
n
grids
ε
u
(
m
)
(
2
)
|
2
n
grids
)
.
For the simplicity in computation, we have chosen τ
k
=
τ
=
constant, for k
=
1
(
1
)
n, and computed the geometric mesh as follows:
(49)
x
0
=
α
,
h
1
=
{
(
b
-
a
)
(
1
-
τ
)
(
1
-
τ
n
+
1
)
,
τ
<
1
,
(
b
-
a
)
(
τ
-
1
)
(
τ
n
+
1
-
1
)
,
τ
>
1
.
The subsequent mesh spacing is determined by h
k
+
1
=
τ
h
k, k
=
1
(
1
)
n. If the boundary value problems exhibit layer behaviour near the left boundary, the solution value can be captured by choosing τ
>
1. If the layer occurs at the right boundary, we choose τ
<
1. If the layer occurs in the interior region, then mesh in the first half of the interval may be arranged by choosing τ
>
1 and second half of the interval by choosing τ
<
1 [27]. All the numerical computations are performed using long double length arithmetic in C under Linux operating system with 2 GB operational memory.Example 1 (see [28]).
Consider the linear nonsingular convection diffusion problem(50)
u
(
2
)
(
x
)
-
ξ
u
(
1
)
(
x
)
=
0
,
0
<
x
<
1
,
u
(
0
)
=
1
,
u
(
1
)
=
0
.
The analytical solution is u
(
x
)
=
(
e
ξ
-
1
)
-
1
(
e
ξ
-
e
ξ
x
). The numerical results are computed for ξ = 10, 100 and 500 in Tables 1, 2, and 3 for various values of n.Table 1
Error estimates and order of convergence for Example1 at ξ
=
10.
n
τ
=
0.97
,
θ
=
0.70
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
1.84
e
-
05
9.72
e
-
06
—
2.65
e
-
05
1.16
e
-
05
—
40
1.15
e
-
06
6.90
e
-
07
4.0
1.81
e
-
06
7.88
e
-
07
3.9
80
1.40
e
-
07
9.27
e
-
08
3.0
1.18
e
-
07
5.13
e
-
08
3.9Table 2
Error estimates and order of convergence for Example1 at ξ
=
100.
n
τ
=
0.94
,
θ
=
0.85
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
2.04
e
-
02
4.68
e
-
03
—
8.79
e
-
02
1.97
e
-
02
—
40
1.81
e
-
04
5.38
e
-
05
6.0
1.44
e
-
02
2.33
e
-
03
2.6
80
2.76
e
-
06
1.35
e
-
06
6.0
1.26
e
-
03
1.68
e
-
04
3.5Table 3
Error estimates and order of convergence for Example1 at ξ
=
500.
n
τ
=
0.80
,
θ
=
0.47
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
1.51
e
-
03
4.88
e
-
04
—
6.04
e
-
01
1.69
e
-
01
—
40
3.26
e
-
05
1.38
e
-
05
5.5
3.74
e
-
01
6.38
e
-
02
0.7
80
3.09
e
-
05
9.39
e
-
06
0.1
1.47
e
-
01
1.67
e
-
02
1.3Example 2 (see [29]).
Consider the linear nonsingular fourth order stiff two-point boundary value problem:(51)
u
(
4
)
(
x
)
-
(
1
+
λ
)
u
(
2
)
(
x
)
+
λ
u
(
x
)
=
λ
2
x
2
+
1
,
0
<
x
<
1
.
The theoretical solution is given by u
(
x
)
=
1
+
(
x
2
/
2
)
+
sinh
⁡
(
x
). We know that ±
1 and ±
λ are the eigenvalues of this equation and hence the problem is stiff for large values of λ. We have solved the problem for small as well as large values of λ and the behavior of solution is sufficiently smooth for λ
<
1
0
10 in case of both uniform mesh and geometric mesh. We have computed the approximate solution with λ
=
1
0
12. Tables 4-5 present the accuracy of solution and computational order of convergence in case of uniform mesh and geometric mesh.Table 4
Error estimates and order of convergence for Example2 at λ
=
1
0
12
,
τ
=
1
,
θ
=
1.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
1.33
e
-
09
1.33
e
-
09
9.81
e
-
10
9.83
e
-
10
—
—
40
9.20
e
-
11
9.14
e
-
11
6.68
e
-
11
6.64
e
-
11
3.9
3.9
80
6.04
e
-
12
6.14
e
-
12
4.36
e
-
12
4.35
e
-
12
3.9
3.9Table 5
Error estimates for Example2 at λ
=
1
0
12
,
τ
=
1.001
,
θ
=
0.8.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
8.61
e
-
11
8.66
e
-
11
6.23
e
-
11
6.23
e
-
11
—
—
40
4.05
e
-
12
4.05
e
-
12
2.80
e
-
12
2.80
e
-
12
4.4
4.4
80
9.15
e
-
14
1.65
e
-
13
5.16
e
-
14
7.00
e
-
14
5.5
4.6Example 3 (see [30]).
Consider the linear nonsingular sixth order stiff two-point boundary value problem:(52)
u
(
6
)
(
x
)
-
(
1
+
λ
)
u
(
4
)
(
x
)
+
λ
u
(
2
)
(
x
)
-
λ
x
=
0
,
0
<
x
<
1
.
The theoretical solution is given by u
(
x
)
=
1
+
(
x
3
/
6
)
+
sinh
⁡
(
x
). We know that ±
1 and ±
λ are the eigenvalues of this equation and hence the problem is stiff for large values of λ. We have solved the problem for small as well as large values of λ and the behavior of solution is sufficiently smooth in case of both uniform mesh and geometric mesh. We have computed the approximate solution with λ
=
1
0
8. Tables 6-7 present the accuracy of solution and computational order of convergence in case of uniform mesh and geometric mesh.Table 6
Error estimates for Example3 at θ
=
1
,
τ
=
1
,
λ
=
1
0
8.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
2
20
1.33
e
-
09
1.33
e
-
09
1.33
e
-
09
9.81
e
-
10
9.81
e
-
10
9.81
e
-
10
—
40
9.19
e
-
11
9.19
e
-
11
9.19
e
-
11
6.67
e
-
11
6.67
e
-
11
6.67
e
-
11
3.9
80
6.04
e
-
12
6.04
e
-
12
6.04
e
-
12
4.35
e
-
12
4.35
e
-
12
4.35
e
-
12
3.9Table 7
Error estimates for Example3 at θ
=
0.999
,
τ
=
0.998
,
λ
=
1
0
8.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
2
20
1.31
e
-
09
1.31
e
-
09
6.79
e
-
09
9.71
e
-
10
9.71
e
-
10
4.83
e
-
09
—
40
8.40
e
-
11
8.40
e
-
11
9.87
e
-
10
6.17
e
-
11
6.17
e
-
11
6.67
e
-
10
4.0
80
3.81
e
-
12
3.81
e
-
12
1.45
e
-
10
2.78
e
-
12
2.78
e
-
12
8.27
e
-
11
4.5Example 4 (see [21]).
Consider the linear singular second order problem:(53)
∇
2
u
(
x
)
≡
(
d
2
d
x
2
+
λ
x
d
d
x
-
λ
x
2
)
u
(
x
)
=
f
(
x
)
,
0
<
x
<
1
.
The above problem for λ
=
0,1
,
2 represents Cartesian, cylindrical, and spherical symmetry, respectively. The unknown function f
(
x
) is chosen in such a manner that u
(
x
)
=
e
x
4 is the theoretical solution. The accuracy of numerical and theoretical solutions is obtained for λ
=
1,2 with different values of n. The results are reported with geometric mesh as well as with uniform mesh. Tables 8-9 present the maximum absolute errors, root mean square errors, and computational order of convergence computed for various values of n.Table 8
Error estimates and order of convergence for Example4 at λ
=
1.
n
τ
=
0.991
,
θ
=
0.42
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
3.72
e
-
05
2.65
e
-
05
—
7.20
e
-
05
5.17
e
-
05
—
40
1.35
e
-
06
9.77
e
-
07
4.8
5.03
e
-
06
3.58
e
-
06
3.8
80
1.26
e
-
08
8.52
e
-
09
6.7
3.32
e
-
07
2.34
e
-
07
3.9Table 9
Error estimates and order of convergence for Example4 at λ
=
2.
n
τ
=
0.981
,
θ
=
0.93
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
3.18
e
-
05
2.52
e
-
05
—
8.51
e
-
05
6.10
e
-
05
—
40
1.43
e
-
06
7.63
e
-
07
4.5
5.95
e
-
06
4.21
e
-
06
3.8
80
3.14
e
-
07
8.84
e
-
08
2.2
3.92
e
-
07
2.76
e
-
07
3.9Example 5 (see [31]).
Consider the linear singular fourth order problem:(54)
∇
4
u
(
x
)
≡
(
d
2
d
x
2
+
λ
x
d
d
x
)
2
u
(
x
)
=
(
1
+
2
λ
x
+
λ
(
λ
-
2
)
x
2
-
λ
(
λ
-
2
)
x
3
)
e
x
,
hhhhhhhhhhhhhhhhhh
0
<
x
<
1
.
For =
0,1
,
2, the above problem represents Cartesian, cylindrical, and spherical symmetry, respectively. The theoretical solution is given by u
(
x
)
=
e
x. The accuracy of numerical and theoretical solutions is obtained for λ = 2 with different values of n. The results are reported with geometric mesh as well as with uniform mesh. Tables 10-11 present the maximum absolute errors, root mean square errors, and computational order of convergence. The improvement in the accuracy of solutions is evident with geometric mesh as compared to spline method with uniform mesh.Table 10
Error estimates for Example5 at λ
=
2
,
θ
=
1
,
τ
=
1.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
1.26
e
-
07
7.22
e
-
06
9.33
e
-
08
2.06
e
-
06
—
—
40
8.83
e
-
09
9.17
e
-
07
6.48
e
-
09
1.88
e
-
07
3.8
3.0
80
5.86
e
-
10
1.14
e
-
07
4.28
e
-
10
1.66
e
-
08
3.8
3.0Table 11
Error estimates for Example5 at λ
=
2
,
θ
=
0.92
,
τ
=
0.991.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
7.25
e
-
08
4.08
e
-
06
5.27
e
-
08
1.16
e
-
06
—
—
40
1.48
e
-
09
1.67
e
-
07
9.43
e
-
10
3.90
e
-
08
5.6
4.6
80
5.86
e
-
10
1.22
e
-
07
4.15
e
-
10
1.59
e
-
08
1.3
0.5Example 6 (see [32]).
Consider the sixth order linear singular problem:(55)
∇
6
u
(
x
)
≡
(
d
2
d
x
2
+
λ
x
d
d
x
)
3
u
(
x
)
=
(
1
+
3
λ
(
λ
-
2
)
x
2
+
3
λ
(
λ
-
2
)
(
λ
-
4
)
x
4
)
sinh
⁡
(
x
)
+
λ
x
(
3
+
(
λ
-
2
)
(
λ
-
7
)
x
2
h
h
h
h
l
h
+
3
(
λ
-
2
)
(
λ
-
4
)
x
4
)
cosh
⁡
(
x
)
,
0
<
x
<
1
.
The theoretical solution is given by u
(
x
)
=
sinh
⁡
(
x
). The problem has been solved numerically in the spherical polar coordinate (λ
=
2, x
→
r). The errors estimates for the various values of n are reported in Tables 12-13 for uniform mesh (
τ
=
1
) as well as for geometric mesh (
τ
≠
1
).Table 12
Error estimates for Example6 at θ
=
1
,
τ
=
1
,
λ
=
2.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
2
20
5.89
e
-
09
7.24
e
-
08
4.90
e
-
06
4.24
e
-
09
5.38
e
-
08
1.34
e
-
06
—
40
4.15
e
-
10
5.08
e
-
09
6.55
e
-
07
2.95
e
-
10
3.73
e
-
09
1.28
e
-
07
3.8
80
2.76
e
-
11
3.37
e
-
10
8.45
e
-
08
1.95
e
-
11
2.46
e
-
10
1.18
e
-
08
3.9Table 13
Error estimates for Example6 at θ
=
0.28
,
τ
=
0.999
,
λ
=
2.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
2
20
1.14
e
-
09
8.16
e
-
08
3.59
e
-
05
6.91
e
-
10
5.46
e
-
08
9.64
e
-
06
—
40
9.93
e
-
11
4.64
e
-
09
4.99
e
-
06
6.34
e
-
11
3.01
e
-
09
9.67
e
-
07
3.5
80
9.28
e
-
12
1.88
e
-
10
6.88
e
-
07
6.42
e
-
12
1.27
e
-
10
9.54
e
-
08
3.4Example 7 (see [33]).
Consider the nonlinear nonsingular system that describes heat and mass transfer and an exothermal chemical reaction in a flat plate:(56)
u
(
2
)
(
x
)
+
4
α
(
1
-
u
(
x
)
)
exp
⁡
(
β
γ
u
(
x
)
1
+
β
u
(
x
)
)
=
f
(
x
)
,
0
<
x
<
1
,
u
(
0
)
=
u
(
1
)
=
0
.
The test solution is taken as u
(
x
)
=
sin
(
π
x
) and it has been observed that the behavior of the solution with uniform mesh deteriorates; however, the geometric mesh shows smooth behaviour of the solution. The accuracy of the solutions is computed for α
=
10, β
=
30, and γ
=
2 in Table 14.Table 14
Error estimates and order of convergence for Example7.
n
τ
=
0.98
,
θ
=
0.95
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
2.77
e
-
07
1.80
e
-
07
—
Overflow
Overflow
—
40
3.53
e
-
08
1.95
e
-
08
3.0
Overflow
Overflow
—
80
8.96
e
-
09
3.59
e
-
09
2.0
Overflow
Overflow
—Example 8 (see [17]).
Consider the nonlinear nonsingular boundary value problems related to Fisher-Kolmogorov and Swift-Hohenberg equation:(57)
u
(
4
)
(
x
)
=
λ
u
(
2
)
(
x
)
+
x
2
u
(
x
)
-
2
(
x
+
1
)
(
u
(
x
)
)
3
+
f
(
x
)
,
0
<
x
<
M
.
The analytical solution is u
(
x
)
=
sin
(
π
x
/
5
). The errors estimates are reported in Tables 15-16 for M
=
2, λ
=
2 and for various values of n using uniform mesh and geometric mesh, respectively.Table 15
Error estimates for Example8 at θ
=
1
,
τ
=
1.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
4.89
e
-
09
3.51
e
-
09
3.61
e
-
09
2.32
e
-
09
—
—
40
3.36
e
-
10
2.41
e
-
10
2.45
e
-
10
1.57
e
-
10
3.9
3.9
80
2.21
e
-
11
1.59
e
-
11
1.60
e
-
11
1.03
e
-
11
3.9
3.9Table 16
Error estimates for Example8 at θ
=
0.92
,
τ
=
0.99.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
8.83
e
-
09
1.20
e
-
09
6.57
e
-
09
8.36
e
-
10
—
—
40
9.27
e
-
10
1.72
e
-
10
6.69
e
-
10
9.99
e
-
11
3.3
2.8
80
1.26
e
-
10
4.63
e
-
11
8.54
e
-
11
2.89
e
-
11
2.9
2.0Example 9 (see [34]).
Consider the sixth order nonlinear nonsingular problem:(58)
u
(
6
)
(
x
)
+
(
u
(
x
)
)
2
=
f
(
x
)
,
0
<
x
<
1
.
The theoretical solution is given by u
(
x
)
=
(
1
-
x
2
)
sin
(
x
). The accuracy of solutions has been reported in Tables 17-18 for various values of n with uniform mesh and geometric mesh, respectively.Table 17
Error estimates for Example9 at θ
=
1
,
τ
=
1.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
0
20
6.32
e
-
08
1.04
e
-
07
1.35
e
-
07
4.64
e
-
08
7.65
e
-
08
9.90
e
-
08
—
40
4.35
e
-
09
7.17
e
-
09
9.27
e
-
09
3.16
e
-
09
5.20
e
-
09
6.73
e
-
09
3.9
80
2.86
e
-
10
4.70
e
-
10
6.08
e
-
10
2.06
e
-
10
3.39
e
-
10
4.39
e
-
10
3.9Table 18
Error estimates for Example9 at θ
=
0.99996
,
τ
=
0.99995.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
4
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
ε
u
(
4
)
(
2
)
Θ
0
20
5.84
e
-
08
9.49
e
-
08
1.32
e
-
07
4.30
e
-
08
6.98
e
-
08
9.74
e
-
08
—
40
3.02
e
-
09
4.57
e
-
09
8.52
e
-
09
2.19
e
-
09
3.32
e
-
09
6.21
e
-
09
4.3
80
6.71
e
-
11
2.19
e
-
10
4.02
e
-
10
4.79
e
-
11
1.57
e
-
10
2.95
e
-
10
5.5Example 10 (see [28]).
Consider nonlinear singular Burger’s equations:(59)
R
e
-
1
∇
2
u
(
x
)
≡
R
e
-
1
(
d
2
d
x
2
+
λ
x
d
d
x
-
λ
x
2
)
u
(
x
)
=
u
(
x
)
u
(
1
)
(
x
)
+
f
(
x
)
,
0
<
x
<
1
.
The analytical solution is given by u
(
x
)
=
x
2
cosh
⁡
(
x
). The errors estimates for various values of n are reported in Tables 19-20 for uniform mesh (
τ
=
1
) as well as for geometric mesh (
τ
≠
1
) with high Reynold’s number R
e
=
1000.Table 19
Error estimates and order of convergence for Example10 at λ
=
1
,
R
e
=
1000.
n
τ
=
0.989
,
θ
=
0.91
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
2.16
e
-
04
1.22
e
-
04
—
3.17
e
-
04
1.69
e
-
04
—
40
8.34
e
-
06
4.66
e
-
06
4.7
3.73
e
-
05
1.67
e
-
05
3.1
80
9.36
e
-
07
4.45
e
-
07
3.2
3.39
e
-
06
1.30
e
-
06
3.5Table 20
Error estimates and order of convergence for Example10 at λ
=
2
,
R
e
=
1000.
n
τ
=
0.992
,
θ
=
0.92
τ
=
1.0
,
θ
=
1.0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
ε
u
(
∞
)
ε
u
(
2
)
Θ
0
20
2.50
e
-
04
1.40
e
-
04
—
3.17
e
-
04
1.68
e
-
04
—
40
1.59
e
-
05
8.06
e
-
06
4.0
3.72
e
-
05
1.67
e
-
05
3.1
80
3.70
e
-
07
8.75
e
-
08
5.4
3.38
e
-
06
1.30
e
-
06
3.5Example 11 (see [35]).
Consider the fourth order nonlinear singular problem which arises in the study of spreading of oil on water:(60)
(
1
-
x
2
)
u
(
4
)
(
x
)
=
4
x
u
(
3
)
(
x
)
-
3
u
(
1
)
(
x
)
u
(
2
)
(
x
)
-
u
(
x
)
u
(
3
)
(
x
)
+
f
(
x
)
,
0
<
x
<
1
.
The theoretical solution is given by u
(
x
)
=
e
x. The errors estimates are reported in Tables 21-22 for various values of n with constant mesh and geometric mesh, respectively. The computational order of convergence has very good agreement with the theoretical order of convergence.Table 21
Error estimates for Example11 at θ
=
1
,
τ
=
1.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
5.89
e
-
06
6.08
e
-
05
4.28
e
-
06
4.34
e
-
05
—
—
40
7.77
e
-
07
8.02
e
-
06
5.57
e
-
07
5.65
e
-
06
3.0
3.0
80
9.91
e
-
08
1.02
e
-
06
7.06
e
-
08
7.18
e
-
07
3.0
3.0Table 22
Error estimates for Example11 at θ
=
0.2
,
τ
=
0.998.
n
ε
u
(
∞
)
ε
u
(
2
)
(
∞
)
ε
u
(
2
)
ε
u
(
2
)
(
2
)
Θ
0
Θ
2
20
2.08
e
-
06
2.68
e
-
05
1.49
e
-
06
1.96
e
-
05
—
—
40
8.58
e
-
08
1.20
e
-
06
6.04
e
-
08
8.80
e
-
07
4.6
4.5
80
5.19
e
-
10
4.33
e
-
08
3.92
e
-
10
1.16
e
-
08
7.4
4.8


## 6. Conclusion

We have developed anO
(
h
k
3
)-accurate algorithm based on cubic spline finite difference approximations using geometric mesh for the numerical treatment of two-point boundary value problems of even higher orders. In this particular case, if the geometric mesh ratio parameter τ is unity, the method achieves fourth order of accuracy. The proposed method shows its unique importance in the equation governing heat and mass transfer and exothermal chemical reaction, where the classical uniform mesh fails to converge, whereas with small rearrangement of grids it shows superiority in terms of order and accuracy. The numerical results confirm that the proposed geometric mesh cubic spline finite difference scheme converges uniformly and is applicable to both singular and nonsingular differential equations. The optimum mesh ratio parameter near unity may be obtained by simulations. Applications of the proposed scheme to the parabolic problems are an open problem.



---

*Source: 527924-2014-07-23.xml*