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

Authors: Navnit Jha; R. K. Mohanty; Vinod Chauhan
Journal: Advances in Numerical Analysis
 (2013)

---

## 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.


---

## Body

## 1. Introduction

Consider the fourth-order boundary value problem:(1)u(4)(r)=f(r,u(r),u(1)(r),u(2)(r),u(3)(r)),hhhhhhhhhhhhhhhhhhhhhhhhhhα<r<β
subject to the necessary boundary conditions:
(2)u(α)=α0,u(2)(α)=α1,u(β)=β0,u(2)(β)=β0,wheref∈C4[α,β] and α0, α1, β0, and β1 are real constants and -∞<α<r<β<∞.Or equivalently(3)u(2)(r)=v(r),v(2)(r)=f(r,u(r),u(1)(r),v(r),v(1)(r)),hhhhhhhhhhhhhhhhhhhhhhhhα<r<βsubject to the natural boundary conditions:(4)u(α)=α0,v(α)=α1,u(β)=β0,v(β)=β1.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 Section2, 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 α=r0<r1<⋯<rn<rn+1=β. Let hk=rk-rk-1 and k=1(1)n+1 be the nonuniform step size, and let σk=hk+1/hk>0 be the geometric mesh ratio. Let Uk=u(rk) and Vk=v(rk) denote the exact solution values of  u and v at the mesh rk and uk, vk be their approximate solutions, respectively.Consider the following three-point geometric mesh discretizations foru(2)(r)=v(r) and v(2)(r)=f(r):
(5)Uk+1-(1+σk)Uk+σkUk-1+h212(pkVk+1+qkVk+rkVk-1)=O(hk5),k=1(1)n,(6)Vk+1-(1+σk)Vk+σkVk-1+h212(pkfk+1+qkfk+rkfk-1)=Ek,k=1(1)n,
where
(7)Ek=hk5360σk(2σk+1)(σk+2)(1-σk2)fk(5)(ξ)Ek=+O(hk6),rk-1<ξ<rk+1,pk=1-σk-σk2,qk=-(σk+1)(σk2+3σk+1),rk=σk(σk2-σk-1).Define(8)U~k′=Uk+1-(1-σk2)Uk-σk2Uk-1hkσk(1+σk),V~k′=Vk+1-(1-σk2)Vk-σk2Vk-1hkσk(1+σk),(9)U~k+1′=(1+2σk)Uk+1-(1+σk)2Uk+σk2Uk-1hkσk(1+σk),V~k+1′=(1+2σk)Vk+1-(1+σk)2Vk+σk2Vk-1hkσk(1+σk),U~k-1′=-Uk+1+(1+σk)2Uk-σk(2+σk)Uk-1hkσk(1+σk),V~k-1′=-Vk+1+(1+σk)2Vk-σk(2+σk)Vk-1hkσk(1+σk),f~k±1=f(rk±1,Uk±1,U~k±1′,Vk±1,V~k±1′).It is easy to verify that(10)f~k+1=fk+1-hk26σk(σk+1)(GkUk(3)+HkVk(3))+O(hk3),f~k-1=fk-1-hk26(σk+1)(GkUk(3)+HkVk(3))+O(hk3),
where  Gk=(∂f/∂U(1))rk, Hk=(∂f/∂V(1))rk, and so forth.Now, let(11)U^k′=U~k′+γhk(Vk+1-Vk-1),V^k′=V~k′+δhk(f~k+1-f~k-1),
where γ and δ are free parameters to be determined.With the help of (8) and (11), it follows that
(12)U^k′=Uk(1)+hk26(σk+6γ(1+σk))Uk(3)+O(hk3),V^k′=Vk(1)+hk26(σk+6δ(1+σk))Vk(3)+O(hk3).Further, we define(13)f^k=f(rk,Uk,U^k′,Vk,V^k′).With the help of (12)-(13), it follows that
(14)f^k=fk+hk26(σk+6γ(1+σk))GkUk(3)f^k=+hk26(σk+6δ(1+σk))HkVk(3)+O(hk3).With the help of approximations (10), (14), the difference equation (6) at each internal mesh rk, is approximated as follows:
(15)Vk+1-(1+σk)Vk+σkVk-1+hk212(pkf~k+qkf^k+rkf~k-1)=E~k,k=1(1)n,
where
(16)E~k=-hk472Uk(3)Gk(1+σk)(σk36γ+(1+24γ)σk(1+σk)+(1+6γ)σk3)-hk472Vk(3)Hk×(1+σk)(6δ+(1+24δ)σk+(1+6δ)σk3)+O(hk5),σk≠1.The difference scheme (15) to be of O(hk3), the coefficients of hk4 in (16) must be zero, and hence we obtain
(17)6γ+(1+24γ)σk(1+σk)+(1+6γ)σk3=0,6δ+(1+24δ)σk(1+σk)+(1+6δ)σk3=0.Thus the values associated with (17) are given by γ=δ=-σk(σk2+σk+1)/(6(1+σk)(σk2+3σk+1)), and the local truncation errors given by (16) becomes  E~k=O(hk5), σk≠1.  However, for σk≠1, the error reduces to E~k=O(hk6).


## 3. Application to Singular Problem

Consider the singular fourth-order linear differential equation in cylindrical polar coordinates:(18)Δ4u(r)≡(d2dr2+λrddr)2u(r)=g(r),0<r<1
or equivalently
(19)u(2)(r)=v(r),(20)v(2)(r)=a(r)v(1)(r)+b(r)v(r)+c(r)u(1)(r)+g(r),wherea(r)=-2λ/r,   b(r)=-λ(λ-2)/r2, and c(r)=λ(λ-2)/r3 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 k=1(1)n:
(21)[sbk11sbk12sbk21sbk22][Uk-1Vk-1]+[dgk11dgk12dgk21dgk22][UkVk]+[spk11spk12spk21spk22][Uk+1Vk+1]=[Rk11Rk21],
where
(22)sbk11=σk,spk11=1,dgk11=-1-σk,sbk12=hk212rk,spk12=hk212pk,dgk12=hk212qk,Rk11=0,sbk21=hk12(1+σk)×(r52γhkqkak(σkck+1+(2+σk)(ak-1+ck-1))-σk(qkck-pkck+1)+(2+σk)rk(ak-1+ck-1)),spk21=hk12σk(1+σk)×(γhkqkak((2σk+1)ck+1+ck-1)+(2σk+1)ck+1pk+ckqk-rkck-1),dgk21=hk12σk(-γhkqkak(1+σk)(ck+1+ck-1)σk(qkck-pkck+1+rkck-1)-ck+1pk+rkck-1-ckqk-γhkqkak(1+σk)(ck+1+ck-1)),sbk22=σk-γhk312qk(akbk-1+ck)+hk12(1+σk)(hk(γqkakak+1σk+(1+σk)rkbk-1)+σk(pkak+1-akqk)),spk22=1+γhk312qk(akbk+1+ck)+hk12σk(1+σk)×(f22hk(γqkak((2σk+1)ak+1+ak-1)+σk(σk+1)pkbk+1)+(1+2σk)pkak+1-rkak-1+akqkf22),dgk22=-1-σk+hk12σk((σk(qkak-pkak+1+rkak-1)-pkak+1+rkak-1-akqk)-hkqk(γak(σk+1)(ak+1+ak-1)-bkσk)),Rk21=-γhk312akqk(gk+1-gk-1)-hk212(qkgk+pkgk+1+rkgk-1).Note that the scheme (21) fails when the solution is to be determined at k=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 r=0. We consider the following approximations:
(23)ak+1=ak+hkσkak(1)+hk2σk22!ak(2)+hk3σk33!ak(3)+hk4σk44!ak(4)+O(hk5),ak-1=ak-hkak(1)+hk22!ak(2)-hk33!ak(3)+hk44!ak(4)+O(hk5).Using the similar approximations ofbk±1, ck±1 and gk±1 and neglecting O(hk5) terms, we can rewrite (21) with the help of compact operators as
(24)[sbk11sbk12sb~k21sb~k22][Uk-1Vk-1]+[dgk11dgk12dg~k21dg~k22][UkVk]+[spk11spk12sp~k21sp~k22][Uk+1Vk+1]=[Rk11R~k21],where(25)sb~k21=((12γqkakσk(σk+1)2hk2+6(σk+1)σk×(-rkσk+σkpk-qkσk-2rk)hkhk2)ck+(-6γqkakσk(σk+1)2hk3+6(σk+1)×σk(pkσk2+rkσk+2rk)hk2)ck(1)-3(σk+1)×σk(pkσk2+rk)hk3ck(2)-3σk(σk+2)×(2γqkakak(2)+ak(2)rk)hk3+6σk(σk+2)(σk+1)(rkak(1)+γak2qk)hk2-6rkakσk(σk+2)(σk+1)hkfv2)×(72σk(σk+1)2)-1,spk21=(6qkakγσk(σk+2)hk3ak(1)+(12akγqk(σk+1)2hk2+6(σk+1)×((2σk+1)pk-rk+qk)hkf2)ck+6hk2(σk+1)(qkakγσk(σk+1)hk+σk(2σk+1)pk+rk)ck(1)+3(σk+1)σk(pkσk2+rk)hk3ck(2)+3rkak(2)σk×(σk+2)hk3)×(72σk(1+σk)2)-1,dg~k21=(σk(σk+2)(σk+1)(ak(3)rk+3γqkakak(2))hk4-3σk(σk+2)(σk+1)(2γqkakak(1)+ak(2)rk)hk3-6(σk+1)3(pkσkck(1)+2γqkakck+rkck(1))hk2-6ck(σk+1)2(-rkσk+σkpk-qkσk-rk+qk+pk)×hkak(3))×(72σk(σk+1)2)-1,sb~k22=1+(6γqkσk2(σk+1)hk2ak2-6hkσk1+h×(γhk2qk((σk2-1+σk)ak(1)+bk(σk+1)2)-σk(σk+1)(pk-qk)f3)ak+6pkσk31+h×(σk+1)hk2ak(1)-3σk(pkσk3+pkσk2-rk)1+h×hk3ak(2)-6σk(σk+1)(pkσk2+rk)hk3bk(1)1+h+6hk2σk(σk+1)2(rkbk-γqkckhk))×(72σk2(σk+1)2)-1,sp~k22=1+(12qkγ(σk+1)2hk2ak2+(6qkγσk(σk2-1+σk)hk3ak(1)1+6qkγbkσk(σk+1)2hk31+6(σk+1)(2σkpk-rk+qk+pk)hkhk3)ak1+6(σk+1)(2pkσk2+σkpk+rk)×hk2ak(1)+3σk(pkσk3+pkσk2-rk)hk3ak(2)+6(σk+1)σk(pkσk2+rk)hk3bk(1)+6qkckγσk(σk+1)2hk3+6pkσkbk×(σk+1)2hk2)×(72σk(σk+1)2)-1,dg~k22=-1-σk+(g6-12γ(σk+1)3qkhk2ak2+(6γσk(σk+2)(σk+1)qkhk3ak(1)-3γσk(σk+2)(σk+1)qkhk4ak(2)+6γσk(σk+1)3qkhk4bk(1)+(6(σk-1)(σk+1)2qk+6(σk+1)3×(rk-pk)h63)hkh63)ak-6(σk+1)3×(σkpk+rk)hk2ak(1)+3σk(σk+2)(σk+1)×rkhk3ak(2)+6(σk+1)2σk(σkpk-rk)hk3bk(1)+3σk(σk+1)2(σk2pk+rk)hk4bk(2)-ak(3)σk(σk+2)(σk+1)rkhk4+6σkbk(σk+1)2qkhk2)×(72σk(σk+1)2)-1,R~k21=-hk424(2(σk+1)γqkakgk(1)+(σk2pk+rk)gk(2))-hk312(σkpk-rk)gk(1)-hk212(qk+rk+pk)gk.The modified scheme (24) is free from the terms 1/(k±1), hence easily solved for k=1(1)n. The difference equation (24) along with the boundary conditions (4) gives a 2n×2n linear system of equations for the unknowns Uk, Vk, k=1(1)n. 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)u(6)(r)=f(r,u(r),u(1)(r),u(2)(r),u(3)(r),u(4)(r),u(6)(r)=u(5)(r)),α<r<β
subject to the necessary boundary conditions:
(27)u(α)=α0,u(2)(α)=α1,u(4)(α)=α2,u(β)=β0,u(2)(β)=β1,u(4)(β)=β2or equivalently,(28)u(2)(r)=v(r),v(2)(r)=w(r),w(2)(r)=f(r,u,u(1)(r),v(r),v(1)(r),w(r),w(1)(r)),hhhhhhhhhhhhhhhhhhhhhhhhhα<r<β,
subject to the natural boundary conditions:
(29)u(α)=α0,v(α)=α1,w(α)=α2,u(β)=β0,v(β)=β1,w(β)=β2.We outline the similar algorithm for (28) as follows:
(30)U~k′=Uk+1-(1-σk2)Uk-σk2Uk-1hkσk(1+σk),V~k′=Vk+1-(1-σk2)Vk-σk2Vk-1hkσk(1+σk),W~k′=Wk+1-(1-σk2)Wk-σk2Wk-1hkσk(1+σk),U~k+1′=(1+2σk)Uk+1-(1+σk)2Uk+σk2Uk-1hkσk(1+σk),V~k+1′=(1+2σk)Vk+1-(1+σk)2Vk+σk2Vk-1hkσk(1+σk),W~k+1′=(1+2σk)Wk+1-(1+σk)2Wk+σk2Wk-1hkσk(1+σk),U~k-1′=-Uk+1+(1+σk)2Uk-σk(2+σk)Uk-1hkσk(1+σk),V~k-1′=-Vk+1+(1+σk)2Vk-σk(2+σk)Vk-1hkσk(1+σk),W~k-1′=-Wk+1+(1+σk)2Wk-σk(2+σk)Wk-1hkσk(1+σk),f~k±1=f(rk±1,Uk±1,U~k±1′,Vk±1,V~k±1′,Wk±1,W~k±1′),U^k′=U~k′+γhk(Vk+1-Vk-1),V^k′=V~k′+γhk(Wk+1-Wk-1),W^k′=W~k′+γhk(f~k+1-f~k-1),
where the values of γ are same as obtained in Section 2.Then, theO(hk3)-approximations for (26) or (28) can be obtained by the following relations for k=1(1)n:
(31)Uk+1-(1+σk)Uk+σkUk-1+h212(pkVk+1+qkVk+rkVk-1)=O(hk5),Vk+1-(1+σk)Vk+σkVk-1+h212(pkWk+1+qkWk+rkWk-1)=O(hk5),Wk+1-(1+σk)Wk+σkWk-1+hk212(pkf~k+qkf^k+rkf~k-1)=O(hk5).The boundary conditions (27) are used to obtain values at k±1 for k=1 and n, respectively. The numerical scheme may be implemented by neglecting O(hk5) 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)u(4)(r)=a(r)u(r)+g(r),α<r<β
or equivalently
(33)u(2)(r)=v(r),v(2)(r)=a(r)u(r)+g(r),α<r<βalong with the boundary conditions (2).Ifa(r),g(r)∈C2[α,β] and a(r)>0 on [α,β], then the boundary value problem (33) has a unique solution. Under these assumptions (see, [27]), there are constants a_>0 and g_ such that a_=mink{ak,ak(1),ak(2)},  g_=mink{gk,gk(1),gk(2)}.For the convergence, the coefficientspk, qk, and rk associated with (5) and (15) must be negative (see, [22]), from which we obtain the condition |σk-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 a(r) and g(r), we obtain the following system of difference equations:
(34)PkZk-1+QkZk+RkZk+1=Sk+Tk(hk),hhhhhhhhhhhhhhhhhhhhhhhk=1(1)n,
where
(35)Pk=[-σkhk212σk(1+σk-σk2)hk212σk(1+σk-σk2)×(ak+hkak(1))-σk],Qk=(σk+1)×[1hk212(σk2+3σk+1)hk224(ak(2)2(σk2+3σk+1)ak+4hkσk(σk-1)ak(1)+hk2σk(σk2-σk+1σk2)ak(2))1],Rk=[-1hk212(σk2+σk-1)hk212(ak(1)(σk2+σk-1)ak+hkσk(σk2-σk+1)ak(1))-1],Sk=[0,hk212((σk3-σk2-σk)gk(2)-(σk2+σk-1)gk(1)-(σk3+4σk2+4σk+1)gkgk(2))hk212]T,Zk=[Uk,Vk]T,Tk(hk)=O(hk5).
Incorporating the boundary values U0=α0, V0=α1, Un+1=β0, and Vn+1=β1, the system of difference equations (34) in the matrix-vector form can be written as
(36)MZ=J+Tk(hk),
where M=[PkQkRk] is the block tridiagonal matrix,  J=[S1-P1α,S2,…,Sn-1,Sn-Rnβ]T, α=[α0,α1]T, and β=[β0,β1]T.Letzk=[uk,vk]T, k=1(1)n, and  z=[z1,…,zn]≅Z, which satisfies
(37)Mz=J.Letε=[ε1,ε2,…,εn]T be the discretization error vector, and let εk=Zk-zk, k=1(1)n be the discretization errors at the node rk. Subtracting (37) from (36), we obtain the error equation:
(38)Mε=Tk(hk).
Also, we obtain
(39)∥Pk∥∞=max2≤k≤n{σk+hk212σk|1+σk-σk2|,σk+hk212σk|1+σk-σk2|·|ak+hkak(1)|,∥Rk∥∞=max1≤k≤n-1{1+hk212|σk2+σk-1|,1+hk212|(σk2+σk-1)ak+hkσk(σk2-σk+1)ak(1)|.Thus for sufficiently smallhk or equivalently as hk→0, we obtain the relations ∥Pk∥∞=σk, k=2(1)n  and  ∥Rk∥∞=1, k=1(1)n-1. Hence, the graph G(M) of the matrix M is strongly connected, and thus the matrix M is irreducible (see, [28]).Further, letΣl be the sum of the lth row sum of the matrix M; then we have the following.Forl=1,
(40)Σl=σl+hl212σl(σl2+5σl+5),Σl+1=σl+hl212alσl(σl2+5σl+5)+hl312al(1)σl(3σl2-σl-1)+hl424al(2)σl(σl3+1).Forl=3(2)2n-3,
(41)Σl=hl22σl(σl+1),Σl+1=hl22alσl(σl+1)+hl36al(1)σl(σl2-1)+hl424al(2)σl(σl3+1).Forl=2n-1,
(42)Σl=1+hl212(5σl2+5σl+1),Σl+1=1+hl212al(5σl2+5σl+1)+hl312al(1)σl(σl2+σl-3)+hl424al(2)σl(σl3+1).This implies the following.Forl=1,
(43)Σl=σl+hl212σl(σl2+5σl+5),Σl+1≥σl+hl212a_σl(σl2+5σl+5)+hl312a_σl(3σl2-σl-1)+hl424a_σl(σl3+1).Forl=3(2)2n-3,
(44)Σl=hl22σl(σl+1),Σl+1≥hl22a_σl(σl+1)+hl36a_σl(σl2-1)+hl424a_σl(σl3+1).Forl=2n-1,
(45)Σl=1+hl212(5σl2+5σl+1),Σl+1≥1+hl212a_(5σl2+5σl+1)+hl312a_σl(σl2+σl-3)+hl424a_σl(σl3+1),For sufficiently small value ofhk, that is, in the limiting case as hk→0, we obtain
(46)Σl≥σl>0,l=1,2,Σl≥0,l=3(1)2n-2,Σl≥1>0,l=2n-1,2n.Hence we find thatM is monotone (see, [29, 30]). Consequently M-1 exists, and M-1≥0.LetMi,l-1 be the (i, l)th element of M-1, and we define
(47)∥M-1∥=max1≤i≤2n∑l=12n|Mi,l-1|,∥T∥=max1≤l≤2n∑l=12n|Tl(hl)|=O(hl3).From the theory of matrix, we know that(48)∑l=12nMi,l-1Σl=1,1≤i≤2n.Thus the following bounds can be estimated with the help of series expansions.Forl=1,
(49)Mi,l-1≤Σl-1=1σl-hl212σl(σl2+5σl+5)+O(hl4),Mi,l+1-1≤Σl+1-1≤1σl-hl2a_12σl(σl2+5σl+5)-hl3a_12σl(3σl2-σl-1)+O(hl4).Forl=3(2)2n-3,
(50)Mi,l-1≤minlΣl-1=2hl2σl(σl+1),Mi,l+1-1≤minlΣl+1-1≤2hl2a_σl(σl+1)-2(σl-1)3hla_σl(σl+1)+σl2-5σl+118a_σl(σl+1)+O(hl).Forl=2n-1,
(51)Mi,l-1≤Σl-1=1-hl212(5σl2+5σl+1)+O(hl4),Mi,l+1-1≤Σl+1-1≤1-hl2a_12(5σl2+5σl+1)-hl3a_12σl(σl2+σl-3)+O(hl4).With the help of  (48), we obtain the following bounds:
(52)∥M-1∥≤2(a_+1)hl2a_σl(σl+1)-2(σl-1)3hla_σl(σl+1)+(σl2-5σl+1+36a_(σl+1)2)18a_σl(σl+1)+O(hl).From (38) and (52), we obtain the following error estimates:
(53)∥ε∥≤∥M-1∥·∥Tl(hl)∥≤O(hl3).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 hk and |σk-5/2|<1/2, σk≠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−15. For the simplicity in computation, we choose σk=σ=constant, for k=1(1)n and define the geometric mesh as follows ([24]):
(54)x0=α,h1={(β-α)(1-σ)(1-σn+1),σ<1,(β-α)(σ-1)(σn+1-1),σ>1.The subsequent mesh spacing is determined byhk+1=σhk, k=1(1)n. 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 inC under Linux operating system with 2 GB operational memory.Example 1.
Consider the fourth-order linear problem (see, [31]) in the polar form:
(55)Δ4u(r)≡(d2dr2+λrddr)2u(r)=(1+2λr+λ(λ-2)r2-λ(λ-2)r3)er,hhhhhhhhhhhhhhhhhhhhh0<r<1.
The analytical solution isu(r)=er. The errors estimates for various values of  n are reported in Tables 1 and 2 for uniform mesh (σ=1) and geometric mesh (σ≠1), respectively.Table 1
The maximum absolute errors for Example1.
n
λ
=
1
λ
=
2
σ
1
0.9937
1
0.996759
8
0.866
e
-
04
0.772
e
-
04
0.290
e
-
05
0.271
e
-
05
16
0.223
e
-
04
0.176
e
-
04
0.190
e
-
05
0.170
e
-
05
32
0.123
e
-
04
0.730
e
-
05
0.640
e
-
08
0.246
e
-
08
64
0.615
e
-
06
0.348
e
-
08
0.365
e
-
09
0.882
e
-
10Table 2
The root mean square errors for Example1.
n
λ
=
1
λ
=
2
σ
1
0.9937
1
0.996759
8
0.664
e
-
04
0.590
e
-
04
0.219
e
-
05
0.205
e
-
05
16
0.166
e
-
04
0.131
e
-
04
0.139
e
-
05
0.125
e
-
05
32
0.877
e
-
05
0.517
e
-
05
0.419
e
-
08
0.239
e
-
08
64
0.417
e
-
06
0.214
e
-
08
0.214
e
-
09
0.483
e
-
10Example 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
(56)u(4)(r)=λu(r)u(2)(r)-λ(r2-1)×(1+4r+r2)e2r-(11+8r+r2)er,hhhhhhhhhhhhhhhhhhhhhh0<r<1.
The analytical solution isu(r)=(1-r2)er. The errors estimates are reported in Table 3 for various values of  n and λ=103.Table 3
The maximum absolute errors and root mean square errors for Example2.
n
σ
1
0.998
1
0.998
8
0.302
e
-
05
0.255
e
-
05
0.235
e
-
05
0.199
e
-
05
16
0.238
e
-
06
0.169
e
-
06
0.179
e
-
06
0.128
e
-
06
32
0.167
e
-
07
0.751
e
-
08
0.124
e
-
07
0.557
e
-
08
64
0.111
e
-
08
0.919
e
-
10
0.819
e
-
09
0.649
e
-
10Example 3.
Consider the sixth-order linear singular problem:(57)Δ6u(r)≡(d2dr2+λrddr)3u(r)=(1+3λ(λ-2)r2+3λ(λ-2)(λ-4)r4)sinh(r)=+λr(3+(λ-2)(λ-7)r2+3(λ-2)(λ-4)r4)cosh(r),0<r<1.
The analytical solution isu(r)=sinh(r). The errors estimates for various values of  n are reported in Tables 4 and 5 for uniform mesh (σ=1) and geometric mesh (σ≠1), respectively.Table 4
The maximum absolute errors for Example3.
n
λ
=
1
λ
=
2
σ
1
1.11
1
0.994
8
0.338
e
-
03
0.189
e
-
03
0.322
e
-
06
0.268
e
-
06
16
0.102
e
-
03
0.292
e
-
04
0.261
e
-
07
0.183
e
-
07
32
0.337
e
-
04
0.481
e
-
05
0.186
e
-
08
0.847
e
-
09
64
0.111
e
-
04
0.662
e
-
06
0.125
e
-
09
0.126
e
-
10Table 5
The root mean square errors for Example3.
n
λ
=
1
λ
=
2
σ
1
1.11
1
0.994
8
0.255
e
-
03
0.140
e
-
03
0.245
e
-
06
0.202
e
-
06
16
0.742
e
-
04
0.197
e
-
04
0.193
e
-
07
0.134
e
-
07
32
0.240
e
-
04
0.270
e
-
05
0.136
e
-
08
0.589
e
-
09
64
0.782
e
-
05
0.243
e
-
06
0.900
e
-
10
0.834
e
-
11Example 4.
Consider the nonlinear problem(58)u(6)(r)=λ(u(r)u(2)(r)+u(2)(r)u(4)(r)+u(r)u(4)(r))+sinh(r)(1-3λsinh(r)),hhhhhhhhhhhhhhhhhhhhhhhhh0<r<1.
The analytical solution isu(r)=sinh(r). The errors estimates are reported in Table 6 for various values of  n and λ=103.Table 6
The maximum absolute errors and root mean square errors for Example4.
n
σ
1
0.9991
1
0.9991
8
0.380
e
-
07
0.302
e
-
07
0.291
e
-
07
0.230
e
-
07
16
0.301
e
-
08
0.186
e
-
08
0.222
e
-
08
0.134
e
-
08
32
0.210
e
-
09
0.587
e
-
10
0.153
e
-
09
0.373
e
-
10
64
0.127
e
-
10
0.913
e
-
11
0.912
e
-
11
0.654
e
-
11


## 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.



---

*Source: 614508-2013-10-24.xml*