Geometric Mesh Three-Point Discretization for Fourth-Order Nonlinear Singular Differential Equations in Polar System

Abstract

Numerical method based on three geometric stencils has been proposed for the numerical solution of nonlinear singular fourth-order ordinary differential equations. The method can be easily extended to the sixth-order differential equations. Convergence analysis proves the third-order convergence of the proposed scheme. The resulting difference equations lead to block tridiagonal matrices and can be easily solved using block Gauss-Seidel algorithm. The computational results are provided to justify the usefulness and reliability of the proposed method.

1. Introduction

Consider the fourth-order boundary value problem: 𝑢⁽⁴⁾(𝑟)=𝑓𝑟,𝑢(𝑟),𝑢⁽¹⁾(𝑟),𝑢⁽²⁾(𝑟),𝑢⁽³⁾(𝑟),𝛼<𝑟<𝛽 (1) subject to the necessary boundary conditions: 𝑢(𝛼)=𝛼₀,𝑢⁽²⁾(𝛼)=𝛼₁,𝑢𝛽=𝛽₀,𝑢⁽²⁾𝛽=𝛽₀, (2)

where 𝑓∈𝐶⁴[𝛼,𝛽] and 𝛼₀, 𝛼₁, 𝛽₀, and 𝛽₁ are real constants and −∞<𝛼<𝑟<𝛽<∞.

Or equivalently 𝑢⁽²⁾(𝑟)(𝑟),𝑣=𝑣₍₂₎(𝑟)=𝑓𝑟,𝑢(𝑟),𝑢⁽¹⁾(𝑟),𝑣(𝑟),𝑣⁽¹⁾(𝑟),𝛼<𝑟<𝛽 (3)

subject to the natural boundary conditions: 𝑢(𝛼)=𝛼₀,𝑣(𝛼)=𝛼₁,𝑢𝛽=𝛽₀,𝑣𝛽=𝛽₁. (4)

Fourth-order differential equations occur in various areas of mathematics such as viscoelastic and inelastic flows, beam theory, Lifshitz point in phase transition physics (e.g., nematic liquid crystal, crystals, and ferroelectric crystals) [1], the rolls in a Rayleigh-Benard convection cell (two parallel plates of different temperature with a liquid in between) [2], spontaneous pattern formation in second-order materials (e.g., polymeric fibres) [3], the waves on a suspension bridge [4, 5], geological folding of rock layers [6], buckling of a strut on a nonlinear elastic foundation [7], traveling water waves in a shallow channel [8], pulse propagation in optical fibers [9], system of two reaction diffusion equation [10], and so forth.

The existence and uniqueness of the solution for the fourth and higher-order boundary value problems have been discussed in [11–14]. In the recent past, the numerical solution of fourth-order differential equations has been developed using multiderivative, finite element method, Ritz method, spline collocation, and finite difference method [15–18]. The determination of eigen values of self adjoint fourth-order differential equations was developed in [19] using finite difference scheme. The motivation of variable mesh technique for differential equations arises from the theory of electrochemical reaction-convection-diffusion problems in one-dimensional space geometry [20]. The geometric mesh method for self-adjoint singular perturbation problems using finite difference approximations was discussed in [21]. The use of geometric mesh in the context of boundary value problems was studied extensively in [22–24]. In this paper, we derive a geometric mesh finite difference method for the solution of fourth- and sixth-order differential boundary value problems with order of accuracy being three. The simplicity of the proposed method lies in its three-point discretization without any use of fictitious nodes. The scheme is compact and applicable to both singular and nonsingular problems. The resulting difference equations are solved using block Gauss-Seidel algorithm for linear case, and corresponding Newton’s method has been applied to nonlinear problems.

The paper is outlined in the following manner: in Section 2, the derivation of the method is discussed in detail. In Section 3, we define the procedure for numerical solution to singular problems in such a way that the method retains the order and accuracy even in the vicinity of singularity. In Section 4, algorithmic details are provided for the numerical solution of sixth-order differential equations. The convergence property has been discussed briefly in Section 4. The numerical illustrations based on geometric mesh as well as uniform mesh were provided in Section 5. The paper is concluded in the last section with future development and remarks.

2. Derivation of the Numerical Scheme

We discretize the solution region [𝛼,𝛽] such that 𝛼=𝑟₀<𝑟₁<⋯<𝑟_𝑛<𝑟_𝑛+1=𝛽. Let ℎ_𝑘=𝑟_𝑘−𝑟_𝑘−1 and 𝑘=1(1)𝑛+1 be the nonuniform step size, and let 𝜎_𝑘=ℎ_𝑘+1/ℎ_𝑘>0 be the geometric mesh ratio. Let 𝑈_𝑘=𝑢(𝑟_𝑘) and 𝑉_𝑘=𝑣(𝑟_𝑘) denote the exact solution values of 𝑢 and 𝑣 at the mesh 𝑟_𝑘 and 𝑢_𝑘, 𝑣_𝑘 be their approximate solutions, respectively.

Consider the following three-point geometric mesh discretizations for 𝑢⁽²⁾(𝑟)=𝑣(𝑟) and 𝑣⁽²⁾(𝑟)=𝑓(𝑟): 𝑈_𝑘+1−1+𝜎_𝑘𝑈_𝑘+𝜎_𝑘𝑈_𝑘−1+ℎ²𝑝12^𝑘𝑉_𝑘+1+𝑞_𝑘𝑉_𝑘+𝑟_𝑘𝑉_𝑘−1=𝑂ℎ⁵_𝑘(1)𝑉,𝑘=1𝑛,_𝑘+1−1+𝜎_𝑘𝑉_𝑘+𝜎_𝑘𝑉_𝑘−1+ℎ²𝑝12^𝑘𝑓_𝑘+1+𝑞_𝑘𝑓_𝑘+𝑟_𝑘𝑓_𝑘−1=𝐸_𝑘,𝑘=1(1)𝑛, (5)(6) where 𝐸_𝑘=ℎ⁵_𝑘𝜎360^𝑘2𝜎_𝑘+1𝜎_𝑘+21−𝜎²_𝑘𝑓_𝑘⁽⁵⁾(𝜉)+𝑂ℎ⁶_𝑘,𝑟_𝑘−1<𝜉<𝑟_𝑘+1,𝑝_𝑘=1−𝜎_𝑘−𝜎²_𝑘,𝑞_𝑘=−𝜎_𝑘+1𝜎²_𝑘+3𝜎_𝑘,𝑟+1_𝑘=𝜎_𝑘𝜎²_𝑘−𝜎_𝑘−1. (7)

Define 𝑈^′_𝑘=𝑈_𝑘+1−1−𝜎²_𝑘𝑈_𝑘−𝜎²_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘=𝑉_𝑘+1−1−𝜎²_𝑘𝑉_𝑘−𝜎²_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑈^′_𝑘+1=1+2𝜎_𝑘𝑈_𝑘+1−1+𝜎_𝑘²𝑈_𝑘+𝜎²_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘+1=1+2𝜎_𝑘𝑉_𝑘+1−1+𝜎_𝑘²𝑉_𝑘+𝜎²_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑈^′_𝑘−1=−𝑈_𝑘+1+1+𝜎_𝑘²𝑈_𝑘−𝜎_𝑘2+𝜎_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘−1=−𝑉_𝑘+1+1+𝜎_𝑘²𝑉_𝑘−𝜎_𝑘2+𝜎_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑓_𝑘±1=𝑓𝑟_𝑘±1,𝑈_𝑘±1,𝑈^′_𝑘±1,𝑉_𝑘±1,𝑉^′_𝑘±1. (8)(9)

It is easy to verify that 𝑓_𝑘+1=𝑓_𝑘+1−ℎ²_𝑘6𝜎_𝑘𝜎_𝑘+1𝐺_𝑘𝑈_𝑘⁽³⁾+𝐻_𝑘𝑉_𝑘⁽³⁾+𝑂ℎ³_𝑘,𝑓_𝑘−1=𝑓_𝑘−1−ℎ²_𝑘6𝜎_𝑘+1𝐺_𝑘𝑈_𝑘⁽³⁾+𝐻_𝑘𝑉_𝑘⁽³⁾+𝑂ℎ³_𝑘, (10) where 𝐺_𝑘=(𝜕𝑓/𝜕𝑈⁽¹⁾)_𝑟𝑘, 𝐻_𝑘=(𝜕𝑓/𝜕𝑉⁽¹⁾)_𝑟𝑘, and so forth.

Now, let 𝑈^′_𝑘=𝑈^′_𝑘+𝛾ℎ_𝑘𝑉_𝑘+1−𝑉_𝑘−1,𝑉^′_𝑘=𝑉^′_𝑘+𝛿ℎ_𝑘𝑓_𝑘+1−𝑓_𝑘−1, (11) where 𝛾 and 𝛿 are free parameters to be determined.

With the help of (8) and (11), it follows that 𝑈^′_𝑘=𝑈_𝑘⁽¹⁾+ℎ²_𝑘6𝜎_𝑘+6𝛾1+𝜎_𝑘𝑈_𝑘⁽³⁾+𝑂ℎ³_𝑘,𝑉^′_𝑘=𝑉_𝑘⁽¹⁾+ℎ²_𝑘6𝜎_𝑘+6𝛿1+𝜎_𝑘𝑉_𝑘⁽³⁾+𝑂ℎ³_𝑘. (12)

Further, we define 𝑓_𝑘=𝑓𝑟_𝑘,𝑈_𝑘,𝑈^′_𝑘,𝑉_𝑘,𝑉^′_𝑘. (13)

With the help of (12)-(13), it follows that 𝑓_𝑘=𝑓_𝑘+ℎ²_𝑘6𝜎_𝑘+6𝛾1+𝜎_𝑘𝐺_𝑘𝑈_𝑘⁽³⁾+ℎ²_𝑘6𝜎_𝑘+6𝛿1+𝜎_𝑘𝐻_𝑘𝑉_𝑘⁽³⁾+𝑂ℎ³_𝑘. (14)

With the help of approximations (10), (14), the difference equation (6) at each internal mesh 𝑟_𝑘, is approximated as follows: 𝑉_𝑘+1−1+𝜎_𝑘𝑉_𝑘+𝜎_𝑘𝑉_𝑘−1+ℎ²_𝑘𝑝12^𝑘𝑓_𝑘+𝑞_𝑘𝑓_𝑘+𝑟_𝑘𝑓_𝑘−1=𝐸_𝑘,𝑘=1(1)𝑛, (15) where 𝐸_𝑘ℎ=−^4𝑘𝑈72^𝑘(3)𝐺_𝑘1+𝜎_𝑘6𝛾+1+24𝛾𝜎_𝑘1+𝜎_𝑘+1+6𝛾𝜎³_𝑘−ℎ⁴_𝑘𝑉72^𝑘(3)𝐻_𝑘×1+𝜎_𝑘6𝛿+(1+24𝛿)𝜎_𝑘+(1+6𝛿)𝜎³_𝑘+𝑂ℎ⁵_𝑘,𝜎_𝑘≠1. (16)

The difference scheme (15) to be of 𝑂(ℎ³_𝑘), the coefficients of ℎ⁴_𝑘 in (16) must be zero, and hence we obtain 6𝛾+1+24𝛾𝜎_𝑘1+𝜎_𝑘+1+6𝛾𝜎³_𝑘()𝜎=0,6𝛿+1+24𝛿_𝑘1+𝜎_𝑘+(1+6𝛿)𝜎³_𝑘=0. (17)

Thus the values associated with (17) are given by 𝛾=𝛿=−𝜎_𝑘(𝜎²_𝑘+𝜎_𝑘+1)/(6(1+𝜎_𝑘)(𝜎²_𝑘+3𝜎_𝑘+1)), and the local truncation errors given by (16) becomes 𝐸_𝑘=𝑂(ℎ⁵_𝑘), 𝜎_𝑘≠1. However, for 𝜎_𝑘≠1, the error reduces to 𝐸_𝑘=𝑂(ℎ⁶_𝑘).

3. Application to Singular Problem

Consider the singular fourth-order linear differential equation in cylindrical polar coordinates: Δ⁴𝑢(𝑟)≡𝑑²𝑑𝑟²+𝜆𝑟𝑑𝑑𝑟²𝑢(𝑟)=𝑔(𝑟),0<𝑟<1 (18) or equivalently 𝑢⁽²⁾(𝑟)(𝑟),𝑣=𝑣₍₂₎(𝑟)=𝑎(𝑟)𝑣⁽¹⁾(𝑟)+𝑏(𝑟)𝑣(𝑟)+𝑐(𝑟)𝑢⁽¹⁾(𝑟)+𝑔(𝑟), (19)(20)

where 𝑎(𝑟)=−2𝜆/𝑟, 𝑏(𝑟)=−𝜆(𝜆−2)/𝑟², and 𝑐(𝑟)=𝜆(𝜆−2)/𝑟³ are the singular coefficients. The associated boundary conditions are given by (4). For 𝜆=0, 1, or 2, the differential equation (18) shows planar, cylindrical, or spherical geometries (see, [25, 26]).

Applying the difference schemes (5) and (15) to (19) and (20), respectively, we obtain a system of coupled difference equations for 𝑘=1(1)𝑛: 𝑠𝑏_𝑘¹¹𝑠𝑏_𝑘¹²𝑠𝑏_𝑘²¹𝑠𝑏_𝑘²²𝑈_𝑘−1𝑉_𝑘−1+𝑑𝑔_𝑘¹¹𝑑𝑔_𝑘¹²𝑑𝑔_𝑘²¹𝑑𝑔_𝑘²²𝑈_𝑘𝑉_𝑘+𝑠𝑝_𝑘¹¹𝑠𝑝_𝑘¹²𝑠𝑝_𝑘²¹𝑠𝑝_𝑘²²𝑈_𝑘+1𝑉_𝑘+1=𝑅_𝑘¹¹𝑅_𝑘²¹, (21) where 𝑠𝑏_𝑘¹¹=𝜎_𝑘,𝑠𝑝_𝑘¹¹=1,𝑑𝑔_𝑘¹¹=−1−𝜎_𝑘,𝑠𝑏_𝑘¹²=ℎ²_𝑘𝑟12^𝑘,𝑠𝑝_𝑘¹²=ℎ²_𝑘𝑝12^𝑘,𝑑𝑔_𝑘¹²=ℎ²_𝑘𝑞12^𝑘,𝑅_𝑘¹¹=0,𝑠𝑏_𝑘²¹=ℎ_𝑘121+𝜎_𝑘×𝛾ℎ_𝑘𝑞_𝑘𝑎_𝑘𝜎_𝑘𝑐_𝑘+1+2+𝜎_𝑘𝑎_𝑘−1+𝑐_𝑘−1−𝜎_𝑘𝑞_𝑘𝑐_𝑘−𝑝_𝑘𝑐_𝑘+1+2+𝜎_𝑘𝑟_𝑘𝑎_𝑘−1+𝑐_𝑘−1,𝑠𝑝_𝑘²¹=ℎ_𝑘12𝜎_𝑘1+𝜎_𝑘×𝛾ℎ_𝑘𝑞_𝑘𝑎_𝑘2𝜎_𝑘+1𝑐_𝑘+1+𝑐_𝑘−1+2𝜎_𝑘+1𝑐_𝑘+1𝑝_𝑘+𝑐_𝑘𝑞_𝑘−𝑟_𝑘𝑐_𝑘−1,𝑑𝑔_𝑘²¹=ℎ_𝑘12𝜎_𝑘𝜎_𝑘𝑞_𝑘𝑐_𝑘−𝑝_𝑘𝑐_𝑘+1+𝑟_𝑘𝑐_𝑘−1−𝑐_𝑘+1𝑝_𝑘+𝑟_𝑘𝑐_𝑘−1−𝑐_𝑘𝑞_𝑘−𝛾ℎ_𝑘𝑞_𝑘𝑎_𝑘1+𝜎_𝑘𝑐_𝑘+1+𝑐_𝑘−1,𝑠𝑏_𝑘²²=𝜎_𝑘−𝛾ℎ³_𝑘𝑞12^𝑘𝑎_𝑘𝑏_𝑘−1+𝑐_𝑘+ℎ_𝑘121+𝜎_𝑘ℎ_𝑘𝛾𝑞_𝑘𝑎_𝑘𝑎_𝑘+1𝜎_𝑘+1+𝜎_𝑘𝑟_𝑘𝑏_𝑘−1+𝜎_𝑘𝑝_𝑘𝑎_𝑘+1−𝑎_𝑘𝑞_𝑘,𝑠𝑝_𝑘²²=1+𝛾ℎ³_𝑘𝑞12^𝑘𝑎_𝑘𝑏_𝑘+1+𝑐_𝑘+ℎ_𝑘12𝜎_𝑘1+𝜎_𝑘×ℎ_𝑘𝛾𝑞_𝑘𝑎_𝑘2𝜎_𝑘+1𝑎_𝑘+1+𝑎_𝑘−1+𝜎_𝑘𝜎_𝑘+1𝑝_𝑘𝑏_𝑘+1+1+2𝜎_𝑘𝑝_𝑘𝑎_𝑘+1−𝑟_𝑘𝑎_𝑘−1+𝑎_𝑘𝑞_𝑘,𝑑𝑔_𝑘²²=−1−𝜎_𝑘+ℎ_𝑘12𝜎_𝑘𝜎_𝑘𝑞_𝑘𝑎_𝑘−𝑝_𝑘𝑎_𝑘+1+𝑟_𝑘𝑎_𝑘−1−𝑝_𝑘𝑎_𝑘+1+𝑟_𝑘𝑎_𝑘−1−𝑎_𝑘𝑞_𝑘−ℎ_𝑘𝑞_𝑘𝛾𝑎_𝑘𝜎_𝑘+1𝑎_𝑘+1+𝑎_𝑘−1−𝑏_𝑘𝜎_𝑘,𝑅_𝑘²¹=−𝛾ℎ³_𝑘𝑎12^𝑘𝑞_𝑘𝑔_𝑘+1−𝑔_𝑘−1−ℎ²_𝑘𝑞12^𝑘𝑔_𝑘+𝑝_𝑘𝑔_𝑘+1+𝑟_𝑘𝑔_𝑘−1. (22)

Note that the scheme (21) fails when the solution is to be determined at 𝑘=1. We overcome this difficulty by modifying the scheme in such a way that the solutions retain order and accuracy even in the vicinity of singularity 𝑟=0. We consider the following approximations: 𝑎_𝑘+1=𝑎_𝑘+ℎ_𝑘𝜎_𝑘𝑎_𝑘⁽¹⁾+ℎ²_𝑘𝜎²_𝑘𝑎2!^𝑘(2)+ℎ³_𝑘𝜎³_𝑘𝑎3!^𝑘(3)+ℎ⁴_𝑘𝜎⁴_𝑘𝑎4!^𝑘(4)+𝑂ℎ⁵_𝑘,𝑎_𝑘−1=𝑎_𝑘−ℎ_𝑘𝑎_𝑘⁽¹⁾+ℎ²_𝑘𝑎2!^𝑘(2)−ℎ³_𝑘𝑎3!^𝑘(3)+ℎ⁴_𝑘𝑎4!^𝑘(4)+𝑂ℎ⁵_𝑘. (23)

Using the similar approximations of 𝑏_𝑘±1, 𝑐_𝑘±1 and 𝑔_𝑘±1 and neglecting 𝑂(ℎ⁵_𝑘) terms, we can rewrite (21) with the help of compact operators as 𝑠𝑏_𝑘¹¹𝑠𝑏_𝑘¹²𝑠𝑏_𝑘²¹𝑠𝑏_𝑘²²𝑈_𝑘−1𝑉_𝑘−1+𝑑𝑔_𝑘¹¹𝑑𝑔_𝑘¹²𝑑𝑔_𝑘²¹𝑑𝑔_𝑘²²𝑈_𝑘𝑉_𝑘+𝑠𝑝_𝑘¹¹𝑠𝑝_𝑘¹²𝑠𝑝_𝑘²¹𝑠𝑝_𝑘²²𝑈_𝑘+1𝑉_𝑘+1=𝑅_𝑘¹¹𝑅_𝑘²¹, (24)

where 𝑠𝑏_𝑘²¹=12𝛾𝑞_𝑘𝑎_𝑘𝜎_𝑘𝜎_𝑘+1²ℎ²_𝑘+6𝜎_𝑘+1𝜎_𝑘×−𝑟_𝑘𝜎_𝑘+𝜎_𝑘𝑝_𝑘−𝑞_𝑘𝜎_𝑘−2𝑟_𝑘ℎ_𝑘𝑐_𝑘+−6𝛾𝑞_𝑘𝑎_𝑘𝜎_𝑘𝜎_𝑘+1²ℎ³_𝑘+6𝜎_𝑘+1×𝜎_𝑘𝑝_𝑘𝜎²_𝑘+𝑟_𝑘𝜎_𝑘+2𝑟_𝑘ℎ²_𝑘𝑐_𝑘⁽¹⁾−3𝜎_𝑘+1×𝜎_𝑘𝑝_𝑘𝜎²_𝑘+𝑟_𝑘ℎ³_𝑘𝑐_𝑘⁽²⁾−3𝜎_𝑘𝜎_𝑘×+22𝛾𝑞_𝑘𝑎_𝑘𝑎_𝑘⁽²⁾+𝑎_𝑘⁽²⁾𝑟_𝑘ℎ³_𝑘+6𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1𝑟_𝑘𝑎_𝑘⁽¹⁾+𝛾𝑎²_𝑘𝑞_𝑘ℎ²_𝑘−6𝑟_𝑘𝑎_𝑘𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1ℎ_𝑘×72𝜎_𝑘(𝜎_𝑘+1)²⁻¹,𝑠𝑝_𝑘²¹=6𝑞_𝑘𝑎_𝑘𝛾𝜎_𝑘𝜎_𝑘+2ℎ³_𝑘𝑎_𝑘⁽¹⁾+12𝑎_𝑘𝛾𝑞_𝑘𝜎_𝑘+1²ℎ²_𝑘+6𝜎_𝑘×+12𝜎_𝑘+1𝑝_𝑘−𝑟_𝑘+𝑞_𝑘ℎ_𝑘𝑐_𝑘+6ℎ²_𝑘𝜎_𝑘+1𝑞_𝑘𝑎_𝑘𝛾𝜎_𝑘𝜎_𝑘+1ℎ_𝑘+𝜎_𝑘2𝜎_𝑘+1𝑝_𝑘+𝑟_𝑘𝑐_𝑘⁽¹⁾+3𝜎_𝑘+1𝜎_𝑘𝑝_𝑘𝜎²_𝑘+𝑟_𝑘ℎ³_𝑘𝑐_𝑘⁽²⁾+3𝑟_𝑘𝑎_𝑘⁽²⁾𝜎_𝑘×𝜎_𝑘+2ℎ³_𝑘×72𝜎_𝑘(1+𝜎_𝑘)²⁻¹,𝑑𝑔_𝑘²¹=𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1𝑎_𝑘⁽³⁾𝑟_𝑘+3𝛾𝑞_𝑘𝑎_𝑘𝑎_𝑘⁽²⁾ℎ⁴_𝑘−3𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+12𝛾𝑞_𝑘𝑎_𝑘𝑎_𝑘⁽¹⁾+𝑎_𝑘⁽²⁾𝑟_𝑘ℎ³_𝑘−6𝜎_𝑘+1³𝑝_𝑘𝜎_𝑘𝑐_𝑘⁽¹⁾+2𝛾𝑞_𝑘𝑎_𝑘𝑐_𝑘+𝑟_𝑘𝑐_𝑘⁽¹⁾ℎ²_𝑘−6𝑐_𝑘𝜎_𝑘+1²−𝑟_𝑘𝜎_𝑘+𝜎_𝑘𝑝_𝑘−𝑞_𝑘𝜎_𝑘−𝑟_𝑘+𝑞_𝑘+𝑝_𝑘×ℎ_𝑘×72𝜎_𝑘(𝜎_𝑘+1)²⁻¹,𝑠𝑏_𝑘²²=1+6𝛾𝑞_𝑘𝜎²_𝑘𝜎_𝑘+1ℎ²_𝑘𝑎²_𝑘−6ℎ_𝑘𝜎_𝑘×𝛾ℎ²_𝑘𝑞_𝑘𝜎²_𝑘−1+𝜎_𝑘𝑎_𝑘⁽¹⁾+𝑏_𝑘𝜎_𝑘+1²−𝜎_𝑘𝜎_𝑘+1𝑝_𝑘−𝑞_𝑘𝑎_𝑘+6𝑝_𝑘𝜎³_𝑘×𝜎_𝑘+1ℎ²_𝑘𝑎_𝑘⁽¹⁾−3𝜎_𝑘𝑝_𝑘𝜎³_𝑘+𝑝_𝑘𝜎²_𝑘−𝑟_𝑘×ℎ³_𝑘𝑎_𝑘⁽²⁾−6𝜎_𝑘𝜎_𝑘+1𝑝_𝑘𝜎²_𝑘+𝑟_𝑘ℎ³_𝑘𝑏_𝑘⁽¹⁾+6ℎ²_𝑘𝜎_𝑘𝜎_𝑘+1²𝑟_𝑘𝑏_𝑘−𝛾𝑞_𝑘𝑐_𝑘ℎ_𝑘×72𝜎²_𝑘(𝜎_𝑘+1)²⁻¹,𝑠𝑝_𝑘²²=1+12𝑞_𝑘𝛾𝜎_𝑘+1²ℎ²_𝑘𝑎²_𝑘+6𝑞_𝑘𝛾𝜎_𝑘𝜎²_𝑘−1+𝜎_𝑘ℎ³_𝑘𝑎_𝑘⁽¹⁾+6𝑞_𝑘𝛾𝑏_𝑘𝜎_𝑘𝜎_𝑘+1²ℎ³_𝑘+6𝜎_𝑘+12𝜎_𝑘𝑝_𝑘−𝑟_𝑘+𝑞_𝑘+𝑝_𝑘ℎ_𝑘𝑎_𝑘+6𝜎_𝑘+12𝑝_𝑘𝜎²_𝑘+𝜎_𝑘𝑝_𝑘+𝑟_𝑘×ℎ²_𝑘𝑎_𝑘⁽¹⁾+3𝜎_𝑘𝑝_𝑘𝜎³_𝑘+𝑝_𝑘𝜎²_𝑘−𝑟_𝑘ℎ³_𝑘𝑎_𝑘⁽²⁾+6𝜎_𝑘+1𝜎_𝑘𝑝_𝑘𝜎²_𝑘+𝑟_𝑘ℎ³_𝑘𝑏_𝑘⁽¹⁾+6𝑞_𝑘𝑐_𝑘𝛾𝜎_𝑘𝜎_𝑘+1²ℎ³_𝑘+6𝑝_𝑘𝜎_𝑘𝑏_𝑘×𝜎_𝑘+1²ℎ²_𝑘×72𝜎_𝑘(𝜎_𝑘+1)²⁻¹,𝑑𝑔_𝑘²²=−1−𝜎_𝑘+−12𝛾𝜎_𝑘+1³𝑞_𝑘ℎ²_𝑘𝑎²_𝑘+6𝛾𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1𝑞_𝑘ℎ³_𝑘𝑎_𝑘⁽¹⁾−3𝛾𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1𝑞_𝑘ℎ⁴_𝑘𝑎_𝑘⁽²⁾+6𝛾𝜎_𝑘𝜎_𝑘+1³𝑞_𝑘ℎ⁴_𝑘𝑏_𝑘⁽¹⁾+6𝜎_𝑘−1𝜎_𝑘+1²𝑞_𝑘+6𝜎_𝑘+1³×𝑟_𝑘−𝑝_𝑘ℎ_𝑘𝑎_𝑘−6𝜎_𝑘+1³×𝜎_𝑘𝑝_𝑘+𝑟_𝑘ℎ²_𝑘𝑎_𝑘⁽¹⁾+3𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1×𝑟_𝑘ℎ³_𝑘𝑎_𝑘⁽²⁾+6𝜎_𝑘+1²𝜎_𝑘𝜎_𝑘𝑝_𝑘−𝑟_𝑘ℎ³_𝑘𝑏_𝑘⁽¹⁾+3𝜎_𝑘𝜎_𝑘+1²𝜎²_𝑘𝑝_𝑘+𝑟_𝑘ℎ⁴_𝑘𝑏_𝑘⁽²⁾−𝑎_𝑘⁽³⁾𝜎_𝑘𝜎_𝑘+2𝜎_𝑘+1𝑟_𝑘ℎ⁴_𝑘+6𝜎_𝑘𝑏_𝑘𝜎_𝑘+1²𝑞_𝑘ℎ²_𝑘×72𝜎_𝑘𝜎_𝑘+1²⁻¹,𝑅_𝑘²¹ℎ=−^4𝑘2𝜎24^𝑘+1𝛾𝑞_𝑘𝑎_𝑘𝑔_𝑘⁽¹⁾+𝜎²_𝑘𝑝_𝑘+𝑟_𝑘𝑔_𝑘⁽²⁾−ℎ³_𝑘𝜎12^𝑘𝑝_𝑘−𝑟_𝑘𝑔_𝑘⁽¹⁾−ℎ²_𝑘𝑞12^𝑘+𝑟_𝑘+𝑝_𝑘𝑔_𝑘. (25)

The modified scheme (24) is free from the terms 1/(𝑘±1), hence easily solved for 𝑘=1(1)𝑛. The difference equation (24) along with the boundary conditions (4) gives a 2𝑛×2𝑛 linear system of equations for the unknowns 𝑈_𝑘, 𝑉_𝑘, 𝑘=1(1)𝑛. The resulting block tridiagonal system can be easily solved using block Gauss-Seidel algorithm.

4. Extension to Sixth-Order Differential Equations

The proposed method can be easily extended to the sixth-order differential equations: 𝑢⁽⁶⁾(𝑟)=𝑓𝑟,𝑢(𝑟),𝑢⁽¹⁾(𝑟),𝑢⁽²⁾(𝑟),𝑢⁽³⁾(𝑟),𝑢⁽⁴⁾(𝑟),𝑢⁽⁵⁾(𝑟),𝛼<𝑟<𝛽 (26) subject to the necessary boundary conditions: 𝑢(𝛼)=𝛼₀,𝑢⁽²⁾(𝛼)=𝛼₁,𝑢⁽⁴⁾(𝛼)=𝛼₂,𝑢𝛽=𝛽₀,𝑢⁽²⁾𝛽=𝛽₁,𝑢⁽⁴⁾𝛽=𝛽₂ (27)

or equivalently, 𝑢⁽²⁾(𝑟)(𝑟),𝑣=𝑣₍₂₎(𝑟)(𝑟),𝑤=𝑤₍₂₎(𝑟)=𝑓𝑟,𝑢,𝑢⁽¹⁾(𝑟),𝑣(𝑟),𝑣⁽¹⁾(𝑟),𝑤(𝑟),𝑤⁽¹⁾(𝑟),𝛼<𝑟<𝛽, (28) subject to the natural boundary conditions: 𝑢(𝛼)=𝛼₀,𝑣(𝛼)=𝛼₁,𝑤(𝛼)=𝛼₂,𝑢𝛽=𝛽₀,𝑣𝛽=𝛽₁,𝑤𝛽=𝛽₂. (29)

We outline the similar algorithm for (28) as follows: 𝑈^′_𝑘=𝑈_𝑘+1−1−𝜎²_𝑘𝑈_𝑘−𝜎²_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘=𝑉_𝑘+1−1−𝜎²_𝑘𝑉_𝑘−𝜎²_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑊^′_𝑘=𝑊_𝑘+1−1−𝜎²_𝑘𝑊_𝑘−𝜎²_𝑘𝑊_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑈^′_𝑘+1=1+2𝜎_𝑘𝑈_𝑘+1−1+𝜎_𝑘²𝑈_𝑘+𝜎²_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘+1=1+2𝜎_𝑘𝑉_𝑘+1−1+𝜎_𝑘²𝑉_𝑘+𝜎²_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑊^′_𝑘+1=1+2𝜎_𝑘𝑊_𝑘+1−1+𝜎_𝑘²𝑊_𝑘+𝜎²_𝑘𝑊_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑈^′_𝑘−1=−𝑈_𝑘+1+1+𝜎_𝑘²𝑈_𝑘−𝜎_𝑘2+𝜎_𝑘𝑈_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑉^′_𝑘−1=−𝑉_𝑘+1+1+𝜎_𝑘²𝑉_𝑘−𝜎_𝑘2+𝜎_𝑘𝑉_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑊^′_𝑘−1=−𝑊_𝑘+1+1+𝜎_𝑘²𝑊_𝑘−𝜎_𝑘2+𝜎_𝑘𝑊_𝑘−1ℎ_𝑘𝜎_𝑘1+𝜎_𝑘,𝑓_𝑘±1=𝑓𝑟_𝑘±1,𝑈_𝑘±1,𝑈^′_𝑘±1,𝑉_𝑘±1,𝑉^′_𝑘±1,𝑊_𝑘±1,𝑊^′_𝑘±1,𝑈^′_𝑘=𝑈^′_𝑘+𝛾ℎ_𝑘𝑉_𝑘+1−𝑉_𝑘−1,𝑉^′_𝑘=𝑉^′_𝑘+𝛾ℎ_𝑘𝑊_𝑘+1−𝑊_𝑘−1,𝑊^′_𝑘=𝑊^′_𝑘+𝛾ℎ_𝑘𝑓_𝑘+1−𝑓_𝑘−1, (30) where the values of 𝛾 are same as obtained in Section 2.

Then, the 𝑂(ℎ³_𝑘)-approximations for (26) or (28) can be obtained by the following relations for 𝑘=1(1)𝑛: 𝑈_𝑘+1−1+𝜎_𝑘𝑈_𝑘+𝜎_𝑘𝑈_𝑘−1+ℎ²𝑝12^𝑘𝑉_𝑘+1+𝑞_𝑘𝑉_𝑘+𝑟_𝑘𝑉_𝑘−1=𝑂ℎ⁵_𝑘,𝑉_𝑘+1−1+𝜎_𝑘𝑉_𝑘+𝜎_𝑘𝑉_𝑘−1+ℎ²𝑝12^𝑘𝑊_𝑘+1+𝑞_𝑘𝑊_𝑘+𝑟_𝑘𝑊_𝑘−1=𝑂ℎ⁵_𝑘,𝑊_𝑘+1−1+𝜎_𝑘𝑊_𝑘+𝜎_𝑘𝑊_𝑘−1+ℎ²_𝑘𝑝12^𝑘𝑓_𝑘+𝑞_𝑘𝑓_𝑘+𝑟_𝑘𝑓_𝑘−1=𝑂ℎ⁵_𝑘. (31)

The boundary conditions (27) are used to obtain values at 𝑘±1 for 𝑘=1 and 𝑛, respectively. The numerical scheme may be implemented by neglecting 𝑂(ℎ⁵_𝑘) terms from (31).

5. Convergence Analysis

In this section, we derive the difference scheme of singular problem and investigate the convergence property of the proposed scheme. Consider the model problem 𝑢⁽⁴⁾(𝑟)=𝑎(𝑟)𝑢(𝑟)+𝑔(𝑟),𝛼<𝑟<𝛽 (32) or equivalently 𝑢⁽²⁾(𝑟)(𝑟),𝑣=𝑣₍₂₎(𝑟)=𝑎(𝑟)𝑢(𝑟)+𝑔(𝑟),𝛼<𝑟<𝛽 (33)

along with the boundary conditions (2).

If 𝑎(𝑟),𝑔(𝑟)∈𝐶²[𝛼,𝛽] and 𝑎(𝑟)>0 on [𝛼,𝛽], then the boundary value problem (33) has a unique solution. Under these assumptions (see, [27]), there are constants 𝑎>0 and 𝑔 such that 𝑎=min_𝑘{𝑎_𝑘,𝑎_𝑘⁽¹⁾,𝑎_𝑘⁽²⁾}, 𝑔=min_𝑘{𝑔_𝑘,𝑔_𝑘⁽¹⁾,𝑔_𝑘⁽²⁾}.

For the convergence, the coefficients 𝑝_𝑘, 𝑞_𝑘, and 𝑟_𝑘 associated with (5) and (15) must be negative (see, [22]), from which we obtain the condition |𝜎_𝑘−√5/2|<1/2.

Now applying the methods (5) and (15) to (33) and using the similar technique discussed in Section 3 for singular coefficients 𝑎(𝑟) and 𝑔(𝑟), we obtain the following system of difference equations: 𝐏_𝑘𝐙_𝑘−1+𝐐_𝑘𝐙_𝑘+𝐑_𝑘𝐙_𝑘+1=𝐒_𝑘+𝐓_𝑘ℎ_𝑘,𝑘=1(1)𝑛, (34) where 𝐏_𝑘=−𝜎_𝑘ℎ²_𝑘𝜎12^𝑘1+𝜎_𝑘−𝜎²_𝑘ℎ²_𝑘𝜎12^𝑘1+𝜎_𝑘−𝜎²_𝑘×𝑎_𝑘+ℎ_𝑘𝑎_𝑘⁽¹⁾−𝜎_𝑘,𝐐_𝑘=𝜎_𝑘×1ℎ+1_2𝑘𝜎12^2𝑘+3𝜎_𝑘ℎ+1_2𝑘2𝜎24^2𝑘+3𝜎_𝑘+1𝑎_𝑘+4ℎ_𝑘𝜎_𝑘𝜎_𝑘−1𝑎_𝑘⁽¹⁾+ℎ²_𝑘𝜎_𝑘𝜎²_𝑘−𝜎_𝑘+1𝑎_𝑘⁽²⁾1,𝐑_𝑘=ℎ−1^2𝑘𝜎12^2𝑘+𝜎_𝑘ℎ−1_2𝑘𝜎12^2𝑘+𝜎_𝑘−1𝑎_𝑘+ℎ_𝑘𝜎_𝑘𝜎²_𝑘−𝜎_𝑘+1𝑎_𝑘⁽¹⁾,𝐒−1_𝑘=ℎ0,^2𝑘𝜎12^3𝑘−𝜎²_𝑘−𝜎_𝑘𝑔_𝑘⁽²⁾−𝜎²_𝑘+𝜎_𝑘−1𝑔_𝑘⁽¹⁾−𝜎³_𝑘+4𝜎²_𝑘+4𝜎_𝑘+1𝑔_𝑘^𝑇,𝐙_𝑘=𝑈_𝑘,𝑉_𝑘^𝑇,𝐓_𝑘ℎ_𝑘=𝑂ℎ⁵_𝑘. (35) Incorporating the boundary values 𝑈₀=𝛼₀, 𝑉₀=𝛼₁, 𝑈_𝑛+1=𝛽₀, and 𝑉_𝑛+1=𝛽₁, the system of difference equations (34) in the matrix-vector form can be written as 𝐌𝐙=𝐉+𝐓_𝑘ℎ_𝑘, (36) where 𝐌=[𝐏_𝑘𝐐_𝑘𝐑_𝑘] is the block tridiagonal matrix, 𝐉=[𝐒₁−𝐏₁𝜶,𝐒₂,…,𝐒_𝑛−1,𝐒_𝑛−𝐑_𝑛𝜷]^𝑇, 𝜶=[𝛼₀,𝛼₁]^𝑇, and 𝜷=[𝛽₀,𝛽₁]^𝑇.

Let 𝑧_𝑘=[𝑢_𝑘,𝑣_𝑘]^𝑇, 𝑘=1(1)𝑛, and 𝐳=[𝑧₁,…,𝑧_𝑛]≅𝐙, which satisfies 𝐌𝐳=𝐉. (37)

Let 𝜺=[𝜀₁,𝜀₂,…,𝜀_𝑛]^𝑇 be the discretization error vector, and let 𝜀_𝑘=𝑍_𝑘−𝑧_𝑘, 𝑘=1(1)𝑛 be the discretization errors at the node 𝑟_𝑘. Subtracting (37) from (36), we obtain the error equation: 𝐌𝜺=𝐓_𝑘ℎ_𝑘. (38) Also, we obtain 𝐏_𝑘∞=max_{2≤𝑘≤𝑛}𝜎_𝑘+ℎ²_𝑘𝜎12^𝑘1+𝜎_𝑘−𝜎²_𝑘,𝜎_𝑘+ℎ²_𝑘𝜎12^𝑘1+𝜎_𝑘−𝜎²_𝑘⋅𝑎_𝑘+ℎ_𝑘𝑎_𝑘⁽¹⁾,𝐑_𝑘∞=max_{1≤𝑘≤𝑛−1}ℎ1+^2𝑘𝜎12^2𝑘+𝜎_𝑘,ℎ−11+^2𝑘𝜎12^2𝑘+𝜎_𝑘−1𝑎_𝑘+ℎ_𝑘𝜎_𝑘𝜎²_𝑘−𝜎_𝑘+1𝑎_𝑘⁽¹⁾. (39)

Thus for sufficiently small ℎ_𝑘 or equivalently as ℎ_𝑘→0, we obtain the relations ‖𝐏_𝑘‖_∞=𝜎_𝑘, 𝑘=2(1)𝑛 and ‖𝐑_𝑘‖_∞=1, 𝑘=1(1)𝑛−1. Hence, the graph 𝐆(𝐌) of the matrix 𝐌 is strongly connected, and thus the matrix 𝐌 is irreducible (see, [28]).

Further, let Σ_𝑙 be the sum of the 𝑙th row sum of the matrix 𝐌; then we have the following.

For 𝑙=1, Σ_𝑙=𝜎_𝑙+ℎ²_𝑙𝜎12^𝑙𝜎²_𝑙+5𝜎_𝑙,Σ+5_𝑙+1=𝜎_𝑙+ℎ²_𝑙𝑎12^𝑙𝜎_𝑙𝜎²_𝑙+5𝜎_𝑙+ℎ+5_3𝑙𝑎12^𝑙(1)𝜎_𝑙3𝜎²_𝑙−𝜎_𝑙+ℎ−1^4𝑙𝑎24^𝑙(2)𝜎_𝑙𝜎³_𝑙+1. (40)

For 𝑙=3(2)2𝑛−3, Σ_𝑙=ℎ²_𝑙2𝜎_𝑙𝜎_𝑙,Σ+1_𝑙+1=ℎ²_𝑙2𝑎_𝑙𝜎_𝑙𝜎_𝑙+ℎ+1^3𝑙6𝑎_𝑙⁽¹⁾𝜎_𝑙𝜎²_𝑙+ℎ−1_4𝑙𝑎24^𝑙(2)𝜎_𝑙𝜎³_𝑙+1. (41)

For 𝑙=2𝑛−1, Σ_𝑙ℎ=1+^2𝑙125𝜎²_𝑙+5𝜎_𝑙,Σ+1_𝑙+1ℎ=1+^2𝑙𝑎12^𝑙5𝜎²_𝑙+5𝜎_𝑙+ℎ+1^3𝑙𝑎12^𝑙(1)𝜎_𝑙𝜎²_𝑙+𝜎_𝑙+ℎ−3_4𝑙𝑎24^𝑙(2)𝜎_𝑙𝜎³_𝑙+1. (42)

This implies the following.

For 𝑙=1, Σ_𝑙=𝜎_𝑙+ℎ²_𝑙𝜎12^𝑙𝜎²_𝑙+5𝜎_𝑙,Σ+5_𝑙+1≥𝜎_𝑙+ℎ²_𝑙𝑎12𝜎_𝑙𝜎²_𝑙+5𝜎_𝑙+ℎ+5_3𝑙𝑎12𝜎_𝑙3𝜎²_𝑙−𝜎_𝑙+ℎ−1^4𝑙𝑎24𝜎_𝑙𝜎³_𝑙+1. (43)

For 𝑙=3(2)2𝑛−3, Σ_𝑙=ℎ²_𝑙2𝜎_𝑙𝜎_𝑙,Σ+1_𝑙+1≥ℎ²_𝑙2𝑎𝜎_𝑙𝜎_𝑙+ℎ+1^3𝑙6𝑎𝜎_𝑙𝜎²_𝑙+ℎ−1_4𝑙𝑎24𝜎_𝑙𝜎³_𝑙+1. (44)

For 𝑙=2𝑛−1, Σ_𝑙ℎ=1+^2𝑙125𝜎²_𝑙+5𝜎_𝑙,Σ+1_𝑙+1ℎ≥1+^2𝑙𝑎125𝜎²_𝑙+5𝜎_𝑙+ℎ+1_3𝑙𝑎12𝜎_𝑙𝜎²_𝑙+𝜎_𝑙+ℎ−3^4𝑙𝑎24𝜎_𝑙𝜎³_𝑙+1, (45)

For sufficiently small value of ℎ_𝑘, that is, in the limiting case as ℎ_𝑘→0, we obtain Σ_𝑙≥𝜎_𝑙Σ>0,𝑙=1,2,_𝑙(1)Σ≥0,𝑙=32𝑛−2,_𝑙≥1>0,𝑙=2𝑛−1,2𝑛. (46)

Hence we find that 𝐌 is monotone (see, [29, 30]). Consequently 𝐌⁻¹ exists, and 𝐌⁻¹≥0.

Let 𝐌⁻¹_𝑖,𝑙 be the (𝑖, 𝑙)th element of 𝐌⁻¹, and we define 𝐌⁻¹=max_{1≤𝑖≤2𝑛}^2𝑛_𝑙=1𝐌⁻¹_𝑖,𝑙,‖𝐓‖=max_{1≤𝑙≤2𝑛}^2𝑛_𝑙=1𝐓_𝑙ℎ_𝑙=𝑂ℎ³_𝑙. (47)

From the theory of matrix, we know that ^2𝑛_𝑙=1𝐌⁻¹_𝑖,𝑙Σ_𝑙=1,1≤𝑖≤2𝑛. (48)

Thus the following bounds can be estimated with the help of series expansions.

For 𝑙=1, 𝐌⁻¹_𝑖,𝑙≤Σ_𝑙⁻¹=1𝜎_𝑙−ℎ²_𝑙12𝜎_𝑙𝜎²_𝑙+5𝜎_𝑙+5+𝑂ℎ⁴_𝑙,𝐌⁻¹_𝑖,𝑙+1≤Σ⁻¹_𝑙+1≤1𝜎_𝑙−ℎ²_𝑙𝑎12𝜎_𝑙𝜎²_𝑙+5𝜎_𝑙−ℎ+5_3𝑙𝑎12𝜎_𝑙3𝜎²_𝑙−𝜎_𝑙−1+𝑂ℎ⁴_𝑙. (49)

For 𝑙=3(2)2𝑛−3, 𝐌⁻¹_𝑖,𝑙≤min_𝑙Σ_𝑙⁻¹=2ℎ²_𝑙𝜎_𝑙𝜎_𝑙,𝐌+1_{−1𝑖,𝑙+1}≤min_𝑙Σ⁻¹_𝑙+1≤2ℎ²_𝑙𝑎𝜎_𝑙𝜎_𝑙−2𝜎+1_𝑙−13ℎ_𝑙𝑎𝜎_𝑙𝜎_𝑙+𝜎+1^2𝑙−5𝜎_𝑙+118𝑎𝜎_𝑙𝜎_𝑙ℎ+1+𝑂_𝑙. (50)

For 𝑙=2𝑛−1, 𝐌⁻¹_𝑖,𝑙≤Σ_𝑙⁻¹ℎ=1−^2𝑙125𝜎²_𝑙+5𝜎_𝑙+1+𝑂ℎ⁴_𝑙,𝐌⁻¹_𝑖,𝑙+1≤Σ⁻¹_𝑙+1ℎ≤1−^2𝑙𝑎125𝜎²_𝑙+5𝜎_𝑙−ℎ+1_3𝑙𝑎𝜎12^𝑙𝜎²_𝑙+𝜎_𝑙−3+𝑂ℎ⁴_𝑙. (51)

With the help of (48), we obtain the following bounds: 𝐌⁻¹≤2𝑎+1ℎ²_𝑙𝑎𝜎_𝑙𝜎_𝑙−2𝜎+1^𝑙−13ℎ_𝑙𝑎𝜎_𝑙𝜎_𝑙+𝜎+1_2𝑙−5𝜎_𝑙+1+36𝑎𝜎_𝑙+1²18𝑎𝜎_𝑙𝜎_𝑙ℎ+1+𝑂_𝑙. (52)

From (38) and (52), we obtain the following error estimates: ‖𝜺‖≤𝐌⁻¹⋅𝐓_𝑙ℎ_𝑙≤𝑂ℎ³_𝑙. (53)

This proves the third-order convergence of the proposed method. We generalize the above results in the following theorem.

Theorem 1. The method given by (15) for the numerical solution of fourth-order singular differential equation (1) with sufficiently small ℎ_𝑘 and |𝜎_𝑘−√5/2|<1/2, 𝜎_𝑘≠1 gives a third-order convergent solution.

6. Computational Illustrations

To illustrate the geometric mesh finite difference method, we have solved both linear and nonlinear problems. The boundary conditions may be obtained from the analytical solution as a test procedure. The numerical accuracy of results are tested using maximum absolute errors and root mean square errors with the error tolerance being ≤10⁻¹⁵. For the simplicity in computation, we choose 𝜎_𝑘=𝜎=constant, for 𝑘=1(1)𝑛 and define the geometric mesh as follows ([24]): 𝑥₀ℎ=𝛼,₁=𝛽−𝛼(1−𝜎)1−𝜎^𝑛+1(),𝜎<1,𝛽−𝛼𝜎−1𝜎^𝑛+1−1,𝜎>1. (54)

The subsequent mesh spacing is determined by ℎ_𝑘+1=𝜎ℎ_𝑘, 𝑘=1(1)𝑛. If the boundary value problems exhibit layer behaviour near the left boundary (see, [21]), 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.

All the numerical computations are performed using long double length arithmetic in 𝐶 under Linux operating system with 2 GB operational memory.

Example 1. Consider the fourth-order linear problem (see, [31]) in the polar form: Δ⁴𝑢(𝑟)≡𝑑²𝑑𝑟²+𝜆𝑟𝑑𝑑𝑟²𝑢(𝑟)=1+2𝜆𝑟+𝜆(𝜆−2)𝑟²−𝜆(𝜆−2)𝑟³𝑒^𝑟,0<𝑟<1. (55)
The analytical solution is 𝑢(𝑟)=𝑒^𝑟. The errors estimates for various values of 𝑛 are reported in Tables 1 and 2 for uniform mesh (𝜎=1) and geometric mesh (𝜎≠1), respectively.

Table 1 

The maximum absolute errors for Example 1.

Table 2 

The root mean square errors for Example 1.

Example 2. Consider the boundary value problems that arise from time-dependent Navier-Stokes equation (see, [32]) for axis symmetric flow of an incompressible fluid contained between infinite disks 𝑢⁽⁴⁾(𝑟)=𝜆𝑢(𝑟)𝑢⁽²⁾(𝑟)−𝜆𝑟²×−11+4𝑟+𝑟²𝑒^2𝑟−11+8𝑟+𝑟²𝑒^𝑟,0<𝑟<1. (56)
The analytical solution is 𝑢(𝑟)=(1−𝑟²)𝑒^𝑟. The errors estimates are reported in Table 3 for various values of 𝑛 and 𝜆=10³.

Table 3 

The maximum absolute errors and root mean square errors for Example 2.

Example 3. Consider the sixth-order linear singular problem: Δ⁶𝑢(𝑟)≡𝑑²𝑑𝑟²+𝜆𝑟𝑑𝑑𝑟³𝑢(𝑟)=()1+3𝜆𝜆−2𝑟²+3𝜆(𝜆−2)(𝜆−4)𝑟⁴(𝑟)+𝜆sinh𝑟()3+𝜆−2)(𝜆−7𝑟²+3(𝜆−2)(𝜆−4)𝑟⁴cosh(𝑟),0<𝑟<1. (57)
The analytical solution is 𝑢(𝑟)=sinh(𝑟). The errors estimates for various values of 𝑛 are reported in Tables 4 and 5 for uniform mesh (𝜎=1) and geometric mesh (𝜎≠1), respectively.

Table 4 

The maximum absolute errors for Example 3.

Table 5 

The root mean square errors for Example 3.

Example 4. Consider the nonlinear problem 𝑢⁽⁶⁾(𝑟)=𝜆𝑢(𝑟)𝑢⁽²⁾(𝑟)+𝑢⁽²⁾(𝑟)𝑢⁽⁴⁾(𝑟)+𝑢(𝑟)𝑢⁽⁴⁾(𝑟)+sinh(𝑟)(1−3𝜆sinh(𝑟)),0<𝑟<1. (58)
The analytical solution is 𝑢(𝑟)=sinh(𝑟). The errors estimates are reported in Table 6 for various values of 𝑛 and 𝜆=10³.

Table 6 

The maximum absolute errors and root mean square errors for Example 4.

7. Conclusion

The numerical results confirm that the proposed geometric mesh finite difference scheme converges and applicable to both singular and nonsingular differential equations. The numerical accuracy obtained using geometric mesh shows superiority over corresponding uniform mesh. The optimum mesh ratio parameter within the specified convergent region may be obtained by simulations. We have employed block Gauss-Seidel method to solve the block matrix systems. The method can be extended to general even-order nonlinear differential equations. Application to the proposed scheme to nonlinear singular elliptic problems is an open problem.

Acknowledgment

The authors are thankful to the Professor Rüdiger Weiner and unknown referee for their valuable suggestions which improve the quality of paper.