# An Exponential Spline Difference Scheme for Solving a Class of Boundary Value Problems of Second-Order Ordinary Differential Equations

Authors: Dunqian Cao
Journal: Discrete Dynamics in Nature and Society
 (2020)

---

## Abstract

In this paper, we mainly study an exponential spline function space, construct a basis with local supports, and present the relationship between the function value and the first and the second derivative at the nodes. Using these relations, we construct an exponential spline-based difference scheme for solving a class of boundary value problems of second-order ordinary differential equations (ODEs) and analyze the error and the convergence of this method. The results show that the algorithm is high accurate and conditionally convergent, and an accuracy of1/240h6 was achieved with smooth functions.


---

## Body

## 1. Introduction

In physics, chemistry, biology, sociology, and many other disciplines, there are tremendous problems that can be described by differential equations (DEs), but it is difficult to get their explicit expressions. So, people began to seek the numerical solutions of these problems, which can also be applied to scientific research and engineering practice if their accuracy satisfies the needs. Especially the advent of computers makes it possible to quickly carry out a large number of calculations, which also makes the numerical solution method of DEs become one of the most important branches of computational mathematics. Due to its high smoothness, low power, and easy calculation, the spline function has been widely used in computer graphics, data interpolation and fitting, shape control, and numerical solutions of DEs. There are two main schemes in numerical solutions of DEs using spline functions: the spline finite element method and the spline difference method. The first has a wide range of application and can be applied to many types of equations, but it requires a large amount of calculation. While the second is simple process with a small amount of calculation and high accuracy, but it can only be applied to specific types of equations.In this paper, we mainly focus a class of second-order ordinary differential equations (ODEs):(1)u″+qxu′+pxu=gx,x∈a,b,which meets one of the following boundary conditions.(1)
First boundary condition:(2)ua=μ0,ub=μ1.(2)
Second boundary condition:(3)u1a=μ0,u1b=μ1.(3)
Third boundary condition:(4)ua=μ0,u1b=μ1.(4)
Fourth boundary condition:(5)u1a=μ0,ub=μ1.where μ0 and μ1 are constants and px, qx, and gx are continuous functions in the interval of a,b.Many scholars have been studying such two-point boundary value problems. Albasiny and Hoskins [1] and Raghavarao et al. [2] used cubic polynomial splines to solve such problems. Blue [3] used quintic polysplines to solve such problems. Caglar et al. [4] used cubic B-splines in their solving scheme. Chawla and Shiva Kumar[5] extended the problem to semi-infinite regions.In recent years, with the deepening of research, people have begun to use nonpolynomial splines to solve such problems. Zahra [6], Rao and Kumar [7], Tirmizi et al. [8], Ramadan et al. [9], Surla and Stojanović [10], Jha [11, 12], and Kadalbajoo and Patidar [13] have carried out a lot of research in this area and achieved very high computational accuracy.However, there are still many theoretical problems to be broken in the study of nonpolynomial splines. Due to the diversity of nonpolynomial splines, it is crucial to choose the basis and parameters in solving the problem. However, there is still no reference in this regard. In this paper, a selected set of spline basis functions was used to deduce the relationship between the derivative and the function value and then to obtain the second-order difference scheme for solving second-order ODEs, which provides a method for solving such problems.


## 2. Exponential Spline Function Space

Exponential spline refers to a type of spline in which the nonpolynomial factors of spline basis functions contain only exponential functions. The exponential spline in this sense is not very specific; it can contain many forms of exponential spline, which can produce substantially different splines, and is inconvenient to study. Therefore, the exponential spline refers to that with a specific form in the rest of this work.Next, we define an exponential spline function space. Let(6)six=ai+bix−xi+cieτ1x−xi+dieτ2x−xi,i=1,2,⋯,n,where ai, bi, ci, and di are coefficients and τ1 and τ2 are parameters with τ1≠τ2.Definition 1.
The following function space(7)ℰ3rΔn=s:sx=six,x∈Ii,i=1,2,⋯,nsx∈Cra,b,r<3,is called the cubic r-order exponential spline function space.
Obviously, the functionsx in ℰ32Δn must meet(8)sirxi=si+1rxi,r=0,1,2;i=1,2,⋯,n−1.
The dimension ofℰ32Δn is n+3. Then, we find a set of basic functions with local supports for ℰ32Δn.
Assumesx∈ℰ32Δn, for given j, 2≤j≤n−2; let(9)sj−1rxj−2=0,r=0,1,2,sj+2rxj+2=0,r=0,1,2,sj+krxj+k=sj+k+1rxj+k,k=−1,0,1,r=0,1,2,sjxj=1.
We can obtain(10)ai−1=A1A,bi−1=A2A,ci−1=A3A,di−1=A4A,ai=A5A,bi=A6A,ci=A7A,di=A8A,ai+1=A9A,bi+1=A10A,ci+1=A11A,di+1=A12A,ai+2=A13A,bi+2=A14A,ci+2=A15A,di+2=A16A,where(11)A=A1′A2′,A1′=1−hi−1e−θi−1e−ηi−1000001τ1e−θi−1τ2e−ηi−1000000τ12e−θi−1τ22e−ηi−100001011−1hi−e−θi−e−ηi01τ1τ20−1−τ1e−θi−τ2e−ηi00τ12τ2200−τ12e−θi−τ22e−ηi00001011000001τ1τ2000000τ12τ2200000000000000000000000000000000000000000000000000001011,A2′=000000000000000000000000000000000000000000000000−1hi+1−e−θi+1−e−ηi+100000−1−τ1e−θi+1−τ2e−ηi+1000000−τ12e−θi+1−τ22e−ηi+100001011−1hi+2−e−θi+2−e−ηi+201τ1τ20−1−τ1e−θi+2−τ2e−ηi+200τ12τ2200−τ12e−θi+2−τ22e−ηi+200001011000001τ1τ2000000τ12τ2200000000.
Akk=1,2,…,16 is the matrix obtained by replacing the k column of A with 0,…,0,1T, θi=τ1hi, and ηi=τ2hi.
For even splitting, i.e.,hi=h,i=1,2,…,n, and τ2=−τ1=τ≠0, the results obtained are(12)aj−1=−ωθσ,bj−1=−ωθhσ,cj−1=12ω−coshθ+sinhθσ,dj−1=12ωcoshθ+sinhθσ,aj=θ−4sinh3θ+sinh4θ−4sinhθ+6sinh2θσ,bj=θ−3sinh3θ+2sinh2θ+sinh4θ+sinhθhσ,cj=−12σ−5+4cosh2θ−2sinh2θ−4cosh3θ+3sinh3θ+cosh4θ−sinh4θ+4coshθ−sinhθ,dj=−12σcosh4θ+sinh4θ−5−4cosh3θ−3sinh3θ+4coshθ+sinhθ+4cosh2θ+2sinh2θ,aj+1=−ωθσ,bj+1=−θ−3sinh3θ+2sinh2θ+sinh4θ+sinhθhσ,cj+1=14σ−5+4cosh2θ+10sinh2θ−4cosh3θ+cosh4θ−sinh4θ+4coshθ−16sinhθ,dj+1=14σcosh4θ+sinh4θ−5−4cosh3θ+4coshθ+16sinhθ+4cosh2θ−10sinh2θ,aj+2=0,bj+2=ωθhσ,cj+2=ω2σ,dj+2=−ω2σ,where(13)ω=5sinhθ−4sinh2θ+sinh3θ,σ=5+θsinh4θ−4θsinhθ+6θsinh2θ−4θsinh3θ+4cosh3θ−4coshθ−cosh4θ−4cosh2θ,θ=τh.
Define a function(14)Bjx=sj−1x,x∈Ij−1,sjx,x∈Ij,sj+1x,x∈Ij+1,sj+2x,x∈Ij+2,0,otherwise,j=2,3,…,n−2,where the coefficients of sj−1x,sjx,sj+1x,andsj+2x are given by the solution of function (9). Besides, for j=−1,0,1,n−1,n,n+1, an interval expansion will be conducted, i.e., x0,x1,…,xn will be extended to be(15)x−3,x−2,x−1,x0,x1,…,xn,xn+1,xn+2,xn+3.
Thus,j in (11) can take the values of −1,0,1,…,n,n+1. Of course, the function domain is still x0,xn.
For the space basis ofℰ32Δn, we have the following proposition.Proposition 1.
Whenτ1≠τ2, the function set of Bjxj=−1n+1 is a set of basis of ℰ32Δn.For special cases, we can prove the following.Theorem 1.
Whenτ2=−τ1=τ≠0 and hi=h,i=−2,−1,…,n+3, the function set of Bjxj=−1n+1 is a set of basis of ℰ32Δn.Proof.
Obviously,Bjx∈ℰ32Δn,j=−1,0,1,…,n,n+1, and the number of functions equals the dimension of ℰ32Δn. Therefore, we need to show that Bjxj=−1n+1 is linearly independent.
Forx=xi,i=−1,0,…,n,n+1, let B=βij=Bjxin+3×n+3, if matrix B is invertible, then Bjxj=−1n+1 is linearly independent. It is easy to know that B is a tridiagonal matrix, and(16)bi,i−1=ai−1+ci−1+di−1,bi,i=ai+ci+di=1,bi,i+1=ai+1+ci+1+di+1,where i=−1,0,…,n+1. Because(17)bi,i−1+bi,i+1=ai−1+ci−1+di−1+ai+1+ci+1+di+1=8σcoshθ−121−coshθ2+θsinhθ=1−coshθ2+θsinhθ1−coshθ2+θsinhθcoshθ≤12<1=bii.
B is a strictly tridiagonal matrix. Thus, B is invertible. This proves that Bjxj=−1n+1 is linearly independent.
Whenhi=h,i=−2,−1,…,n+3, Bjxj=−1n+1 has the following properties.Proposition 2.
For any ,x∈a,b(18)∑j=−1n+1Bjx=C,where C is not related to x, and(19)C=τ1τ2hτ2−τ1e−τ1h+e−τ2h−e−hτ1+τ2−1τ22−τ12e−hτ1+τ2+τ12+τ22+τ1τ22h−τ12τ2he−τ1h−τ12+τ22−τ1τ22h+τ12τ2he−τ2h+τ12−τ22.The basis functionBjx has a local supporting set of xi−2,xi+2. Figure 1 shows all functions in the function set of Bjxj=−1n+1 at the same coordinate.Figure 1
Basis function setBjxj=−1n+1 of the exponential spline function space.


## 3. Properties of Exponential Spline Functions

The relationship between the function value and the first derivative and that between the function value and the second derivative are commonly used in numerically solving DEs. Next, we will derive some properties of the functions inℰ32Δn, and these relationships will be used in numerically calculating DEs.Letsx∈ℰ32Δn; then,(20)sx=six,x∈Ii,i=1,2,…,n.Denote(21)Si=sxi,Di=s′xi,Mi=s″,i=1,2,…,n,since sixi=si+1xi,si″xi=si+1″xi, the factors of ai,bi,ci, and di can be written as(22)ai=1τ12τ22ρiτ12−τ22Mi−1+τ22e−τ2hi−τ12e−τ1hiMi+τ12τ22ρiSi,bi=1hiτ12τ22ρiτ22e−τ1hi−τ12e−τ2hi+τ12−τ22Mi−1+τ12e−τ1hi−τ2hi−τ22e−τ1hi−τ2hi−τ12e−τ1hi+τ22e−τ2hiMi+τ12τ22ρiSi+τ12τ22e−τ2hiSi−1−τ12τ22e−τ1hiSi−1,ci=1τ12ρiMi−1−e−τ2hiMi,di=−1τ22ρiMi−1−e−τ1hiMi,where ρi=e−τ1hi−e−τ2hi,i=1,2,…,n (the same applies hereinafter). Using the continuous condition of the first derivative(23)si′xi=si+1′xi,i=1,2,…,n−1,we obtain(24)α1iMi−1+α2iMi+α3iMi+1−1hiSi−1+hi+1+hihihi+1Si−1hi+1Si+1=0,where i=1,2,…,n−1,(25)α1i=−hiτ1τ22−hiτ12τ2−τ12e−τ2hi+τ22e−τ1hi+τ12−τ22τ12τ22ρihi,α2i=−1τ12τ22ρiρi+1hihi+1hi+1τ12−hi+1τ22e−τ1hi−τ2hi−τ1hi+1+hi+1τ22−hi+1τ12e−τ1hi−τ2hi−τ2hi+1+−hiτ12+hiτ22e−τ1hi+hiτ12−hiτ22e−τ2hi+hi+1τ12+hiτ12e−τ2hi+1−τ1hi+hiτ22+hi+1τ22e−τ1hi+1−τ2hi+−hi+1τ22−hiτ12+hihi+1τ1τ22−hihi+1τ12τ2e−τ2hi+1+hi+−hi+1τ12−hiτ22−hihi+1τ1τ22+hihi+1τ12τ2e−τ1hi+1+hi,α3i=τ22e−τ2hi+1−τ12e−τ1hi+1+hi+1τ12τ2−hi+1τ1τ22+τ12−τ22e−hi+1τ1+τ2τ12τ22ρi+1hi+1.Ifhi=h>0,i=1,2,…,n and τ2=−τ1=τ, then the corresponding coefficients become(26)α1i=e2τh−2τheτh−1hτ2e2τh−1,α2i=2−e2τh+τhe2τh+1+τhhτ2e2τh−1,α3i=e2τh−2τheτh−1hτ2e2τh−1,α4i=−1h,α5i=2h,α6i=−1h.Letτ⟶0; then, (24) is equivalent to(27)h2Mi−1+4Mi+Mi+1−6Si−1+12Si−6Si+1=0.In a similar way, we can obtain the relationship between the function value and the first derivative:(28)β1iDi−1+β2iDi+β3iDi+1+β4iSi−1+β5iSi+β6iSi+1=0,where(29)β1i=1KiKi+1τ1−τ22−τ1+τ1τ2hi−τ2e−hi+1τ1+τ2−τ22τ1−τ2e−τ1hi+1−τ2hi+1−τ1hi+τ12τ1−τ2e−τ1hi+1−τ2hi+1−τ2hi−τ22−τ1+τ1τ2hi+1+τ2e−τ1hi+1+hi−τ12τ1τ2hi+1+τ1−τ2e−τ2hi+1+hi+τ22τ2−τ1e−τ1hi−τ12τ2−τ1e−τ2hi+τ1−τ2−τ1+τ1τ2hi−τ2−τ1+τ1τ2hi+1+τ2e−τ1hi+1−τ1−τ2−τ1+τ1τ2hi−τ2τ1τ2hi+1+τ1−τ2e−τ2hi+1+τ12−τ1+τ1τ2hi+1+τ2e−τ1hi+1−τ2hi+τ22τ1τ2hi+1+τ1−τ2e−τ1hi−τ2hi+1+τ1−τ22−τ1+τ1τ2hi−τ2,(30)β2i=1KiKi+1−τ12τ1−τ2−1+τ2hi+1e−τ1hi+1−τ1hi−τ2hi+τ22τ1hi+1−1τ1−τ2e−τ2hi+1−τ1hi−τ2hi−τ12τ1−τ2−1+τ2hie−τ1hi+1−τ2hi+1−τ1hi+τ22τ1−τ2τ1hi−1e−τ1hi+1−τ2hi+1−τ2hi−−τ23−τ12τ2+τ13τ22Zhihi+1−τ13τ2hi−τ1τ23hi+1−τ1τ23hi+τ12τ22hi+1−τ12τ23hihi+1−τ13τ2hi+1+τ12τ22hi+τ13+τ1τ22e−τ1hi+hi+1−τ1+τ2τ1−τ22e−τ1+τ2hi+hi+1+−τ12τ2+τ1τ22−τ12τ23hihi+1+τ13+τ13τ22hihi+1+τ1τ23hi−τ23−τ12τ22hi+τ13τ2hi−τ12τ22hi+1+τ13τ2hi+1+τ1τ23hi+1e−τ2hi+hi+1−τ22τ2−τ1τ1hi+1e−τ1hi−τ12τ2hi+1τ1−τ2e−τ2hi+τ22τ1hi+1+1τ1−τ2e−τ1hi+1+τ12τ2−τ1τ2hi+1+1e−τ2hi+1−τ12τ22hi+hi+1e−τ1hi+1−τ2hi−τ12τ22hi+hi+1e−τ1hi−τ2hi+1+τ1+τ2τ1−τ22,(31)β3,i=1KiKi+1τ1−τ22τ1τ2hi+1+τ1+τ2e−hi+1τ1+τ2+τ12τ2−τ1e−τ1hi+1−τ1hi−τ2hi−τ22τ2−τ1e−τ2hi+1−τ1hi−τ2hi−τ2−τ1τ2+τ1τ2hi−τ1τ1τ2hi+1+τ1+τ2e−τ1hi+1−τ2hi+1−τ1hi−τ1−τ2τ1τ2hi+τ1−τ2τ1τ2hi+1+τ1+τ2e−τ1hi+1−τ2hi+1−τ2hi−τ12τ2+τ1τ2hi−τ1e−τ1hi+1+hi+τ1−τ22τ1τ2hi+1+τ1+τ2e−τ1+τ2hi+1+hi−τ22τ1τ2hi+τ1−τ2e−τ2hi+1+hi−τ12τ1−τ2e−τ1hi+1+τ22τ1−τ2e−τ2hi+1+τ12τ1τ2hi+τ1−τ2e−τ1hi+1−τ2hi+τ22τ2+τ1τ2hi−τ1e−τ1hi−τ2hi+1,(32)β4,i=τ1τ2KiKi+1τ1−τ22e−hi+1τ1+τ2+τ2τ1−τ2e−τ1hi+1−τ2hi+1−τ1hi−τ1τ1−τ2e−τ1hi+1−τ2hi+1−τ2hi+τ2τ1τ2hi+1+τ2−τ1e−τ1hi+1+hi+τ1τ1+τ1τ2hi+1−τ2e−τ2hi+1+hi+τ2τ1−τ2e−τ1hi−τ1τ1−τ2e−τ2hi+τ1−τ2τ1τ2hi+1+τ2−τ1e−τ1hi+1−τ1−τ2τ1+τ1τ2hi+1−τ2e−τ2hi+1−τ1τ1τ2hi+1+τ2−τ1e−τ1hi+1−τ2hi−τ2τ1+τ1τ2hi+1−τ2e−τ1hi−τ2hi+1+τ1−τ22,(33)β6,i=τ1τ2KiKi+1−τ1−τ22e−hi+1τ1+τ2+τ1τ1−τ2e−τ1hi+1−τ1hi−τ2hi−τ2τ1−τ2e−τ2hi+1−τ1hi−τ2hi−τ1−τ2−τ1+τ1τ2hi+τ2e−τ1hi+1−τ2hi+1−τ1hi+τ1−τ2−τ2+τ1τ2hi+τ1e−τ1hi+1−τ2hi+1−τ2hi+τ1−τ1+τ1τ2hi+τ2e−τ1hi+1+hi−τ1−τ22e−τ1+τ2hi+1+hi+τ2−τ2+τ1τ2hi+τ1e−τ2hi+1+hi+τ1τ1−τ2e−τ1hi+1−τ2τ1−τ2e−τ2hi+1−τ1−τ2+τ1τ2hi+τ1e−τ1hi+1−τ2hi−τ2−τ1+τ1τ2hi+τ2e−τ1hi−τ2hi+1,(34)Ki=−τ1+τ1τ2hi+τ2e−τ1hi+−τ1−τ1τ2hi+τ2e−τ2hi+τ1−τ2e−hiτ1+τ2+τ1−τ2.Ifτ2=−τ1=τ and hi=h,i=1,2,…,n, then(35)β1i=−2τheτh+e2τh−1τ4eτh−2−2e2τh+τhe2τh−τh,β2i=2τhe2τh−e2τh+τh+1τ4eτh−2−2e2τh+τhe2τh−τh,β3i=−2τeτhh+e2τh−1τ4eτh−2−2e2τh+τhe2τh−τh,β4i=−1+eτhτ2−2eτh+τeτhh+2+τh,β5i=0,β6i=−−1+eτhτ2−2eτh+τeτhh+2+τh.Letτ⟶0; then, (28) is equivalent to(36)hDi−1+4Di+Di+1+3Si−1−3Si+1=0.This is consistent with the cubic 2nd-order polynomial spline function relationship.Besides, using(37)sixi=Si,si′xi=Di,si″xi=Mi,si″xi−1=Mi−1,we can obtain(38)ai=τ12−τ22Mi−1+τ22e−τ2hi−τ12e−τ1hiMi+τ12τ22ρiSiτ12τ22ρi,bi=τ1−τ2Mi−1+τ2e−τ2hi−τ1e−τ1hiMi+τ1τ2ρiDiτ1τ2ρi,ci=−−Mi−1+e−τ2hiMiτ12ρi,di=−Mi−1+e−τ1hiMiτ22ρi.Then, using continuous condition ofsixi=si+1xi,si′xi=si+1′xi, we obtain(39)−hi+1τ1τ22−hi+1τ12τ2+τ22e−τ1hi+1−τ12e−τ2hi+1+τ12−τ22τ12τ22ρi+1Mi+τ12e−τ1hi+1−τ22e−τ2hi+1+hi+1τ1τ22e−τ2hi+1−hi+1τ12τ2e−τ1hi+1+τ22e−hi+1τ1+τ2−τ12e−hi+1τ1+τ2τ12τ22ρi+1Mi+1+Si−Si+1+hi+1Di+1=0,(40)−τ2e−τ2hi+1−τ1e−τ1hi+1τ1e−hi+1τ1+τ2−τ2e−hi+1τ1+τ2τ1τ2ρi+1Mi+1−τ2e−τ1hi+1−τ1e−τ2hi+1+τ1−τ2τ1τ2ρi+1Mi+Di−Di+1=0.Solving (39) yields(41)Di=−τ12e−τ2hi−τ22e−τ1hi−hiτ1τ22+hiτ12τ2−τ12+τ22hiτ12τ22ρiMi−1+−hiτ1τ22+τ22e−τ2hi+−τ12+hiτ12τ2e−τ1hi+τ12−τ22e−hiτ2+τ1hiτ12τ22ρiMi−Si−1hi+Sihi.Meanwhile, by eliminatingDi+1 with (39) and (40) and rearranging, we obtain(42)Di=−τ12−hi+1τ12τ2e−τ2hi+1+τ22+hi+1τ1τ22e−τ1hi+1−τ22+τ12hi+1τ12τ22ρi+1Mi−−τ12+hi+1τ1τ22−hi+1τ12τ2+τ22e−hi+1τ1+τ2+τ12e−τ1hi+1−τ22e−τ2hi+1hi+1τ12τ22ρi+1Mi+1,(43)−Sihi+1+Si+1hi+1.We can also usesi+1xi+1=Si+1,si+1xi=Si,si+1′xi+1=Di+1,andsi+1′xi=Di to obtain(44)ai=1ϱiτ1hi−τ2hi+ρiDi−1+hiτ2e−τ2hi−τ1e−τ1hihi−ρiDi+τ2e−τ2hi+τ1−τ2−τ1e−τ1hiSi−1+−τ1−τ1hiτ2e−τ2hi+τ1hiτ2+τ2e−τ1hi+τ1−τ2e−hiτ1+τ2Si,bi=1ϱiτ1−τ2−τ1e−τ2hi+τ2e−τ1hiDi−1+τ2e−τ2hi−τ1e−τ1hi+τ1−τ2e−hiτ1+τ2Di−τ1τ2ρiSi−1+τ1τ2ρiSi,ci=1ϱie−τ2hi−1+τ2hiDi−1+−e−τ2hi+1−hiτ2e−τ2hiDi+−τ2e−τ2hi+τ2Si−1+−τ2+τ2e−τ2hiSi,di=1ϱi−τ1hi−e−τ1hi+1Di−1+τ1e−τ1hihi+e−τ1hi−1Di+−τ1+τ1e−τ1hiSi−1+−τ1e−τ1hi+τ1Si,where ϱi=−τ1hiτ2+τ2−τ1e−τ2hi+τ1hiτ2+τ2−τ1e−τ1hi+τ1−τ2e−hiτ1+τ2+τ1−τ2. Fast Hermite interpolation can achieved by using this set of relations.


## 4. Exponential Spline Difference Method

The spline difference method uses the relationship between the spline function and its derivative to construct the differential expression to numerically solve DEs, by which the numerical solution at nodes can be obtained, and that within the subintervals can also be calculated by using the spline function expressions. It is the advantage of this method compared with the general difference schemes. In fact, it can be said that the approximate analytical solution using the splines is obtained.


### 4.1. Differential Expression

The following presents a spline difference method for solving (1) which satisfies one of the boundary conditions (2)–(5) for the boundary value problems of ODEs. Due to the limitation of this method, we only consider the case of qx=0 in this section. For convenience of description, we first consider the boundary condition (2). From (1), we can obtain(45)u″=gx−pxu,x∈a,b.By discretization the above equation, we obtain(46)ui″=gi−piui,i=0,1,2,…,n,where gi=gxiandpi=pxi.SubstituteS with u and M with u″ in (24), and we obtain(47)ti=α1iui−1″+α2iui″+α3iui−1″−ui−1hi+hi+hi+1uihihi+1−ui+1hi+1,where i=1,…,n−1,ui=uxiandui″=u″xi. ti is the local truncation error at xi. By substituting (46) into (47) and rearranging, we obtain(48)−1hi−α1ipi−1ui−1+hi+hi+1hihi+1−α2ipiui+−1hi+1−α3ipi+1ui+1+α1igi−1+α2igi+α3igi+1=ti,where i=1,2,…,n−1. Thus, we get n−1 equations about u0,u1,…,un.The following equation can also be derived from the boundary conditions (2):(49)u0=μ0,un=μ1.So, there will ben+1 equations in (48) and (49) in total, and it matrix form can be written as(50)AU=F+T,where A=B−W, with(51)B=100−1h1h1+h2h1h2−1h2−1h2h2+h3h2h3−1h3⋱−1hn−2hn−2+hn−1hn−2hn−1−1hn−1−1hn−1hn−1+hnhn−1hn−1hn001,(52)W=000α11p0α21p1α31p2α12p1α22p2α32p3⋱α1,n−2pn−3α2,n−2pn−2α3,n−2pn−1α1,n−1pn−2α2,n−1pn−1α3,n−1pn000,(53)F=μ0,−α1igi−1+α2igi+α3igi+1i=1n−1,μ1T,T=tii=0nT,U=uii=0nT.IfA is nonsingular, then equation (50) has a unique solution. Solving (50) yields the approximate values of ui in ux at the splitting points xi,i=0,1,…,n. Since A is a tridiagonal matrix, the catch-up method can be used to reduce the calculation in practice.To obtain the spline function expression, we can calculateui″ from (2) after finding ui and substitute them into (22) to obtain ai,bi,ci,anddi so that the approximate solution of the spline over the entire interval will be found.The case of other boundary conditions will be discussed below. For the second boundary conditions, from (41), (42), and (3), we obtain(54)α20M0+α30M1+S0h1−S1h1=−μ0,(55)α1nMn−1+α2nMn−Sn−1hn+Snhn=μ1,where(56)α20=τ12+h1τ12τ2e−τ2h1+−τ22−τ1h1τ22e−τ1h1−τ12+τ22h1τ12τ22ρ1,α30=−τ22e−τ2h1−τ12e−τ1h1+τ12−τ1h1τ22+h1τ12τ2−τ22e−h1τ2+τ1h1τ12τ22ρ1,α1n=−τ12e−τ2hn−τ22e−τ1hn−hnτ1τ22+hnτ12τ2−τ12+τ22hnτ12τ22ρn,α2n=−hnτ1τ22+τ22e−τ2hn+−τ12+hnτ12τ2e−τ1hn+τ12−τ22e−hnτ2+τ1hnτ12τ22ρn.Solving this type of boundary problem only needs to modify the first line and the last row ofB, W, and F in (47), namely,(57)B0=1h1−1h10⋯0,W0=α20p0α30p10⋯0,Bn=0⋯0−1hn1hn,Wn=0⋯0αn1pn−1αn2pn,F0=−μ0−α20g0+α30g1,Fn=μ1−α1ngn−1+α2ngn.And the third and the fourth boundary conditions should also be modified accordingly.


### 4.2. Error Estimation

Supposeux is sufficiently smooth in a,b, by Taylor expanding ui−1,ui+1,ui−1″,ui+1″ in difference equation (47) at xi, we obtain(58)ti=∑j=26ξjuij+Ohi7,where(59)ξ2=α1i+α2i+α3i−12hi+hi+1,ξ3=−α1ihi+α3ihi+1+16hi2−hi+12,ξ4=12α1ihi2+12α3ihi+12−124hi3+hi+13,ξ5=−16α1ihi3+16α3ihi+13+1120hi4+hi+14,ξ6=124α1ihi4+124α3ihi+14−1720hi5+hi+15.Letτ1⟶0 and τ2⟶0, we find ξ2=ξ3=0 and ξ4=1/24hi+13+hi3. When(60)α1i=112hi2+hihi+1−hi+12hi,α2i=112hi3+4hi2hi+1+4hihi+12+hi+13hihi+1,α3i=−112hi2−hihi+1−hi+12hi+1,we obtain ξ2=ξ3=ξ4=0,(61)ξ5=1360−hi+hi+12hi+hi+1hi+2hi+1hi+hi+1.ξ5=0 for an even splitting, i.e., hi=h,i=1,2,…,n, and ξ6=1/240h6. This indicates that the even splitting has higher accuracy if ux is sufficiently smooth.


### 4.3. Convergence Analysis

We mainly discuss the convergence of equation (47) and the differential expression of (49) in the sense of E∞.Lemma 1 (see [14]).
If then-order matrix B satisfies one of the following two conditions:(1)
B is a strictly diagonally dominant matrix(2)
B is a second half strong diagonally dominant matrix
Then,B is nonsingular and ρI−D−1B<1, where D=diagB and I is the identity matrix.By (14),(62)AS=F.IfA is reversible, combining (50), we obtain(63)E=A−1T=B−W−1T=I−B−1W−1B−1T,where E=ei=U−S, and B, W, and T are given by (51)–(53). The reversibility of B can be proved by 1. From the boundary condition (2), we can get e0=en=0. Thus, discussing the convergence of E∞ is consistent with that of eii=1n−1∞. So, in the following convergence discussion, the first and the last row and the first and the last column will be removed from the original matrices of A, B, and W to obtain the n−1-order ones, keeping the subscript value unchanged.If(64)B−1∞W∞<1,then(65)E∞≤B−1∞T∞1−B−1∞W∞≤B−1∞T∞.To calculateB−1, the following lemma is needed:Lemma 2 (see [15]).
Let the square matrixA be an n-order tridiagonal one, with the following expression:(66)A=a1b1c1a2b2c2⋱⋱⋱⋱bn−1cn−1an,and bici≠0,i=1,2,…,n−1, ai≠0,i=1,2,…,n; then, the expression of A−1=aij−1 can be written as(67)aij−1=−1i+j∏k=ij−1bkdetA1,⋯,i−1detAj+1,⋯,ndetA,i<j,−1i+j∏k=ji−11ckdetA1,⋯,i−1detAj+1,⋯,ndetA,i≥j,where Ai1,i2,…,ik denotes the matrix consisting of the elements of i1,i2,…,ik rows crossed with the columns in A. Specifically, A1,…,i−1=1 for i=1.
Using (2), B−1=bij−1 in (65) can be written as (with ∑k=1nhk=b−a)(68)bij−1=1b−a∑k=1ihk∑k=j+1nhk,i<j,1b−a∑k=1jhk∑k=i+1nhk,i≥j,then 0<bij−1≤b−a/4. Thus,(69)B−1∞≤14nb−a≤b−a24hmin.According to (27), we find ti=Ohi2. Thus,(70)B−1∞T∞≤b−a24hmax2hmin.So, ifhmax=Ohmin, then B−1∞W∞⟶0 when hmax⟶0.So far, we have proved thatTheorem 2.
For givenτ1,τ2, and Δ:0=x0<x1<⋯<xn=1, hi=xi−xi−1,i=1,2,…,n, hmax=max1≤i≤nhi,hmin=min1≤i≤nhi, if B−1∞W∞<1, then(71)E∞=Ohmax2hmin,that is, the differential expression of (47) is convergent if hmax=Ohmin.


## 4.1. Differential Expression

The following presents a spline difference method for solving (1) which satisfies one of the boundary conditions (2)–(5) for the boundary value problems of ODEs. Due to the limitation of this method, we only consider the case of qx=0 in this section. For convenience of description, we first consider the boundary condition (2). From (1), we can obtain(45)u″=gx−pxu,x∈a,b.By discretization the above equation, we obtain(46)ui″=gi−piui,i=0,1,2,…,n,where gi=gxiandpi=pxi.SubstituteS with u and M with u″ in (24), and we obtain(47)ti=α1iui−1″+α2iui″+α3iui−1″−ui−1hi+hi+hi+1uihihi+1−ui+1hi+1,where i=1,…,n−1,ui=uxiandui″=u″xi. ti is the local truncation error at xi. By substituting (46) into (47) and rearranging, we obtain(48)−1hi−α1ipi−1ui−1+hi+hi+1hihi+1−α2ipiui+−1hi+1−α3ipi+1ui+1+α1igi−1+α2igi+α3igi+1=ti,where i=1,2,…,n−1. Thus, we get n−1 equations about u0,u1,…,un.The following equation can also be derived from the boundary conditions (2):(49)u0=μ0,un=μ1.So, there will ben+1 equations in (48) and (49) in total, and it matrix form can be written as(50)AU=F+T,where A=B−W, with(51)B=100−1h1h1+h2h1h2−1h2−1h2h2+h3h2h3−1h3⋱−1hn−2hn−2+hn−1hn−2hn−1−1hn−1−1hn−1hn−1+hnhn−1hn−1hn001,(52)W=000α11p0α21p1α31p2α12p1α22p2α32p3⋱α1,n−2pn−3α2,n−2pn−2α3,n−2pn−1α1,n−1pn−2α2,n−1pn−1α3,n−1pn000,(53)F=μ0,−α1igi−1+α2igi+α3igi+1i=1n−1,μ1T,T=tii=0nT,U=uii=0nT.IfA is nonsingular, then equation (50) has a unique solution. Solving (50) yields the approximate values of ui in ux at the splitting points xi,i=0,1,…,n. Since A is a tridiagonal matrix, the catch-up method can be used to reduce the calculation in practice.To obtain the spline function expression, we can calculateui″ from (2) after finding ui and substitute them into (22) to obtain ai,bi,ci,anddi so that the approximate solution of the spline over the entire interval will be found.The case of other boundary conditions will be discussed below. For the second boundary conditions, from (41), (42), and (3), we obtain(54)α20M0+α30M1+S0h1−S1h1=−μ0,(55)α1nMn−1+α2nMn−Sn−1hn+Snhn=μ1,where(56)α20=τ12+h1τ12τ2e−τ2h1+−τ22−τ1h1τ22e−τ1h1−τ12+τ22h1τ12τ22ρ1,α30=−τ22e−τ2h1−τ12e−τ1h1+τ12−τ1h1τ22+h1τ12τ2−τ22e−h1τ2+τ1h1τ12τ22ρ1,α1n=−τ12e−τ2hn−τ22e−τ1hn−hnτ1τ22+hnτ12τ2−τ12+τ22hnτ12τ22ρn,α2n=−hnτ1τ22+τ22e−τ2hn+−τ12+hnτ12τ2e−τ1hn+τ12−τ22e−hnτ2+τ1hnτ12τ22ρn.Solving this type of boundary problem only needs to modify the first line and the last row ofB, W, and F in (47), namely,(57)B0=1h1−1h10⋯0,W0=α20p0α30p10⋯0,Bn=0⋯0−1hn1hn,Wn=0⋯0αn1pn−1αn2pn,F0=−μ0−α20g0+α30g1,Fn=μ1−α1ngn−1+α2ngn.And the third and the fourth boundary conditions should also be modified accordingly.


## 4.2. Error Estimation

Supposeux is sufficiently smooth in a,b, by Taylor expanding ui−1,ui+1,ui−1″,ui+1″ in difference equation (47) at xi, we obtain(58)ti=∑j=26ξjuij+Ohi7,where(59)ξ2=α1i+α2i+α3i−12hi+hi+1,ξ3=−α1ihi+α3ihi+1+16hi2−hi+12,ξ4=12α1ihi2+12α3ihi+12−124hi3+hi+13,ξ5=−16α1ihi3+16α3ihi+13+1120hi4+hi+14,ξ6=124α1ihi4+124α3ihi+14−1720hi5+hi+15.Letτ1⟶0 and τ2⟶0, we find ξ2=ξ3=0 and ξ4=1/24hi+13+hi3. When(60)α1i=112hi2+hihi+1−hi+12hi,α2i=112hi3+4hi2hi+1+4hihi+12+hi+13hihi+1,α3i=−112hi2−hihi+1−hi+12hi+1,we obtain ξ2=ξ3=ξ4=0,(61)ξ5=1360−hi+hi+12hi+hi+1hi+2hi+1hi+hi+1.ξ5=0 for an even splitting, i.e., hi=h,i=1,2,…,n, and ξ6=1/240h6. This indicates that the even splitting has higher accuracy if ux is sufficiently smooth.


## 4.3. Convergence Analysis

We mainly discuss the convergence of equation (47) and the differential expression of (49) in the sense of E∞.Lemma 1 (see [14]).
If then-order matrix B satisfies one of the following two conditions:(1)
B is a strictly diagonally dominant matrix(2)
B is a second half strong diagonally dominant matrix
Then,B is nonsingular and ρI−D−1B<1, where D=diagB and I is the identity matrix.By (14),(62)AS=F.IfA is reversible, combining (50), we obtain(63)E=A−1T=B−W−1T=I−B−1W−1B−1T,where E=ei=U−S, and B, W, and T are given by (51)–(53). The reversibility of B can be proved by 1. From the boundary condition (2), we can get e0=en=0. Thus, discussing the convergence of E∞ is consistent with that of eii=1n−1∞. So, in the following convergence discussion, the first and the last row and the first and the last column will be removed from the original matrices of A, B, and W to obtain the n−1-order ones, keeping the subscript value unchanged.If(64)B−1∞W∞<1,then(65)E∞≤B−1∞T∞1−B−1∞W∞≤B−1∞T∞.To calculateB−1, the following lemma is needed:Lemma 2 (see [15]).
Let the square matrixA be an n-order tridiagonal one, with the following expression:(66)A=a1b1c1a2b2c2⋱⋱⋱⋱bn−1cn−1an,and bici≠0,i=1,2,…,n−1, ai≠0,i=1,2,…,n; then, the expression of A−1=aij−1 can be written as(67)aij−1=−1i+j∏k=ij−1bkdetA1,⋯,i−1detAj+1,⋯,ndetA,i<j,−1i+j∏k=ji−11ckdetA1,⋯,i−1detAj+1,⋯,ndetA,i≥j,where Ai1,i2,…,ik denotes the matrix consisting of the elements of i1,i2,…,ik rows crossed with the columns in A. Specifically, A1,…,i−1=1 for i=1.
Using (2), B−1=bij−1 in (65) can be written as (with ∑k=1nhk=b−a)(68)bij−1=1b−a∑k=1ihk∑k=j+1nhk,i<j,1b−a∑k=1jhk∑k=i+1nhk,i≥j,then 0<bij−1≤b−a/4. Thus,(69)B−1∞≤14nb−a≤b−a24hmin.According to (27), we find ti=Ohi2. Thus,(70)B−1∞T∞≤b−a24hmax2hmin.So, ifhmax=Ohmin, then B−1∞W∞⟶0 when hmax⟶0.So far, we have proved thatTheorem 2.
For givenτ1,τ2, and Δ:0=x0<x1<⋯<xn=1, hi=xi−xi−1,i=1,2,…,n, hmax=max1≤i≤nhi,hmin=min1≤i≤nhi, if B−1∞W∞<1, then(71)E∞=Ohmax2hmin,that is, the differential expression of (47) is convergent if hmax=Ohmin.


## 5. Example

Example 1.
Solve the following singular boundary value problem (see [7]):(72)−εu″+pxu=gx,0≤x≤1,u0=u1=0.where px=1+x1−x and gx=1+x1−x+2ε−x1−x2e−x/ε+2ε−x21−xe−1−x/ε, and its analytical solution is(73)ux=1+x−1e−x/ε−xe−1−x/ε.
Table1 lists the results of the exponential spline difference method for ε=0.1 and that of the method proposed in [7].
Figure2 shows the results for ε=0.1, and n=5 and 10, respectively. To show the difference between the exact solution and the numerical one, we deliberately use fewer split points and simply connect two adjacent solutions with straight lines. A smoother exponential spline function which is closer to the exact solution can be constructed using these obtained numerical solutions.Table 1
The maximum errorE∞ in solving (72) with the exponential spline difference method.
n=32n=64n=128n=256n=512Ref. [7]1.84e − 44.61e − 51.15e − 52.88e − 67.21e − 7Our method (τ=n1/4,−n1/4)8.11e − 62.03e − 65.09e − 71.28e − 74.25e − 8Figure 2
Results of the boundary value problem (72) for ε=0.1, and n=5,10, and 20, respectively.Example 2.
Consider the convection-dominated equation (see [11]):(74)u″−εu′−εu=−ε,0<x<1,u0=u1=0.
The analytical solution is given by(75)ux=1+e1−1+4/λ/2/λ−1e1+1+4/λ/2/λx+1−e1+1+4/λ/2/λe1−1+4/λ/2/λx1+1+4/λ/2/λ−1−1+4/λ/2/λ.
The computational results are shown in Table2 for ε=0.01 and various values of n (Figure 3).Table 2
The maximum errorE∞ and root mean square error for problem (74).
n=10n=20n=40n=80Ref. [11]1.02e − 32.96e − 53.69e − 65.11e − 8Our method (τ=10.11,−0.012)8.57e − 78.18e − 78.03e − 77.97e − 7Figure 3
Results of the boundary value problem (74) for ε=0.01, and n=5,10, and 20, respectively.Example 3.
Consider the convection-dominated equation (see [16–19]) Figure 4:(76)εu″+u=0,0<x<1,u0=0,u1=1,(77)ux=sinx/εsin1/ε,ε≠nπ−2.
The computational results are shown in Table3 for ε=1/2 and various values of n.Figure 4
Results of the boundary value problem (76) for ε=1/32, and n=5,10, and 20, respectively.Table 3
The maximum errorE∞ and root mean square error for problem (76).
εMethodn=16n=32n=64n=128n=256n=512n=10241/2Ref. [16]2.6e − 46.5e − 51.6e − 54.1e − 61.0e − 62.6e − 76.4e − 8Our method(τ=1.18,−12.75)2.34e − 46.37e − 51.65e − 54.20e − 61.06e − 62.72e − 75.21e − 81/4Ref. [16]1.9e − 34.7e − 41.2e − 42.9e − 57.4e − 61.8e − 64.6e − 7Our method(τ=1.91,−11.77)3.04e − 38.19e − 42.11e − 45.38e − 51.35e − 53.36e − 61.94e − 71/8Ref. [16]7.3e − 21.9e − 24.7e − 31.2e − 32.9e − 47.3e − 51.8e − 5Our method()8.23e − 22.39e − 26.13e − 31.54e − 33.93e − 49.66e − 56.74e − 71/16Ref. [16]2.9e − 27.1e − 31.8e − 34.4e − 41.1e − 42.8e − 56.9e − 6Our method(τ=13,−1)1.61e − 24.03e − 31.04e − 32.67e − 46.80e − 51.66e − 53.74e − 61/32Ref. [16]—3.8e—29.8e—32.5e—36.2e—41.5e—43.9e-5Our method(τ=1,−13)7.52e − 21.99e − 25.06e − 31.30e − 33.25e − 47.93e − 51.12e − 51/64Ref. [16]—3.4e − 28.7e − 32.2e − 35.5e − 41.4e − 43.4e − 5Our method(τ=12.51,−1.48)8.54e − 22.33e − 26.39e − 31.63e − 34.10e − 41.03e − 41.85e − 51/128Ref. [16]——2.8e − 27.2e − 31.8e − 34.5e − 41.1e − 4Our method(τ=11.53,−2.46)—6.46e − 21.69e − 24.33e − 31.09e − 32.73e − 45.70e − 51/256Ref. [16]————6.3e − 21.6e − 23.9e − 3Our method(τ=13,−1)————3.21e − 28.09e − 39.13e − 41/512Ref. [16]————4.2e − 21.0e − 22.5e − 3Our method(τ=11.53,−2.46)————2.12e − 25.35e − 36.92e − 41/1024Ref. [16]—————3.4e − 28.5e − 3Our method(τ=3.93,−10)—————1.92e − 22.55e − 3


## 6. Conclusions

Exponential spline, which is a generalization of polynomial spline, is an ideal function approximation tool due to its excellent curve fitting ability. High accuracy can be achieved in solving second-order ODEs using the exponential spline scheme. The spline difference method is an ideal scheme because it can give not only the numerical results but also the spline function expressions by reusing these numerical results at the same time, whereas it is only suitable for solving certain types of equations and does not have generality. And the selection of the appropriate parameters is also needed for this method, but there is no better guideline for the selecting.



---

*Source: 7056254-2020-08-28.xml*