# Numerical Solution of Two-Point Boundary Value Problems by Interpolating Subdivision Schemes

Authors: Ghulam Mustafa; Syeda Tehmina Ejaz
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A numerical interpolating algorithm of collocation is formulated,
based on 8-point binary interpolating subdivision schemes for the generation of curves, to solve the two-point third order boundary value problems.
It is observed that the algorithm produces smooth continuous solutions of the problems. Numerical examples are given to illustrate the algorithm
and its convergence. Moreover, the approximation properties of the collocation algorithm have also been discussed.


---

## Body

## 1. Introduction

Subdivision plays an important role in computer aided geometric design. Particularly, subdivision schemes for curve design consist of repeated refinement of control polygons. In particular, the linear schemes are well studied with rules for defining the control points at the finer level as finite linear combinations of the control points in the coarser level.The concept of subdivision has been initiated by de Rham [1]. Later on, in [2], Dyn et al. studied a family of schemes with mask of size four, indexed by a tension parameter. In [3], Dubuc and Deslauriers studied a family of schemes indexed by the size of the mask 2m and the arity (or base) b; that is, at each subdivision step, the scheme (m, b) uses masks of 2m points to compute b new points corresponding to each old edge. Mustafa and Rehman [4] unified all existing even-point interpolating and approximating schemes by offering general formula for the mask of (2m+4)-point even-ary subdivision scheme. Aslam et al. [5] presented an explicit formula which unified the mask of (2m-1)-point interpolating as well as approximating schemes. Mustafa et al. [6] presented an explicit formula for the mask of odd-points b-ary (for any odd m≥3) interpolating subdivision schemes. Higher arity schemes have been introduced in [7, 8]. Mustafa and Ivrissimitzis [9] have showed that when subdivision is used as a data modelling tool large support, or large arity schemes, which generally produce smoother curves, it may outperform simpler schemes.An interpolating subdivision scheme approach to the construction of approximate solutions of two-point second order boundary value problems was first time introduced in [10]. In this approach a method of collocation was formulated for linear two-point second order boundary value problems. It is proved that the algorithms produce smooth continuous solutions provided the algorithms are chosen appropriately. Later on, in [11], they reformulated the collocation method by subdivision scheme in order to compute numerical solutions for two-point boundary-value problems of differential equations with deviating arguments. They demonstrated that their approach has further development to their previous works for solving various types of two-point second order boundary-value problems.They approximated the derivative boundary conditions by forward difference operator so it is difficult (with no flexibility to improve the solution) to generalize their approach to solve third or higher order boundary value problems. We express derivative boundary conditions at end points by using interpolatory subdivision algorithm; therefore it is easy to generalize our approach to deal with higher order boundary value problems.In this paper, we have reformulated the collocation algorithm by using 8-point interpolating subdivision scheme to compute the approximate solution of two-point third order boundary value problems. Particularly, collocation algorithm with septic precision treatments at the endpoints has the power of approximationO(h2). Our reformulated collocation algorithm treats the following type of two-point third order boundary value problem:
(1)y′′′(x)=a(x)y(x)+b(x),0≤x≤1,hly(0)=yr,y′(0)=0,y(1)=yl,
where a(x) and b(x) are continuous and a(x)≥0 on [0,1].The outline of the paper is as follows. In Section2, we rewrite general form of interpolating subdivision scheme for curve design [12] and some related results. The 8-point binary interpolating scheme and derivatives of its basis function have also been discussed in this section. In Section 3, a numerical interpolating algorithm of collocation to solve (1) is formulated and its boundary treatments are discussed. In Section 4, approximation properties of this algorithm are given. In Section 5, numerical examples are presented. Section 6 is devoted for conclusion and the possible future research directions.


## 2. Interpolating Schemes for Curve Design

A general compact form of symmetric univariate binary interpolating subdivision scheme [12] which maps polygon pk={pik}i∈Z to a refined polygon pk+1={pik+1}i∈Z is defined by
(2)p2ik+1=pik,p2i+1k+1=∑j=0n‍Ln,j(pi-jk+pi+j+1k),
where n is called the degree of the scheme and the constants are given by
(3)Ln,j=((2n+1)!!)22(4n)(2n+1)!(-1)j(2j+1)(2n+1n-j),hhhhhhhhhhhhhhhhhhhhlj=0,1,2,…,n,
where (2n+1n-j) denotes the binomial coefficient.The boundary treatments are necessary to produce smooth curve segments by scheme (2). Normally higher order approximation formulae should be used near the ends of the segments and thus Lagrange formulae of order 2n+1 are recommended.Remark 1.
Letϕ(x) be the limit curve generated from the cardinal data {pi=(i,δ0)T}; that is, ϕ(x) is the fundamental solution of the subdivision scheme (2); then
(4)ϕ(i)={1,i=0,0,i≠0.
Furthermore, ϕ(x) satisfies the following two-scale equation:
(5)ϕ(x)=ϕn(x)=ϕ(2x)+∑j=-nn‍Ln,|j|ϕ(2x-2j+1),x∈R.Lemma 2 (see [12, 13]).
The support of the fundamental solutionϕn(x) to scheme (2) is finite. Explicitly, supp⁡ϕn(x)=(-2n-1,2n+1).Lemma 3 (see [10]).
Given a square matrixA of order n, let the normalized left and right (generalized) eigenvectors of A be denoted by {ηi,ξi}. Then for any vector f∈Rn, there exists following Fourier expansion:
(6)f=∑i=1n‍(fTηi)ξi.Lemma 4 (see [10]).
SupposeF(t), t∈R is a regular and C2n+2 curve in Rm, m≥2. Let P(t), t∈R be the limit curve generated by (2) from the initial data Pi=F(ih), i∈Z, 0<h<1. Then, on any finite interval [a,b], the following estimates hold:
(7)∥F(ht)-p(t)∥∞≤M2n+2(F)(2n+2)!h2n+2=O(h2n+2),∥hjFj(ht)-pj(t)∥∞=O(h2n+2-j),hhhhhhhhhhhlhhj=0,1,2,…,n+22,
where the number M2n+2(F) depends only on the derivatives of F(t) and n.


### 2.1. 8-Point Interpolating Scheme

Forn=3 by (2) and (3), we have the following 8-point binary interpolating subdivision scheme for curve design:
(8)p2ik+1=pikp2i+1k+1=12252048(pik-pi+1k)-2452048(pi-1k-pi+2k)+492048(pi-2k-pi+3k)-52048(pi-3k-pi+4k).
This scheme is C3-continuous in [14] and reproduces polynomial curve of degree seven in [15]. The local subdivision matrix of (8) is denoted and defined by(9)S=(0001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000),where L3,0=1225/2048, L3,1=-245/2048, L3,2=49/2048, and L3,3=-5/2048. The some of its eigenvalues is
(10)λ=1,12,14,18,116,132,164,1128.
For an eigenvalue λ, the eigenvectors ξ and η that satisfy Sξ=λξ and ηST=ηλ are called right and left eigenvectors of the matrix S, respectively. Some of the normalized left and right eigenvectors corresponding to eigenvalues are
(11)ξ0=(1,1,1,1,1,1,1,1,1,1,1,1,1)T,η0=(0,0,0,0,0,0,1,0,0,0,0,0,0)T,ξ1=(-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6)T,η1=(15946360)(-5,1024,13225,-199680,1141695,-4715520,0,4715520,-1141695,199680,-13225,-1024,5)T,ξ2=(36,25,16,9,4,1,0,1,4,9,16,25,36)T,η2=(134546860)(275,-28160,-182613,2607616,-12053651,45634048,-71955030,45634048,-12053651,2607616,-182613,-28160,275)T,ξ3=(-216,-125,-64,-27,-8,-1,0,1,8,27,64,125,216)T,η3=(115039360)(225,-11520,10952,476928,-3047987,4677632,0,-4677632,3047987,-476928,-10952,11520,-225)T.
Since ξiTηj=1 for i=j and 0 otherwise then by using Lemmas 2 and 3, we get the following result.Lemma 5.
The fundamental solution (Cardinal basis)ϕ(x) is thrice continuously differentiable and supported on (-7,7) and its derivatives at integers are given by
(12)ϕ′(i)=2sign⁡(i)e|i|Tη1,ϕ′′(i)=22e|i|Tη2,ϕ′′′(i)=23sign⁡(i)e|i|Tη3,0000000i000-6≤i≤6,
where
(13)e0=(0,0,0,0,0,0,1,0,0,0,0,0,0)T,e1=(0,0,0,0,0,1,0,0,0,0,0,0,0)T,e2=(0,0,0,0,1,0,0,0,0,0,0,0,0)T,e3=(0,0,0,1,0,0,0,0,0,0,0,0,0)T,e4=(0,0,1,0,0,0,0,0,0,0,0,0,0)T,e5=(0,1,0,0,0,0,0,0,0,0,0,0,0)T,e6=(1,0,0,0,0,0,0,0,0,0,0,0,0)T,ϕ′(0)=0,ϕ′(±1)=∓7859249553,ϕ′(±2)=±76113198212,ϕ′(±3)=∓332849553,ϕ′(±4)=±2645594636,ϕ′(±5)=±256743295,ϕ′(±6)=1594636,ϕ′′(0)=-34264341124,ϕ′′(±1)=57042561079505,ϕ′′(±2)=-120536518636040,ϕ′′(±3)=3259521079505,ϕ′′(±4)=-608712878680,ϕ′′(±5)=-704215901,ϕ′′(±6)=551727208,ϕ′′′(0)=0,ϕ′′′(±1)=∓292352117495,ϕ′′′(±2)=±30479871879920,ϕ′′′(±3)=∓331213055,ϕ′′′(±4)=∓1369234990,ϕ′′′(±5)=±162611,ϕ′′′(±6)=∓541776.The graphical representations of the basis limit function defined on cardinal data and its derivatives up to order three forn=3 are shown in Figure 1. Figure 1(a) represents the basis limit function defined in (4). Graphical representations of first, second and third derivatives of basis limit functions obtained from (5) for n=3 are shown in Figures 1(b), 1(c) and 1(d) at i=0,1, −1, respectively. The numeric values of first, second and third derivative are given in (13).Interpolatory basis functionϕ3(i) is shown in (a), first derivative of ϕ3(i) i-e. ϕ3′(i) shown in (b), second derivative of ϕ3(i) i-e. ϕ3′′(i) shown in (c), and third derivative of ϕ3(i) i-e. ϕ3′′′(i) shown in (d), respectively.
(a)
(b)
(c)
(d)


## 2.1. 8-Point Interpolating Scheme

Forn=3 by (2) and (3), we have the following 8-point binary interpolating subdivision scheme for curve design:
(8)p2ik+1=pikp2i+1k+1=12252048(pik-pi+1k)-2452048(pi-1k-pi+2k)+492048(pi-2k-pi+3k)-52048(pi-3k-pi+4k).
This scheme is C3-continuous in [14] and reproduces polynomial curve of degree seven in [15]. The local subdivision matrix of (8) is denoted and defined by(9)S=(0001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000000000L3,3L3,2L3,1L3,0L3,0L3,1L3,2L3,30000000001000),where L3,0=1225/2048, L3,1=-245/2048, L3,2=49/2048, and L3,3=-5/2048. The some of its eigenvalues is
(10)λ=1,12,14,18,116,132,164,1128.
For an eigenvalue λ, the eigenvectors ξ and η that satisfy Sξ=λξ and ηST=ηλ are called right and left eigenvectors of the matrix S, respectively. Some of the normalized left and right eigenvectors corresponding to eigenvalues are
(11)ξ0=(1,1,1,1,1,1,1,1,1,1,1,1,1)T,η0=(0,0,0,0,0,0,1,0,0,0,0,0,0)T,ξ1=(-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6)T,η1=(15946360)(-5,1024,13225,-199680,1141695,-4715520,0,4715520,-1141695,199680,-13225,-1024,5)T,ξ2=(36,25,16,9,4,1,0,1,4,9,16,25,36)T,η2=(134546860)(275,-28160,-182613,2607616,-12053651,45634048,-71955030,45634048,-12053651,2607616,-182613,-28160,275)T,ξ3=(-216,-125,-64,-27,-8,-1,0,1,8,27,64,125,216)T,η3=(115039360)(225,-11520,10952,476928,-3047987,4677632,0,-4677632,3047987,-476928,-10952,11520,-225)T.
Since ξiTηj=1 for i=j and 0 otherwise then by using Lemmas 2 and 3, we get the following result.Lemma 5.
The fundamental solution (Cardinal basis)ϕ(x) is thrice continuously differentiable and supported on (-7,7) and its derivatives at integers are given by
(12)ϕ′(i)=2sign⁡(i)e|i|Tη1,ϕ′′(i)=22e|i|Tη2,ϕ′′′(i)=23sign⁡(i)e|i|Tη3,0000000i000-6≤i≤6,
where
(13)e0=(0,0,0,0,0,0,1,0,0,0,0,0,0)T,e1=(0,0,0,0,0,1,0,0,0,0,0,0,0)T,e2=(0,0,0,0,1,0,0,0,0,0,0,0,0)T,e3=(0,0,0,1,0,0,0,0,0,0,0,0,0)T,e4=(0,0,1,0,0,0,0,0,0,0,0,0,0)T,e5=(0,1,0,0,0,0,0,0,0,0,0,0,0)T,e6=(1,0,0,0,0,0,0,0,0,0,0,0,0)T,ϕ′(0)=0,ϕ′(±1)=∓7859249553,ϕ′(±2)=±76113198212,ϕ′(±3)=∓332849553,ϕ′(±4)=±2645594636,ϕ′(±5)=±256743295,ϕ′(±6)=1594636,ϕ′′(0)=-34264341124,ϕ′′(±1)=57042561079505,ϕ′′(±2)=-120536518636040,ϕ′′(±3)=3259521079505,ϕ′′(±4)=-608712878680,ϕ′′(±5)=-704215901,ϕ′′(±6)=551727208,ϕ′′′(0)=0,ϕ′′′(±1)=∓292352117495,ϕ′′′(±2)=±30479871879920,ϕ′′′(±3)=∓331213055,ϕ′′′(±4)=∓1369234990,ϕ′′′(±5)=±162611,ϕ′′′(±6)=∓541776.The graphical representations of the basis limit function defined on cardinal data and its derivatives up to order three forn=3 are shown in Figure 1. Figure 1(a) represents the basis limit function defined in (4). Graphical representations of first, second and third derivatives of basis limit functions obtained from (5) for n=3 are shown in Figures 1(b), 1(c) and 1(d) at i=0,1, −1, respectively. The numeric values of first, second and third derivative are given in (13).Interpolatory basis functionϕ3(i) is shown in (a), first derivative of ϕ3(i) i-e. ϕ3′(i) shown in (b), second derivative of ϕ3(i) i-e. ϕ3′′(i) shown in (c), and third derivative of ϕ3(i) i-e. ϕ3′′′(i) shown in (d), respectively.
(a)
(b)
(c)
(d)


## 3. Numerical Interpolating Collocation Algorithm

In this section, first we formulate a numerical interpolating collocation algorithm for linear third order two-point boundary value problems. Then we settle down the boundary conditions to get unique solution.


### 3.1. The Collocation Algorithm

LetN be a positive integer (N≥6), h=1/N, and xi=i/N=ih, i=0,1,2,…,N, and set ai=a(xi), bi=b(xi). Let
(14)Z(x)=∑i=-6N+6‍ziϕ(x-xih),0≤x≤1
be the approximate solution to (1), where {zi} are the unknown to be determined by (1). The collocation algorithm, together with the boundary conditions to be discussed, is given by setting
(15)Z′′′(xj)=a(xj)Z(xj)+b(xj),j=0,1,2,…,N,
where
(16)Z′(xj)=1h∑i=-6N+6‍ziϕ′(xj-xih),Z′′(xj)=1h2∑i=-6N+6‍ziϕ′′(xj-xih),Z′′′(xj)=1h3∑i=-6N+6‍ziϕ′′′(xj-xih).
Using (14) and (16) in (15), we get following N+1 system of equations:
(17)∑i=-6N+6‍ziϕ′′′(xj-xih)-h3aj∑i=-6N+6‍ziϕ(xj-xih)=h3bj,hhhhhhhhhhhhhhhhhhhhhhhhlhj=0,1,2,…,N.Now we simplify the above system in the following theorems.Theorem 6.
Forj=0 by (17), one gets
(18)z-6ϕ6′′′+z-5ϕ5′′′+z-4ϕ4′′′+z-3ϕ3′′′+z-2ϕ2′′′+z-1ϕ1′′′+z0q0+z1ϕ-1′′′+z2ϕ-2′′′+z3ϕ-3′′′+z4ϕ-4′′′+z5ϕ-5′′′+z6ϕ-6′′′=h3b0,
where ϕj′′′=ϕ′′′(j) and q0=ϕ0′′′-a0h3.Proof.
Substitutingj=0 in (17), we get
(19){z-6ϕ′′′(x0-x-6h)+z-5ϕ′′′(x0-x-5h)+⋯+zN+5ϕ′′′(x0-xN+5h)+zN+6ϕ′′′(x0-xN+6h)}-a0h3{z-6ϕ(x0-x-6h)+z-5ϕ(x0-x-5h)+⋯+zN+5ϕ(x0-xN+5h)+zN+6ϕ(x0-xN+6h)}=h3b0.
For xj=jh, j=0,1,2…,N, this implies
(20)z-6ϕ′′′(6)+z-5ϕ′′′(5)+⋯+zN+5ϕ′′′(-N-5)+zN+6ϕ′′′(-N-6)-a0h3{z-6ϕ(6)+z-5ϕ(5)+⋯+zN+5ϕ(-N-5)+zN+6ϕ(-N-6)}=h3b0.
Since the support of basis function ϕ(x) is (-7,7), ϕ′(x), ϕ′′(x), and ϕ′′′(x) are zero outside the interval (-7,7); also by (4) and (13), we get
(21)z-6ϕ′′′(6)+z-5ϕ′′′(5)+z-4ϕ′′′(4)+z-3ϕ′′′(3)+z-2ϕ′′′(2)+z-1ϕ′′′(1)+z0ϕ′′′(0)+z1ϕ′′′(-1)+z2ϕ′′′(-2)+z3ϕ′′′(-3)+z4ϕ′′′(-4)+z5ϕ′′′(-5)+z6ϕ′′′(-6)-a0h3z0ϕ(0)=h3b0.
If ϕj′′′=ϕ′′′(j), then
(22)z-6ϕ6′′′+z-5ϕ5′′′+z-4ϕ4′′′+z-3ϕ3′′′+z-2ϕ2′′′+z-1ϕ1′′′+z0(ϕ0′′′-a0h3)+z1ϕ-1′′′+z2ϕ-2′′′+z3ϕ-3′′′+z4ϕ-4′′′+z5ϕ-5′′′+z6ϕ-6′′′=h3b0.
For q0=ϕ0′′′-a0h3, we get (18). This completes the proof.Theorem 7.
Forj=1,2,3,…,N the system (17) is equivalent to
(23)z-6ϕj+6′′′+z-5ϕj+5′′′+⋯+z0ϕj′′′+z1(ϕj-1′′′-ajh3ϕj-1)+z2(ϕj-2′′′-ajh3ϕj-2)+⋯+zN-1(ϕj-N+1′′′-ajh3ϕj-N+1)+zN(ϕj-N′′′-ajh3ϕj-N)+zN+1ϕj-N-1′′′+⋯+zN+6ϕj-N-6′′′=h3bj.Proof.
By expanding (17), we get(24)z-6ϕ′′′(xj-x-6h)+z-5ϕ′′′(xj-x-5h)+⋯+zN+5ϕ′′′(xj-xN+5h)+zN+6ϕ′′′(xj-xN+6h)-ajh3{z-6ϕ(xj-x-6h)+z-5ϕ(xj-x-5h)+⋯+zN+5ϕ(xj-xN+5h)+zN+6ϕ(xj-xN+6h)}=h3bj.
For xj=jh, j=1,2,…,N, we get
(25)z-6ϕ′′′(j+6)+z-5ϕ′′′(j+5)+⋯+zN+5ϕ′′′(j-N-5)+zN+6ϕ′′′(j-N-6)-ajh3{z-6ϕ(j+6)+z-5ϕ(j+5)+⋯+zN+5ϕ(j-N-5)+zN+6ϕ(j-N-6)}=h3bj.
Ifϕj′′′=ϕ′′′(j), for j=1,2,…,N, then
(26)z-6(ϕj+6′′′-ajh3ϕj+6)+z-5(ϕj+5′′′-ajh3ϕj+5)+⋯+zN+5(ϕj-N-5′′′-ajh3ϕj-N-5)+zN+6(ϕj-N-6′′′-ajh3ϕj-N-6)=h3bj.
Since ϕ′(x), ϕ′′(x), and ϕ′′′(x) are zero outside the interval (-7,7) then, by (4) and (13), we get (23).From (18) and (23), we get following undetermined system of (N+1) equations with (N+13) unknowns {zi}:
(27)AZ=D,
where the matrices A, Z, and D of orders (N+1)×(N+13), N+13, and N+1, respectively, are given by(28)A=(ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q0ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′ϕ-6′′′⋯0000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q1ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′⋯00000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q2ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′⋯000000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q3ϕ-1′′′ϕ-2′′′ϕ-3′′′⋯000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000000000000⋯ϕ-5′′′ϕ-6′′′00000000000000⋯ϕ-4′′′ϕ-5′′′ϕ-6′′′),Z=(z-6,z-5,z-4,z-3,z-2,…,zN+6)T, and D=(b0h3,b1h3,b2h3,b3h3,…,bNh3)T, where ϕj′′′=ϕ′′′(j) and qj=ϕ0′′′-ajh3.


### 3.2.  Adjustment of Boundary Conditions

The order of the coefficient matrix (28) is (N+1)×(N+13). In order to get unique solution of the system, we need twelve more conditions. Here we consider only two different cases. In coming section we will show that the approximate solution can be improved by adjusting different boundary conditions.Case 1.
If we assumez0′=0(equivalentlyz0=yr=finite) then two conditions can be achieved by using following given boundary conditions i-e:
(29)z0=yr,z0′=0,zN=yl.
Still we need ten more conditions to get stable system. Since subdivision scheme reproduces seven degree polynomials, we define boundary conditions of order eight for solution of (27). For simplicity only the left end points are discussed and the values of right end points zN+1, zN+2, zN+3, zN+4, zN+5 can be treated similarly.
The valuesz-5, z-4, z-3, z-2, z-1 can be determined by the septic polynomial q(x) interpolating at (xi,zi), 0≤i≤7. Precisely, we have
(30)z-i=q(-xi),i=1,2,3,4,5,
where
(31)q(xi)=∑j=18‍(8j)(-1)j+1Z(xi-j).
Since by (14) Z(xi)=zi for i=1,2,3,4,5 then, by replacing xi by -xi, we have
(32)q(-xi)=∑j=18‍(8j)(-1)j+1z-i+j.
Hence the following boundary conditions can be employed at the left end:
(33)∑j=08‍(8j)(-1)jz-i+j=0,i=5,4,3,2,1.
Similarly, for the right end, we can define zi=q(-xi), i=N+1, N+2, N+3, N+4, N+5, and
(34)q(xi)=∑j=18‍(8j)(-1)j+1zi-j.
So we have the following boundary conditions at the right end:
(35)∑j=08‍(8j)(-1)jzi-j=0,i=N+1,N+2,N+3,N+4,N+5.
Finally, we get the following new system of (N+13) linear equations with (N+13) unknowns {zi}, in which N+1 equations are obtained from (18) and (23), two equations from boundary conditions (29) and ten from boundary conditions (33) and (35):
(36)BZ=R,
where the coefficients matrix B=(B0T,AT,B1T)T, A is defined by (28), and B0 and B1are formed by (29), (33), and (35)(37)B0=(01-828-5670-5628-810000⋯00001-828-5670-5628-81000⋯000001-828-5670-5628-8100⋯0000001-828-5670-5628-810⋯00000001-828-5670-5628-81⋯0000000010000000⋯00),where the first five rows of B0 come from (33) and the sixth row comes from (29) at z0=yr. Consider (38)B1=(00⋯0000000100000000⋯1-828-5670-5628-810000000⋯01-828-5670-5628-81000000⋯001-828-5670-5628-8100000⋯0001-828-5670-5628-810000⋯00001-828-5670-5628-810),where first row of B1 comes from (29) at zN=yl remaining rows come from (35) and the matrices Z and R are defined as
(39)Z=(z-6,z-5,…,zN+5,zN+6)T,R=(0,0,0,0,0,yr,DT,yl,0,0,0,0,0)T.Case 2.
In this case we express the given boundary conditionz0′=0 in the following way.
By using (16) we have
(40)Z′(xj)=1h∑i=-6N+6‍ziϕ′(xj-xih).
As we defined earlier xj=jh if we put j=0 we get x0=0; the above equation can be written as
(41)Z′(0)=1h∑i=-6N+6‍ziϕ′(-i).
Since by boundary condition z0′=Z′(0)=0,
(42)∑i=-6N+6ziϕ′(-i)=0.
By using (13) we can express above equation as
(43)1594636z-6+256743295z-5+2645594636z-4-332849553z-3+76113198212z-2-7859249553z-1+7859249553z1-76113198212z2+332849553z3-2645594636z4-256743295z5-1594636z6=0.
Finally, we get a following new system of(N+13) linear equations with (N+13) unknowns {zi}, in which N+1 equations are obtained from (18) and (23), two equations from boundary conditions (29) and ten from boundary conditions (33) for i=5,4,3,2, (35), and (43):
(44)BZ=R,
where the coefficients matrix B=(B0T,AT,B1T)T, A is defined by (28), and B0 and B1 formed by (29), (33), (35), and (43) are defined as(45)B0=(01-828-5670-5628-8001-828-5670-56280001-828-5670-5600001-828-567015946362567432952645594636-33284955376113198212-7859249553078592495537611319821200000010010000⋯00-81000⋯0028-8100⋯00-5628-810⋯00332849553-2645594636-256743295-15946360⋯0000000⋯00),in B first four rows come from (33) for i=5,4,3,2, fifth row comes from (43), and the last row comes from (29).The matrixB1 is same as defined in Case 1 and the matrices Z and R are defined as
(46)Z=(z-6,z-5,…,zN+5,zN+6)T,R=(0,0,0,0,y′(0),yr,DT,yl,0,0,0,0,0)T.
The nonsingularity of the coefficients matrix B has been discussed in next section.


### 3.3. Existence of the Solution

In this section, we discuss the nonsingularity of the coefficients matrixB. We observe that the coefficients matrix B is neither symmetric nor diagonally dominant. However it can be shown that B is a nonsingular. Since B is almost a band matrix with half band width 7, numerical complexity for solving the linear system using Gaussian elimination is about 49(N+9) multiplications. For large N, the matrix is almost symmetric except the first and last six rows and columns due to the boundary conditions. Therefore we first consider the symmetric part of it, that is, square symmetric matrix C of order N+3 defined as(47)C=(ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′ϕ-6′′′⋯000ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′⋯000ϕ3′′′ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′⋯000ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′⋯00000000000⋯ϕ4′′′ϕ3′′′ϕ2′′′00000000⋯ϕ3′′′ϕ2′′′ϕ1′′′).It can be shown that C is always nonsingular for each value of N. However, B is nonsingular for N⩽1000. We have checked the nonsingularity of matrix B by different methods. In first method we observe that the determinants of matrix B increase as N increases and it is not zero for N⩽1000. So B is nonsingular. The determinants of B at some values of N are shown in Table 1. In second method we observe that for N⩽1000, the eigenvalues of matrix B are nonzero so by [16] matrix B is nonsingular. However for N>1000 matrix B may or may not be nonsingular. Therefore we claim that systems (36) and (44) are stable for N⩽1000.Table 1
Determinants of the matrices.
N
C
B
10
-
8667
/
56
1
/
1048870018371741
50
−177183
1
/
4981270309
100
−552709050
1
/
1964492
500
-
5.033491471916955
×
10
36
4.728852755761116
×
10
21
1000
-
4.477989536166907
×
10
71
42069711017699999
×
10
40


## 3.1. The Collocation Algorithm

LetN be a positive integer (N≥6), h=1/N, and xi=i/N=ih, i=0,1,2,…,N, and set ai=a(xi), bi=b(xi). Let
(14)Z(x)=∑i=-6N+6‍ziϕ(x-xih),0≤x≤1
be the approximate solution to (1), where {zi} are the unknown to be determined by (1). The collocation algorithm, together with the boundary conditions to be discussed, is given by setting
(15)Z′′′(xj)=a(xj)Z(xj)+b(xj),j=0,1,2,…,N,
where
(16)Z′(xj)=1h∑i=-6N+6‍ziϕ′(xj-xih),Z′′(xj)=1h2∑i=-6N+6‍ziϕ′′(xj-xih),Z′′′(xj)=1h3∑i=-6N+6‍ziϕ′′′(xj-xih).
Using (14) and (16) in (15), we get following N+1 system of equations:
(17)∑i=-6N+6‍ziϕ′′′(xj-xih)-h3aj∑i=-6N+6‍ziϕ(xj-xih)=h3bj,hhhhhhhhhhhhhhhhhhhhhhhhlhj=0,1,2,…,N.Now we simplify the above system in the following theorems.Theorem 6.
Forj=0 by (17), one gets
(18)z-6ϕ6′′′+z-5ϕ5′′′+z-4ϕ4′′′+z-3ϕ3′′′+z-2ϕ2′′′+z-1ϕ1′′′+z0q0+z1ϕ-1′′′+z2ϕ-2′′′+z3ϕ-3′′′+z4ϕ-4′′′+z5ϕ-5′′′+z6ϕ-6′′′=h3b0,
where ϕj′′′=ϕ′′′(j) and q0=ϕ0′′′-a0h3.Proof.
Substitutingj=0 in (17), we get
(19){z-6ϕ′′′(x0-x-6h)+z-5ϕ′′′(x0-x-5h)+⋯+zN+5ϕ′′′(x0-xN+5h)+zN+6ϕ′′′(x0-xN+6h)}-a0h3{z-6ϕ(x0-x-6h)+z-5ϕ(x0-x-5h)+⋯+zN+5ϕ(x0-xN+5h)+zN+6ϕ(x0-xN+6h)}=h3b0.
For xj=jh, j=0,1,2…,N, this implies
(20)z-6ϕ′′′(6)+z-5ϕ′′′(5)+⋯+zN+5ϕ′′′(-N-5)+zN+6ϕ′′′(-N-6)-a0h3{z-6ϕ(6)+z-5ϕ(5)+⋯+zN+5ϕ(-N-5)+zN+6ϕ(-N-6)}=h3b0.
Since the support of basis function ϕ(x) is (-7,7), ϕ′(x), ϕ′′(x), and ϕ′′′(x) are zero outside the interval (-7,7); also by (4) and (13), we get
(21)z-6ϕ′′′(6)+z-5ϕ′′′(5)+z-4ϕ′′′(4)+z-3ϕ′′′(3)+z-2ϕ′′′(2)+z-1ϕ′′′(1)+z0ϕ′′′(0)+z1ϕ′′′(-1)+z2ϕ′′′(-2)+z3ϕ′′′(-3)+z4ϕ′′′(-4)+z5ϕ′′′(-5)+z6ϕ′′′(-6)-a0h3z0ϕ(0)=h3b0.
If ϕj′′′=ϕ′′′(j), then
(22)z-6ϕ6′′′+z-5ϕ5′′′+z-4ϕ4′′′+z-3ϕ3′′′+z-2ϕ2′′′+z-1ϕ1′′′+z0(ϕ0′′′-a0h3)+z1ϕ-1′′′+z2ϕ-2′′′+z3ϕ-3′′′+z4ϕ-4′′′+z5ϕ-5′′′+z6ϕ-6′′′=h3b0.
For q0=ϕ0′′′-a0h3, we get (18). This completes the proof.Theorem 7.
Forj=1,2,3,…,N the system (17) is equivalent to
(23)z-6ϕj+6′′′+z-5ϕj+5′′′+⋯+z0ϕj′′′+z1(ϕj-1′′′-ajh3ϕj-1)+z2(ϕj-2′′′-ajh3ϕj-2)+⋯+zN-1(ϕj-N+1′′′-ajh3ϕj-N+1)+zN(ϕj-N′′′-ajh3ϕj-N)+zN+1ϕj-N-1′′′+⋯+zN+6ϕj-N-6′′′=h3bj.Proof.
By expanding (17), we get(24)z-6ϕ′′′(xj-x-6h)+z-5ϕ′′′(xj-x-5h)+⋯+zN+5ϕ′′′(xj-xN+5h)+zN+6ϕ′′′(xj-xN+6h)-ajh3{z-6ϕ(xj-x-6h)+z-5ϕ(xj-x-5h)+⋯+zN+5ϕ(xj-xN+5h)+zN+6ϕ(xj-xN+6h)}=h3bj.
For xj=jh, j=1,2,…,N, we get
(25)z-6ϕ′′′(j+6)+z-5ϕ′′′(j+5)+⋯+zN+5ϕ′′′(j-N-5)+zN+6ϕ′′′(j-N-6)-ajh3{z-6ϕ(j+6)+z-5ϕ(j+5)+⋯+zN+5ϕ(j-N-5)+zN+6ϕ(j-N-6)}=h3bj.
Ifϕj′′′=ϕ′′′(j), for j=1,2,…,N, then
(26)z-6(ϕj+6′′′-ajh3ϕj+6)+z-5(ϕj+5′′′-ajh3ϕj+5)+⋯+zN+5(ϕj-N-5′′′-ajh3ϕj-N-5)+zN+6(ϕj-N-6′′′-ajh3ϕj-N-6)=h3bj.
Since ϕ′(x), ϕ′′(x), and ϕ′′′(x) are zero outside the interval (-7,7) then, by (4) and (13), we get (23).From (18) and (23), we get following undetermined system of (N+1) equations with (N+13) unknowns {zi}:
(27)AZ=D,
where the matrices A, Z, and D of orders (N+1)×(N+13), N+13, and N+1, respectively, are given by(28)A=(ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q0ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′ϕ-6′′′⋯0000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q1ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′⋯00000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q2ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′⋯000000ϕ6′′′ϕ5′′′ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′q3ϕ-1′′′ϕ-2′′′ϕ-3′′′⋯000⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯0000000000000⋯ϕ-5′′′ϕ-6′′′00000000000000⋯ϕ-4′′′ϕ-5′′′ϕ-6′′′),Z=(z-6,z-5,z-4,z-3,z-2,…,zN+6)T, and D=(b0h3,b1h3,b2h3,b3h3,…,bNh3)T, where ϕj′′′=ϕ′′′(j) and qj=ϕ0′′′-ajh3.


## 3.2.  Adjustment of Boundary Conditions

The order of the coefficient matrix (28) is (N+1)×(N+13). In order to get unique solution of the system, we need twelve more conditions. Here we consider only two different cases. In coming section we will show that the approximate solution can be improved by adjusting different boundary conditions.Case 1.
If we assumez0′=0(equivalentlyz0=yr=finite) then two conditions can be achieved by using following given boundary conditions i-e:
(29)z0=yr,z0′=0,zN=yl.
Still we need ten more conditions to get stable system. Since subdivision scheme reproduces seven degree polynomials, we define boundary conditions of order eight for solution of (27). For simplicity only the left end points are discussed and the values of right end points zN+1, zN+2, zN+3, zN+4, zN+5 can be treated similarly.
The valuesz-5, z-4, z-3, z-2, z-1 can be determined by the septic polynomial q(x) interpolating at (xi,zi), 0≤i≤7. Precisely, we have
(30)z-i=q(-xi),i=1,2,3,4,5,
where
(31)q(xi)=∑j=18‍(8j)(-1)j+1Z(xi-j).
Since by (14) Z(xi)=zi for i=1,2,3,4,5 then, by replacing xi by -xi, we have
(32)q(-xi)=∑j=18‍(8j)(-1)j+1z-i+j.
Hence the following boundary conditions can be employed at the left end:
(33)∑j=08‍(8j)(-1)jz-i+j=0,i=5,4,3,2,1.
Similarly, for the right end, we can define zi=q(-xi), i=N+1, N+2, N+3, N+4, N+5, and
(34)q(xi)=∑j=18‍(8j)(-1)j+1zi-j.
So we have the following boundary conditions at the right end:
(35)∑j=08‍(8j)(-1)jzi-j=0,i=N+1,N+2,N+3,N+4,N+5.
Finally, we get the following new system of (N+13) linear equations with (N+13) unknowns {zi}, in which N+1 equations are obtained from (18) and (23), two equations from boundary conditions (29) and ten from boundary conditions (33) and (35):
(36)BZ=R,
where the coefficients matrix B=(B0T,AT,B1T)T, A is defined by (28), and B0 and B1are formed by (29), (33), and (35)(37)B0=(01-828-5670-5628-810000⋯00001-828-5670-5628-81000⋯000001-828-5670-5628-8100⋯0000001-828-5670-5628-810⋯00000001-828-5670-5628-81⋯0000000010000000⋯00),where the first five rows of B0 come from (33) and the sixth row comes from (29) at z0=yr. Consider (38)B1=(00⋯0000000100000000⋯1-828-5670-5628-810000000⋯01-828-5670-5628-81000000⋯001-828-5670-5628-8100000⋯0001-828-5670-5628-810000⋯00001-828-5670-5628-810),where first row of B1 comes from (29) at zN=yl remaining rows come from (35) and the matrices Z and R are defined as
(39)Z=(z-6,z-5,…,zN+5,zN+6)T,R=(0,0,0,0,0,yr,DT,yl,0,0,0,0,0)T.Case 2.
In this case we express the given boundary conditionz0′=0 in the following way.
By using (16) we have
(40)Z′(xj)=1h∑i=-6N+6‍ziϕ′(xj-xih).
As we defined earlier xj=jh if we put j=0 we get x0=0; the above equation can be written as
(41)Z′(0)=1h∑i=-6N+6‍ziϕ′(-i).
Since by boundary condition z0′=Z′(0)=0,
(42)∑i=-6N+6ziϕ′(-i)=0.
By using (13) we can express above equation as
(43)1594636z-6+256743295z-5+2645594636z-4-332849553z-3+76113198212z-2-7859249553z-1+7859249553z1-76113198212z2+332849553z3-2645594636z4-256743295z5-1594636z6=0.
Finally, we get a following new system of(N+13) linear equations with (N+13) unknowns {zi}, in which N+1 equations are obtained from (18) and (23), two equations from boundary conditions (29) and ten from boundary conditions (33) for i=5,4,3,2, (35), and (43):
(44)BZ=R,
where the coefficients matrix B=(B0T,AT,B1T)T, A is defined by (28), and B0 and B1 formed by (29), (33), (35), and (43) are defined as(45)B0=(01-828-5670-5628-8001-828-5670-56280001-828-5670-5600001-828-567015946362567432952645594636-33284955376113198212-7859249553078592495537611319821200000010010000⋯00-81000⋯0028-8100⋯00-5628-810⋯00332849553-2645594636-256743295-15946360⋯0000000⋯00),in B first four rows come from (33) for i=5,4,3,2, fifth row comes from (43), and the last row comes from (29).The matrixB1 is same as defined in Case 1 and the matrices Z and R are defined as
(46)Z=(z-6,z-5,…,zN+5,zN+6)T,R=(0,0,0,0,y′(0),yr,DT,yl,0,0,0,0,0)T.
The nonsingularity of the coefficients matrix B has been discussed in next section.


## 3.3. Existence of the Solution

In this section, we discuss the nonsingularity of the coefficients matrixB. We observe that the coefficients matrix B is neither symmetric nor diagonally dominant. However it can be shown that B is a nonsingular. Since B is almost a band matrix with half band width 7, numerical complexity for solving the linear system using Gaussian elimination is about 49(N+9) multiplications. For large N, the matrix is almost symmetric except the first and last six rows and columns due to the boundary conditions. Therefore we first consider the symmetric part of it, that is, square symmetric matrix C of order N+3 defined as(47)C=(ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′ϕ-6′′′⋯000ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′ϕ-5′′′⋯000ϕ3′′′ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′ϕ-4′′′⋯000ϕ4′′′ϕ3′′′ϕ2′′′ϕ1′′′ϕ0′′′ϕ-1′′′ϕ-2′′′ϕ-3′′′⋯00000000000⋯ϕ4′′′ϕ3′′′ϕ2′′′00000000⋯ϕ3′′′ϕ2′′′ϕ1′′′).It can be shown that C is always nonsingular for each value of N. However, B is nonsingular for N⩽1000. We have checked the nonsingularity of matrix B by different methods. In first method we observe that the determinants of matrix B increase as N increases and it is not zero for N⩽1000. So B is nonsingular. The determinants of B at some values of N are shown in Table 1. In second method we observe that for N⩽1000, the eigenvalues of matrix B are nonzero so by [16] matrix B is nonsingular. However for N>1000 matrix B may or may not be nonsingular. Therefore we claim that systems (36) and (44) are stable for N⩽1000.Table 1
Determinants of the matrices.
N
C
B
10
-
8667
/
56
1
/
1048870018371741
50
−177183
1
/
4981270309
100
−552709050
1
/
1964492
500
-
5.033491471916955
×
10
36
4.728852755761116
×
10
21
1000
-
4.477989536166907
×
10
71
42069711017699999
×
10
40


## 4. Error Estimation

In this section, we discuss the approximation properties of the numerical interpolating collocation algorithm. Since the scheme (8) reproduces polynomial curve of degree seven so by Dyn [14] scheme has approximation order eight. So the collocation algorithm (14) with septic precision treatments at the endpoints has the power of approximation O(h2). Here we present our main result for error estimation.Proposition 8.
Suppose the exact solutiony(x)∈C8[0,1] and {zi} are obtained by solving either (36) or (44) with 8th order boundary condition at the end points; then
(48)∥err⁡(x)∥∞=∥Zj-yj∥∞=O(h2-j),j=0,1,2,
where j denotes the order of derivative.Proof.
Since the order of approximation of subdivision scheme (8) is eight so by using (13), we can write for smooth function y(x) and small h as
(49)y′′′(xj)=2315039360h3×{225y(xj-6h)-11520y(xj-5h)+10952y(xj-4h)+476928y(xj-3h)-3047987y(xj-2h)+4677632y(xj-h)-4677632y(xj+h)+3047987y(xj+2h)-476928y(xj+3h)-10952y(xj+4h)+11520y(xj+5h)-225y(xj+6h)}.
This can be written as
(50)y′′′(xj)=2315039360h3×{225yj-6-11520yj-5+10952yj-4+476928yj-3-3047987yj-2+4677632yj-1-4677632yj+1+3047987yj+2-476928yj+3-10952yj+4+11520yj+5-225yj+6}.
Similarly, we have
(51)Z′′′(xj)=2315039360h3×{225zj-6-11520zj-5+10952zj-4+476928zj-3-3047987zj-2+4677632zj-1-4677632zj+1+3047987zj+2-476928zj+3-10952zj+4+11520zj+5-225zj+6}.
If we define error function e(x)=Z(x)-y(x) and error vectors at the nodes by
(52)e(xj)=Z(xj)-y(xj+jh),-6≤j≤N+6,
or equivalently ej=Zj-y(xj+jh), -6≤j≤N+6, then this implies
(53)ej′=Zj′-y′(x+jh),ej′′=Zj′′-y′′(x+jh),ej′′′=Zj′′′-y′′′(x+jh).
By subtracting (50) from (51), we get
(54)Z′′′(xj)-y′′′(xj)=2315039360h3×{225(zj-6-yj-6)-11520(zj-5-yj-5)+10952(zj-4-yj-4)+476928(zj-3-yj-3)-3047987(zj-2-yj-2)+4677632(zj-1-yj-1)-4677632(zj+1-yj+1)+3047987(zj+2-yj+2)-476928(zj+3-yj+3)-10952(zj+4-yj+4)+11520(zj+5-yj+5)-225(zj+6-yj+6)}.
This implies
(55)e′′′(xj)=2315039360h3×{225ej-6-11520ej-5+10952ej-4+476928ej-3-3047987ej-2+4677632ej-1-4677632ej+1+3047987ej+2-476928ej+3-10952ej+4+11520ej+5-225ej+6}.
By Lemma 5, we get the following expression:
(56)ej′′′=1h3{ϕ6′′′ej-6+ϕ5′′′ej-5+ϕ4′′′ej-4+ϕ3′′′ej-3+ϕ2′′′ej-2+ϕ1′′′ej-1+ϕ0′′′ej+ϕ-1′′′ej+1+ϕ-2′′′ej+2+ϕ-3′′′ej+3+ϕ-4′′′ej+4+ϕ-5′′′ej+5+ϕ-6′′′ej+6}+O(h8),
where j=0,1,2,…,N.
By subtracting (1) from (15), we get
(57)Zj′′′-yj′′′=aj(Zi-Yj).
This implies
(58)ej′′′=ajej,0≤i≤N.
Using (56), we get
(59)ϕ6′′′ej-6+ϕ5′′′ej-5+ϕ4′′′ej-4+ϕ3′′′ej-3+ϕ2′′′ej-2+ϕ1′′′ej-1+qjej+ϕ-1′′′ej+1+ϕ-2′′′ej+2+ϕ-3′′′ej+3+ϕ-4′′′ej+4+ϕ-5′′′ej+5+ϕ-6′′′ej+6=0,
where qj=ϕ0′′′-h3aj and j=0,1,2,…,N.
As0≤x≤1 and xj=jh, j=0,1,2,…,N so e0,e1,…,eN are nonzero while e-6,e-5,…,e-1 and eN+1,eN+2,…,eN+6 are zero because they lie outside the interval [0,1]. Let us define these (the left and right end) error values as
(60)ej={max⁡0≤k≤7{|ek|}O(h8),-6≤j≤0,max⁡N-7≤k≤N{|ek|}O(h8),N≤j≤N+6.
Thus system (59) is equivalent to
(61)(B+O(h6))E=0,
where B+O(h6) is the matrix obtained by deleting the first and last six rows and columns of the matrix B, where
(62)E=(e-6,e-5,e-4,…,eN+4,eN+5,eN+6)T.
By using (7)
(63)(B+O(h6))E=O(h8)∥Z(xj)-y(xj)∥=O(h8)∥E∥=O(h8).
Hence, for small h, the coefficients matrix B+O(h6) will be invertible and thus using the standard result from linear algebra, we have
(64)∥E∥≤(∥B-1∥1-O(h6)O(h8))=O(h2).
This completes the proof.The above discussion suggests that the approximations of the solution computed by the method developed in pervious section are second order accurate approximations. This suggestion is supported by the numerical experiments given in the next section


## 5. Numerical Examples and Discussions

In this section, the numerical collocation algorithm based on 8-point interpolating subdivision scheme described in Section3, with the 8th order boundary conditions at the end points, is tested on the two-point third order boundary value problems. Absolute errors in the analytical solutions are also calculated. For the sake of comparisons, we also tabulated the results in this section.Example 1.
Consider the boundary value problem(65)y′′′(x)=y(x)-3ex,0<x<1,
with boundary conditions y′(0)=0, y(1)=0, y(0)=1. The analytical solution of this problem is
(66)y(x)=(1-x)ex.
By using the collocation algorithm for N=10, we get following solution of the above problem: Zj=∑i=-616‍ziϕ(j-i), where the values of {z-6,z-5,…,z5,z16} by using (36) are
(67){24.5275284967525,0.690493373832105,0.811531107670859,0.901997043612006,0.962835576958380,0.995101658947808,1,0.978925359857389,0.933503956247400,0.865636023693632,0.797539554103408,0.771795250759200,0.651392727341382,0.519777983609805,0.380902186129123,0.199271785305208,0,-0.131140299007903,-0.247662275449558,-0.342309352611581,-0.406997945734322,-0.432827137019005,-0.413625213339484}
and by using (44) are
(68){0.8188703235427,0.8701325409827,0.9139995186898,0.9498836119009,0.9768946175522,0.9940000846665,1,0.9934991816499,0.9728768806892,0.9362530908960,0.8814510666459,0.8059555490326,0.7068662002854,0.5808456281783,0.4240613657868,0.2321210447258,0,-0.2780394553965,-0.6085378480019,-0.9989382153934,-1.457696393548,1.993136716176,-2.54764506053343}.
By using two different boundary treatments presented in Section 3, we obtained two different solutions which are presented in Table 2 along with their absolute errors. The graphical representations of the analytic and approximate solutions of the above problem are shown in Figure 2. Figure 2(a) represents the comparison of analytic and approximate solutions obtained by (36) while analytic and approximate solutions obtained by (44) are shown in Figure 2(b). From this table and figure, we observe that the solution obtained by (44) is significantly better than the solution obtained by (36). So our claim; that is, the approximate solution can be improved by adjusting boundary treatment, is justified. The maximum absolute errors in the solutions obtained by (36) and (44) at step size 10 are 9.755×10-2 and 2.328×10-2, respectively.Table 2
Solutions and error estimation of Example1.
x
j
Analytic solutionY
Approximate solutionZ by (36)
Approximate solutionZ by (44)
Absolute error by (36)
Absolute error by (44)
0.0
1
1
0.0000
0
0
0.1
0.9946538262
0.978925359857389
0.9934991878
0.0157284663
0.00115463837
0.2
0.9771222064
0.933503956247400
0.9728768936
0.0436182502
0.00424531280
0.3
0.9449011656
0.865636023693632
0.9362531105
0.0792651419
0.00864805513
0.4
0.8950948188
0.797539554103408
0.8814510903
0.0975552647
0.01364372850
0.5
0.8243606355
0.771795250759200
0.8059555784
0.0525653847
0.0184050571
0.6
0.7288475200
0.651392727341382
0.7068662336
0.0774547927
0.0219812864
0.7
0.6041258121
0.519777983609805
0.5808456679
0.0843478285
0.0232801441
0.8
0.4451081856
0.380902186129123
0.42406141196
0.06420599959
0.0210467736
0.9
0.2459603111
0.199271785305208
0.23212109670
0.0466885258
0.0138392144
1.0
0
0
0.0000
0.000
0.000Comparison between analytic and approximating solutions.
(a)
(b)Example 2.
Consider the following third order boundary values problem:(69)y′′′(x)=xy(x)+(x3-2x2-5x-3)ex,0<x<1,
with boundary conditions y(0)=0=y(1), y′(0)=1. Its exact solution is
(70)y(x)=x(1-x)ex.
By the homogeneous process of the boundary condition, let u(x)=y(x)-x(1-x); then above problem can be transformed into its equivalent form
(71)u′′′(x)=xu(x)+x2(1-x)+(x3-2x2-5x-3)ex,hhhhhhhhhhhh0<x<1,u(0)=u(1)=u′(0)=0.
The solutions of this problem and their absolute errors obtained by two different boundary treatments are shown in Table 3. The graphical representations of the analytic and approximate solutions are shown in Figure 3. Figure 3(a) represents the comparison of analytic and approximate solutions obtained by (36). Figure 3(b) represents the comparison of analytic and approximate solutions obtained by (44). We observe that the solution obtained by (44) has less absolute error than that of the solution obtained by (36). Again this supports our claim.Table 3
Solutions and error estimation of Example2.
x
j
Analytic solutionY
Approximate solutionZ by (36)
Approximate solutionZ by (44)
Absolute error by (36)
Absolute error by (44)
0.0
0
0
0.0000
0.000
0.000
0.1
0.09946538262
0.1142745006
0.09710744670
0.01480911798
0.00235793592
0.2
0.1954244413
0.2168561360
0.1876024920
0.0214316947
0.0078219493
0.3
0.2834703497
0.3045203260
0.2767792681
0.0210499763
0.0066910816
0.4
0.3580379275
0.3732002884
0.3437107723
0.0151623609
0.0143271552
0.5
0.4121803178
0.4177996505
0.4018316295
0.0056193327
0.0103486883
0.6
0.4373085120
0.4319808639
0.4229368864
0.0053276481
0.0143716256
0.7
0.4228880685
0.4079258312
0.4199575084
0.0149622373
0.0029305601
0.8
0.3560865485
0.3360647276
0.3473442367
0.0200218209
0.0087423118
0.9
0.2213642800
0.2047684968
0.2069998108
0.0165957832
0.0143644692
1.0
0
0
0.0000
0.000
0.000Comparison between analytic and approximating solutions.
(a)
(b)Comparison. The maximum absolute errors in the solutions obtained by (36) and (44) at step size 10 are 2.143×10-2 and 1.437×10-2, respectively. Caglar et al. [17] obtained the same maximum absolute errors but at the step size 32 and 50, respectively. Therefore we conclude that our method is more efficient than that of Caglar et al.


## 6. Conclusion and Future Work

In this work, we present an interpolatory symmetric subdivision algorithm for the numerical solution of third order linear problems. Septic polynomials were used for the adjustment of boundary conditions at the end points. We established collocation method and obtained stable system of linear equations which can be solved by any well-known numerical method. The numerical result showed that the adjustment of boundary conditions at the end points influence the accuracy of approximate solution. That is, the accuracy of the solution can be improved by the proper adjustment of boundary conditions. So our algorithm has flexibility to improve the results by adjusting boundary conditions. The automatic selection and adjustment of the boundary conditions which improve the approximation order of the solution is possible future research direction.



---

*Source: 721314-2014-07-17.xml*