1. Introduction
Most problems arising from scientific and engineering applications, especially applications in geodesics, are boundary value problems (BVPs), which are much more difficult to solve than initial value problems (IVPs). Since it is generally difficult to find closed-form solutions for BVPs, many researchers have attempted to develop methods to find approximate or numerical solutions for BVPs. Well-known methods involve the shooting method [
1], finite difference methods [
2,
3,
4] and spectral methods [
5,
6,
7]. In some real-life situations, the shooting method produces numerically sensitive systems of algebraic equations, which must be solved using other numerical methods [
8].
In the present paper, we consider the following system of linear two-point second-order BVPs:
with boundary conditions
where
. In particular,
and
are given functions, and
and
with
are continuous and sufficiently smooth functions on the interval
. Theorems that systematically list the existence and uniqueness of the problem solutions of (
1) and (
2) have been studied in [
9]. In recent times, applications of linear and non-linear systems of two-point boundary value problems can be found in economics, biology, physics and mathematics. For instance, Nikooeinejad et al. [
10] obtained the approximate solution of two-point BVPS for four applications of differential games in economics and management science using a combined numerical algorithm. In biology, the Shortley–Weller scheme has been implemented for a two-point boundary value problem. This numerical scheme was later applied to investigate tumor growth problems in heterogeneous microenvironments [
11]. On the other hand, the application of two-point boundary value problems has been addressed in the problem of calculating rocket trajectories in the atmosphere [
12].
Several researchers have investigated the linear and non-linear systems of second-order boundary value problems and produced various efficient and accurate numerical methods. These methods include the Laplace homotopy analysis [
13,
14], continuous genetic algorithm method [
15], sinc collocation method [
16,
17], He’s homotopy perturbation method [
18], reproducing kernel space method [
19], multistage optimal homotopy asymptotic method [
20], variational iteration method (VIM) [
21] and Chebyshev finite difference method [
22].
Researchers have been interested in the families of B-splines for their potential to approximate the solution of BVPs accurately and efficiently. B-spline methods have several attractive features and flexibility that make them useful in numerical computation to solve BVPs [
23]. For example, the B-spline is the smoothest interpolation function compared to other piecewise polynomial interpolation functions [
24]. Moreover, B-splines have small local support properties. In recent years, the cubic B-spline collocation method captured the attention of some researchers to solve partial differential equations [
25], fractional differential equations [
26], fractional partial differential equations [
27], etc.
This study focuses on finding the solutions of two-point BVPs using the cubic B-spline method. Bickley was the first to explore cubic splines to approximate the solutions of two-point BVPs [
28]. Later, Albasiny and Hoskins enhanced Bickley’s work by solving the two-point BVPs using a tri-diagonal matrix of coefficients [
29]. Since then, several researchers have earned more interest in employing spline functions for solving BVPs [
30,
31,
32,
33]. Caglar et al. in [
34] evaluated the two-point BVPs solutions using the cubic B-spline basis function. Hamid et al. [
35,
36] considered the ECBM and cubic trigonometric B-spline method for the solution of linear two-point BVPs. Apart from that, Heilat and Ismail [
37] used a hybrid cubic B-spline method to evaluate the solutions of non-linear two-point boundary value problems. Recently, a hybrid cubic B-spline method with an optimized parameter was used by Heilat et al. [
38] to solve linear two-point BVPs.
The linear system of second-order BVPs has gained attention from Caglar and Caglar [
39] and Heilat et al. [
40]. They represented the cubic B-spline method (CBM) and ECBM, respectively. The ECBM involved two parameters in boosting the flexibility of the spline curve. Based on the investigation, the ECBM is the best compared to the CBM, He’s homotopy perturbation method [
18], Laplace homotopy analysis method [
13], reproducing kernel [
19] and sinc-collocation method [
17]. In recent years, Zhang and Niu [
41] found the approximate solution of second-order BVPs using a Lobatto-reproducing kernel and declared the method has high precision accuracy in different spaces.
The new symmetric cubic B-spline method (NCBM) was first studied by Lang and Xu in [
42] to solve non-linear second-order BVPs with two dependent variables. The NCBM is a typical CBM, equipped with a new second-order derivative approximation. Then, Iqbal et al. [
43] explored the NCBM for solving several third-order Emden–Flower type equations. A year after that, Wasim et al. [
44] extended the NCBM and proposed the new extended cubic B-spline method (NECBM) for solving the class of second-order singular BVPs. Moreover, the nonlinear third-order Korteweg–de Vries equations were solved by Abbas et al. [
45] using the NCBM to approximate the solutions. Later, Nazir et al. [
46] improved the method to a new quintic B-spline approximation technique as a method to approximate the numerical solution of the Boussinesq equation. Recently, Nazir et al. [
47] implemented the NCBM for the numerical solutions of coupled viscous Burgers equations.
Thus, motivated by all these works, we aim to figure out whether the proposed method, the NCBM, can perform much better in solving the linear system of two-point second-order BVPs. The rest of this paper is as follows. In
Section 2, the typical definition of cubic B-spline basis functions is described. Then,
Section 3 presents the descriptions of the numerical method. The convergence analysis of the method has been proven in
Section 4. The numerical results and their discussion are summarized in
Section 5. Finally, the paper ends in
Section 6 with a brief conclusion.
2. Preliminary Concepts
This section describes the classical cubic B-spline approximation and the new second-order approximation invented by Lang and Xu [
42]. Let the finite interval
, where
is divided into uniform partitions with a mesh point
using a step size
. The typical cubic B-spline basis function is defined as [
34].
where
. The cubic B-spline holds the geometric invariability, convex hull property and symmetry [
48]. For a sufficiently smooth function
and
, there exist a unique third-degree spline
and
that satisfies the prescribed interpolating conditions given by
and
in which
where
and
are unknown real coefficients to be computed. The values of
and the first and second derivatives
and
at mesh point
are tabulated in
Table 1. From (
4), (
5) and
Table 1, the cubic B-spline approximations
,
,
and
can be simplified as follows:
The second derivatives,
and
can be simplified as
and
, respectively. Subsequently, the new approximation for second-order derivatives can be represented as follows [
42,
49]:
and
We note that this NCBM has up to a fifth-order accuracy [
49].
3. Implementation of the Method
In this section, we extended Caglar’s work [
39] by solving the linear system two-point second-order BVP and adopted the new second-order approximation.
Discretizing (
1) at the knot
gives the following expression:
where
. By substituting (
6)–(
11) into (
12) for
,
, we obtain the following equation:
For
,
By substituting (
6)–(
11) into (
13) for
,
, we obtain the following equation:
For
,
Consequently, we have
linear equations involving
unknowns. Thus, we need four additional equations, which can be obtained from the boundary conditions in (
2) below:
Hence, the above system will have the
dimensional matrix form that can be expressed as:
where matrix
A is given by:
and
The four sub-matrices
,
,
and
are represented as follows:
and
For
Since the matrix
A is a banded matrix, the system of linear equations is solved using a generalization of the Thomas algorithm. This method has been proposed in [
50]. Matlab R2018a running on an Intel(R) CORE(TM) i7-1165G7 CPU 1.30 GHz processor, 8.00 GB RAM, was used to execute the numerical computations.