# A Class of Numerical Methods for the Solution of Fourth-Order Ordinary Differential Equations in Polar Coordinates

Authors: Jyoti Talwar; R. K. Mohanty
Journal: Advances in Numerical Analysis
 (2012)

---

## Abstract

In this piece of work using only three grid points, we propose two sets of numerical methods in a coupled manner for the solution of fourth-order ordinary differential equationuiv(x)=f(x,u(x),u′(x),u′′(x),u′′′(x)), a<x<b, subject to boundary conditions u(a)=A0,  u′(a)=A1, u(b)=B0, and u′(b)=B1, where A0, A1, B0, and B1 are real constants. We do not require to discretize the boundary conditions. The derivative of the solution is obtained as a byproduct of the discretization procedure. We use block iterative method and tridiagonal solver to obtain the solution in both cases. Convergence analysis is discussed and numerical results are provided to show the accuracy and usefulness of the proposed methods.


---

## Body

## 1. Introduction

Consider the fourth-order boundary value problem(1.1)uiv(x)=f(x,u(x),u′(x),u′′(x),u′′′(x)),a<x<b,
subject to the prescribed natural boundary conditions(1.2)u(0)=A0,u′(0)=A1,u(1)=B0,u′(1)=B1,
or equivalently, for u′(x)=v(x),(1.3)uiv(x)=f(x,u(x),v(x),u′′(x),v′′(x)),a<x<b,
subject to the natural boundary conditions
(1.4)u(0)=A0,v(0)=A1,u(1)=B0,v(1)=B1,
where A0, A1, B0, and B1 are real constants and -∞<a≤x≤b<∞.Fourth-order differential equations occur in a number of areas of applied mathematics, such as in beam theory, viscoelastic and inelastic flows, and electric circuits. Some of them describe certain phenomena related to the theory of elastic stability. A classical fourth-order equation arising in the beam-column theory is the following (see Timoshenko [1]):(1.5)EId4udx4+Pd2udx2=q,
where u is the lateral deflection, q is the intensity of a distributed lateral load, P is the axial compressive force applied to the beam, and EI represents the flexural rigidity in the plane of bending. Various generalizations of the equation describing the deformation of an elastic beam with different types of two-point boundary conditions have been extensively studied via a broad range of methods.The existence and uniqueness of solutions of boundary value problems are discussed in the papers and book of Agarwal and Krishnamoorthy, Agarwal and Akrivis (see [2–5]). Several authors have investigated solving fourth-order boundary value problem by some numerical techniques, which include the cubic spline method, Ritz method, finite difference method, multiderivative methods, and finite element methods (see [6–16]). In the 1980s, Usmani et al. (see [17–19]) worked on finite difference methods for solving [p(x)y′′]′′+q(x)y=r(x) and finite difference methods for computing eigenvalues of fourth-order linear boundary value problem. In 1984, Twizell and Tirmizi (see [20]) developed multi-derivative methods for linear fourth-order boundary value problems. In 1984, Agarwal and Chow (see [21]) developed iterative methods for a fourth-order boundary value problem. In 1991, O’Regan (see [13]) worked on the solvability of some fourth-(and higher) order singular boundary value problems. In 1994, Cabada (see [22]) developed the method of lower and upper solutions for fourth- and higher-order boundary value problems. In 2005 Franco et al. (see [23]) dealt with some fourth-order problems with nonlinear boundary conditions. In 2006, Noor and Mohyud-Din (see [12]) used the variational iteration method to solve fourth-order boundary value problems and further developed the homotopy perturbation method for solving fourth order boundary value problems. Some of these methods use transformation in order to reduce the equation into more simple equation or system of equations and some other methods give the solution in a series form which converges to the exact solution. Later, Han and Li (see [24]) worked on some fourth-order boundary value problems.In this paper we present the finite difference methods for the solution of fourth-order boundary value problem. We discretize the interval[a,b] into N+1 subintervals each of width h=(b-a)/(N+1), where N is a positive integer. We seek the solution of (la) or (2a) at the grid points, xk=kh, k=1,2,…,N. Let uk and uk′ denote the approximate solutions, and let Uk=u(xk) and Uk′=u′(xk) be the exact solution values of u(x) and u′(x) at the grid point x=xk, respectively. Also, x0=a and xN+1=b.Using the second-order central differences, we obtain a five-point difference formula for (1.1), which requires the use of fictitious points outside [a,b]. The accuracy of the numerical solution depends upon the boundary approximation used. The finite difference method discussed here is based only on three grid points for second-and fourth-order methods. Therefore, no fictitious points are required for incorporating the boundary conditions. Here, we use a combination of the value of the solution u(x) and its derivative u′(x) to derive the difference scheme using three grid points.Since we need to solve a coupled system of equations at each mesh point, the block successive overrelaxation (BSOR) iterative method is used.


## 2. The Finite Difference Method

The method is described as follows: fork=1(1)N, let(2.1)u̅k′′=(uk+1′-uk-1′)2h,u̅k′′′=(uk+1′-2uk′+uk-1′)h2,f̅k=f(xk,uk,uk′,u̅k′′,u̅k′′′).
Then, the difference method of order two for the given differential (1.1) is given by,(2.2a)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h46f̅k,
and the corresponding difference method for the derivative u′(x) is given by(2.2b)-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=0.Also, let (2.3)u̅k±1′′=(±3uk±1′∓4uk′±uk∓1′)2h,(2.4)u̅̅k′′′=15(uk+1-uk-1)2h3-3(uk+1′+8uk′+uk-1′)2h2,(2.5)u̅̅k±1′′=-(11uk±1-16uk+5uk∓1)2h2±(4uk±1′-uk∓1′)h,(2.6)u̅̅k′′=2(uk+1-2uk+uk-1)h2-(uk+1′-uk-1′)2h,(2.7)u̅k±1′′′=∓(27uk±1-48uk+21uk∓1)2h3+(15uk±1′-9uk∓1′)2h2,(2.8)u̅̅k±1′′′=∓(99uk±1-48uk-51uk∓1)2h3+(39uk±1′+96uk′+15uk∓1′)2h2,
and set(2.9a)f̅k±1=f(xk±1,uk±1,uk±1′,u̅k±1′′,u̅k±1′′′),(2.9b)f̅̅k±1=f(xk±1,uk±1,uk±1′,u̅̅k±1′′,u̅̅k±1′′′),(2.9c)f̅̅k=f(xk,uk,uk′,u̅̅k′′,u̅̅k′′′).Then, the difference method of order four for the differential equation and the corresponding difference method for the derivative uk′ are given by(2.10)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h490(f̅̅k+1+f̅̅k-1+13f̅̅k),(2.11)-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=h460(f̅k+1-f̅k-1).Note thatu0, u0′, uN+1 and uN+1′, are prescribed. It is convenient to express the above finite difference schemes in block tridiagonal matrix form. If the differential equation is linear, the resulting block tridiagonal linear system can be solved using the block Gauss-Seidel (BGS) iterative method. If the differential equation is nonlinear, the system can be solved using the Newton nonlinear block successive overrelaxation (NBSOR) method (see [25, 26]).


## 3. Derivation of the Difference Scheme

For the derivation of the method we follow the approaches given by Chawla [27] and Mohanty [10]. In this section we discuss the derivation of the difference methods and the block iterative methods.At the grid pointxk, the given differential equation can be written as(3.1a)ukiv=f(xk,uk,uk′,uk′′,uk′′′)=fk,k=1(1)N.
Similarly,(3.1b)fk±1=f(xk±1,uk±1,uk±1′,uk±1′′,uk±1′′′),k=1(1)N.Using the Taylor expansion about the grid pointxk, we first obtain(3.2)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h490(fk+1+fk-1+13fk)+T1,
where T1=O(h8).Now, we need theO(h2) approximation for uk±1′′′. Let(3.3)u̅k+1′′′=1h3[a10uk+a11uk+1+a12uk-1]+1h2[b10uk′+b11uk+1′+b12uk-1′].
Expanding each term on the right-hand side of (3.3) in the Taylor series about the point xk and equating the coefficients of hp (p=-3,-2,-1,0, and 1)to zero, we get(3.4)(a10,a11,a12,b10,b11,b12)=(24,-272,-212,0,152,-92).Thus, we obtain(3.5a)u̅k+1′′′=-12h3[27uk+1-48uk+21uk-1]+12h2[15uk+1′-9uk-1′]=uk+1′′′-2h25ukv+O(h3).
Similarly, we obtain(3.5b)u̅k-1′′′=12h3[27uk-1-48uk+21uk+1]+12h2[15uk-1′-9uk+1′]=uk-1′′′-2h25ukv-O(h3).
Further, from (2.3) we obtain:(3.6)u̅k±1′′=uk±1′′-h23ukiv±O(h3).
Hence, we can verify that f̅k±1=fk±1+O(h2).Next, we obtain theO(h4) approximation for uk′′. Let(3.7)u̅̅k′′=1h2[a20uk+a21uk+1+a22uk-1]+1h[b20uk′+b21uk+1′+b22uk-1′].
Using the Taylor series expansion and equating the coefficients of hp (p=-2,-1,0,1,2, and 3) to zero, we get(3.8)(a20,a21,a22,b20,b21,b22)=(-4,2,2,0,-12,12).Therefore,(3.9)u̅̅k′′=2h2[uk+1-2uk+uk-1]-12h[uk+1′-uk-1′]=uk′′+O(h4).
Similarly,
(3.10)u̅̅k′′′=152h3[uk+1-uk-1]-32h2[uk+1′+8uk′+uk-1′]=uk′′′+O(h4),(3.11)u̅̅k±1′′=-12h2[11uk±1-16uk+5uk∓1]±1h[4uk±1′-uk∓1′]=uk±1′′∓h315ukv+O(h4).
Next, we obtain the O(h3) approximation for uk+1′′′. Let(3.12)u̅̅k+1′′′=1h3[a30uk+a31uk+1+a32uk-1]+1h2[b30uk′+b31uk+1′+b32uk-1′].
Equating the coefficients of hp (p=-3,-2,-1,0,1, and 2) to zero, we get(3.13)(a30,a31,a32,b30,b31,b32)=(24,-992,512,48,392,152).Thus, we obtain(3.14)u̅̅k+1′′′=-12h3[99uk+1-48uk-51uk-1]+12h2[39uk+1′+96uk′+15uk-1′]=uk+1′′′-h310ukvi+O(h4).
Similarly,(3.15)u̅̅k-1′′′=12h3[99uk-1-48uk-51uk+1]+12h2[39uk-1′+96uk′+15uk+1′]=uk-1′′′+h310ukvi+O(h4).Let(3.16)αk=∂f∂uk′′,βk=∂f∂uk′′′.From (3.5a), (3.5b), (3.6), and (2.9a) it follows that f̅k±1 provides the O(h2) approximation for fk±1 and(3.17)f̅k±1=fk±1-h23ukivαk-2h25ukvβk±O(h3).Also, from (2.9b), (3.11), (3.14), and (3.15) we obtain the O(h3) approximation for fk±1:(3.18)f̅̅k±1=fk±1∓h315ukvαk∓h310ukviβk+O(h4).
From (2.9c), (3.9), and (3.10) it follows that f̿k provides an O(h4) approximation for fk(3.19)f̅̅k=fk+O(h4).
From (3.18) and (3.19) we get(3.20)f̅̅k+1+f̅̅k-1+13f̅̅k=fk+1+fk-1+13fk+O(h4).
With the help of (3.2) and (3.20), from (2.10), we obtain that the local truncation error associated with the difference scheme (2.10) is of O(h8). Similarly, we can verify that the local truncation error associated with the difference scheme (2.11) is of O(h7).Thus, we have the following result.Theorem 3.1.
Consider the fourth-order nonlinear ordinary differential equation (1.1) along with the boundary conditions (1.2). Let u∈C8[a,b], and the function f is sufficiently differentiable with respect to its arguments. Then, the difference methods (2.10) and (2.11) with the approximations of u′′ and u′′′ listed in (2.3)–(2.8) are of O(h4).If the differential (1.1) is linear, then the difference method (2.10) and (2.11) in the matrix form can be written as(3.21)[A11A12A21A22][uu′]=[d1d2],
where A11, A12, A21, and A22 are the Nth-order tri-diagonal matrices and d1 and d2 are vectors consisting of right-hand side functions and some boundary conditions associated with the block system.The block successive over relaxation (BSOR) method (See Mohanty and Evans [26, 28]) is given by
(3.22)A11u(n+1)=ω[-A12(u′)(n)+d1]+(1-ω)A11u(n),n=0,1,2,…,A22(u′)(n+1)=ω[-A21u(n+1)+d2]+(1-ω)A22(u′)(n),n=0,1,2,…,
where ω is a parameter known as relaxation parameter. With ω=1, the BSOR method reduces to block Gauss-Seidel (BGS) method. If ω>1 or ω<1, we have overrelaxation or underrelaxation, respectively.Iff(x,u(x),v(x),u′(x),v′(x))is nonlinear, the difference equations (2.2a), (2.2b) or (2.10), (2.11) form a coupled nonlinear system. To solve the coupled nonlinear system we apply the Newton NBSOR method.We first write the difference equations (2.2a), (2.2b) or (2.10), (2.11) as(3.23)Φ(u,v)=0,Ψ(u,v)=0,
where(3.24)u=[u1u2⋮uN],v=[v1v2⋮vN],Φ(u,v)=[ϕ1(u,v)ϕ2(u,v)⋮ϕN(u,v)],Ψ(u,v)=[ψ1(u,v)ψ2(u,v)⋮ψN(u,v)].Let(3.25)J=[T11T12T21T22]be the Jacobian ofΦ and Ψ, which is the 2Nth-order block tridiagonal matrix, where
(3.26)T11=∂(ϕ1,ϕ2,…,ϕN)∂(u1,u2,…,uN)=[∂ϕ1∂u1∂ϕ1∂u2∂ϕ2∂u1∂ϕ2∂u2∂ϕ2∂u300∂ϕN∂uN-1∂ϕN∂uN],T12=∂(ϕ1,ϕ2,…,ϕN)∂(v1,v2,…,vN),T21=∂(ψ1,ψ2,…,ψN)∂(u1,u2,…,uN),T22=∂(ψ1,ψ2,…,ψN)∂(v1,v2,…,vN)are theNth-order tridiagonal matrices.In the NBSOR method starting with any initial approximation(u(0),v(0)) of (u(k),v(k)), k=1(1)N, we define(3.27)u(n+1)=u(n)+Δu(n),n=0,1,2,…,v(n+1)=v(n)+Δv(n),n=0,1,2,…,
where Δu and Δv are intermediate values obtained by solving the matrix equation for NBSOR method given by(3.28)[T11T12T21T22][ΔuΔv]=[-Φ-Ψ].
The above system can be solved for Δu(n) and Δv(n) by using the block SOR method (inner iterative method) as follows:(3.29)T11Δu(n+1)=ω[-Φ(u(n),v(n))-T12Δv(n)]+(1-ω)T11Δu(n),n=0,1,2,…,T22Δv(n+1)=ω[-Ψ(u(n),v(n))-T21Δu(n+1)]+(1-ω)T22Δv(n),n=0,1,2,…,
where ω∈(0,2) is a relaxation parameter and n=0,1,2,…. The above system of equations can be solved by using the tridiagonal solver. In order for this method to converge it is sufficient that the initial approximation (u(0),v(0)) be close to the solution.


## 4. Convergence and Stability Analysis

Consider the model problemuiv=f(x), where f is a function of x only. Applying the fourth-order difference methods (2.10) and (2.11) to the above equation, we get(4.1)(uk-1-2uk+uk+1)+h2(vk-1-vk+1)=-h4180(fk+1+fk-1+13fk),3h(uk-1-uk+1)+(vk-1+4vk+vk+1)=h360(fk+1-fk-1),
where we denote u′=v.Let us denote byP=[1,0,1], L=[1,1,1], and M=[1,0,-1] the Nth-order tridiagonal matrices.The system of (4.1) can be written in the block form as
(4.2)[(L-3I)h2M3hM(L+3I)][uv]=[d1d2],
where u and v are two N-dimensional solution vectors and d1, d2 are vectors consisting of right-hand side functions and some boundary values associated with (4.1). The BSOR method for the scheme is(4.3)u(n+1)=(1-ω)u(n)-ωh2(L-3I)-1M(v)(n)+ω(L-3I)-1d1,(v)(n+1)=(1-ω)(v)(n)-3ωh(L+3I)-1Mu(n+1)+ω(L+3I)-1d2,
The associated block SOR and block Jacobi iteration matrices are given by(4.4)Lω=[(1-ω)I-ωh2(L-3I)-1M-3ωh(L+3I)-1M(1-ω)I],B=[0-h2(L-3I)-1M-3h(L+3I)-1M0],
and λ and η are the eigenvalues associated with the corresponding matrices Lω and B, which are related by the equation(4.5)(λ+ω-1)2=λω2η2.
Let [v1v2] be the eigenvector associated with the eigenvalue η so that(4.6)[0-h2(L-3I)-1M-3h(L+3I)-1M0][v1v2]=[v1v2]η.
That is,(4.7)-h2(L-3I)-1Mv2=ηv1,-3h(L+3I)-1Mv1=ηv2.
On eliminating v2, we obtain(4.8)32(L-3I)-1M(L+3I)-1Mv1=η2v1.
The rate of convergence of the BSOR method is given by -log(ρ(Lω)).The rate of convergence of the BSOR method is dependent on the eigenvalues ofB (through relation (4.5)), which are given by (3τ/2)=η2, where τ are the eigenvalues of (L-3I)-1M(L+3I)-1M.Hence, we can determine the optimal parameter as(4.9)ω0=21+1-(3/2)τ̅,
where τ̅=S((L-3I)-1M(L+3I)-1M).Thus, we can determine the convergence factor(4.10)λ̅=ω0-1=1-1-(3/2)τ̅1+1-(3/2)τ̅.
For convergence, we must have |λ̅|<1 to give the range
(4.11)0<τ̅<23Hence, we get the convergence.Thus, we have the following result.Theorem 4.1.
The iteration method of the form (4.3) for the solution of uiv=f(x) converges if 0<τ̅<2/3, where τ̅=S((L-3I)-1M(L+3I)-1M), L=[1,1,1],andM=[1,0,-1] are the N×N tridiagonal matrices.Now, we discuss the stability analysis.An iterative method for (4.1) can be written as (4.12)u(k+1)=12Pu(k)+h4Mv(k)+RHU,v(k+1)=-34hMu(k)-14Pv(k)+RHV,
where u(k), v(k) are solution vectors at the kth iteration and RHU, RHV are right-hand side vectors consisting of boundary and homogenous function values.The above iterative method in matrix form can be written as(4.13)[u(k+1)v(k+1)]=G[u(k)v(k)]+RH,
where(4.14)G=[12Ph4M-34hM-14P],RH=[RHURHV].The eigenvalues ofP and M are 2coskπ/(N+1) and 2icoskπ/(N+1), respectively, where k=1,2,…,N. The characteristic equation of the matrix G is given by(4.15)det⁡[12P-ξINh4M-34hM-14P-ξIN]=0.Thus the eigenvalues ofG are given by
(4.16)det⁡[-14P-ξIN]×det⁡[(12P-ξIN)+316M(-14P-ξIN)-1M]=0.
The proposed iterative method (4.13) is stable as long as the maximum absolute eigenvalues of the iteration matrix are less than or equal to one. It has been verified computationally that all eigenvalues of the system (4.16) are less than one. Hence, the iterative method (4.1) is stable.


## 5. Application to Singular Equation

Consider a singular fourth-order linear ordinary differential equation of the form(5.1)Δ4u≡(d2dr2+1rddr)2u=f(r),0<r<1,
or equivalently(5.2)uiv=b(r)u′′′+c(r)u′′+d(r)u′+f(r),0<r<1,
where
(5.3)b(r)=-2/r,c(r)=1/r2,d(r)=-1/r3.The above equation represents fourth-order ordinary differential equation in cylindrical polar coordinates.The boundary conditions are given by(5.4)u(0)=A0,u′(0)=A1,u(1)=B0,u′(1)=B1,
where A0, A1, B0, and B1 are constants.Applying the difference scheme (2.2a), (2.2b) to the singular equation (5.2), we obtain a second-order difference method(5.5)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h46(bku̅k′′′+cku̅k′′+dkuk′+fk),-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=0,k=1(1)N.
Applying the fourth-order difference scheme (2.10) to the singular equation (5.2), we obtain(5.6)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h490[(bk+1u̅̅k+1′′′+ck+1u̅̅k+1′′+dk+1uk+1′+fk+1)+(bk-1u̅̅k-1′′′+ck-1u̅̅k-1′′+dk-1uk-1′+fk-1)+(13bku̅̅k′′′+13cku̅̅k′′+13dkuk′+13fk)],k=1(1)N,
where
(5.7)bk±1=b(rk±1),ck±1=c(rk±1),dk±1=d(rk±1),fk±1=f(rk±1).Note that the scheme fails when the solution is to be determined atk=1. We overcome the difficulty by modifying the method in such a way that the solutions retain the order and accuracy even in the vicinity of the singularity r=0 (see [29]). We consider the following approximations:(5.8)bk±1=bk±hbk′+h22!bk′′+O(±h3),ck±1=ck±hck′+h22!ck′′+O(±h3),dk±1=dk±hdk′+h22!dk′′+O(±h3),fk±1=fk±hfk′+h22!fk′′+O(±h3).Using the approximation (5.8) in (5.6) and neglecting higher-order terms, we can rewrite (5.6) in compact operator form as(5.9)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h290(-24bk′+18ck-4h2ck′′)δx2uk+h180(45bk-75h2bk′′-6h2ck′)(2μxδx)uk+h490(-45bkh2+75bk′′+6ck′+15dk+h2dk′′)uk′+h2180(15bk+27h2bk′′+6h2ck′+2h2dk)(2μxδx)uk′+h3180(24bk′-3ck+5h2ck′′+2h2dk′)(2μxδx)uk+h490(15fk+h2fk′′),k=1(1)N.Similarly, using the difference scheme (2.11), a fourth-order approximation for the derivative u′ for the singular equation (5.2) in the compact form may be written as(5.10)-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=h60(-24bk)δx2uk+h260(-3bk′)(2μxδx)uk+h360(2ck+3bk′)δx2uk′+h260(12bk+h2(ck′+dk))(2μxδx)uk′+h360(6bk′+2h2dk′)uk′+h460(2hfk′),k=1(1)N,
where δruk=uk+1/2-uk-1/2, μruk=(uk+1/2+uk-1/2)/2.Finite difference equations (5.9) and (5.10) along with the boundary conditions (5.4) gives a 2N×2N linear system of equations for the unknowns u1,u2,…,uN,u1′,u2′,…,uN′. The resulting block tri-diagonal system can be solved using the BGS method. The schemes (5.9) and (5.10) are free from the terms 1/(k±1), hence very easily solved for k=1(1)N in the region (0,1).Consider the coupled nonlinear singular equations(5.11)uIV=a(r)[u′v′′+v′u′′]+f(r),0<r<1,vIV=-a(r)u′u′′+g(r),0<r<1,
where a(r)=1/r, with known boundary conditions u(0), v(0), u′(0), v′(0), u(1), v(1), u′(1), and v′(1). The coupled equations represent model equations of equilibrium for a load symmetrical about the centre (see [30]).The second-order difference scheme for solving the system (5.11) is given by(5.12)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=h46[ak(uk′v̅k′′+vk′u̅k′′)+fk],-2(vk+1-2vk+vk-1)+h(vk+1′-vk-1′)=h46[-akuk′u̅k′′+gk],-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=0,-3(vk+1-vk-1)+h(vk+1′+4vk′+vk-1′)=0,
where ak=a(rk), fk=f(rk), and gk=g(rk).The fourth-order difference scheme for solving the system (5.11) for u, v, u′, and v′ is given by
(5.13)-2(uk+1-2uk+uk-1)+h(uk+1′-uk-1′)=H1[15fk+h2fk′′]+H1a11[uk+1′{-12h2(11vk+1-16vk+5vk-1)+1h(4vk+1′-vk-1′)}+vk+1′{-12h2(11uk+1-16uk+5uk-1)+1h(4uk+1′-uk-1′)}]+H1a12[uk-1′{-12h2(11vk-1-16vk+5vk+1)-1h(4vk-1′-vk+1′)}+vk-1′{-12h2(11uk-1-16uk+5uk+1)-1h(4uk-1′-uk+1′)}]+H1a10[uk′{2h2(vk+1-2vk+vk-1)-12h(vk+1′-vk-1′)}+vk′{2h2(uk+1-2uk+uk-1)-12h(uk+1′-uk-1′)}],-2(vk+1-2vk+vk-1)+h(vk+1′-vk-1′)=H1[15gk+h2gk′′]-H1a11uk+1′[-12h2(11uk+1-16uk+5uk-1)+1h(4uk+1′-uk-1′)]-H1a12uk-1′[-12h2(11uk-1-16uk+5uk+1)-1h(4uk-1′-uk+1′)]-13H1a10uk′[2h2(uk+1-2uk+uk-1)-12h(uk+1′-uk-1′)],-3(uk+1-uk-1)+h(uk+1′+4uk′+uk-1′)=H2(2hfk′)+H2a11[uk+1′{3vk+1′-4vk′+vk-1′2h}+v′k+1{3uk+1′-4uk′+uk-1′2h}]-H2a12[uk-1′{-3vk-1′+4vk′-vk+1′2h}+vk-1′{-3uk-1′+4uk′-uk+1′2h}],-3(vk+1-vk-1)+h(vk+1′+4vk′+vk-1′)=H2(2hgk′)-H2a11uk+1′{3uk+1′-4u′k+uk-1′2h}-H2a12uk-1′{-3uk-1′+4u′k-uk+1′2h},
where(5.14)ak=a(rk),fk=f(rk),gk=g(rk),a11=ak+hak′+h22ak′′,a12=ak-hak′+h22ak′′,a10=ak.The scheme (5.13) is free from the terms 1/k±1, hence very easily solved for k=1(1)N in the region (0,1). The system (5.13) can be solved using the NBSOR method.Consider the boundary value problem [31](5.15)yiv-λyy′′=f(x),y(0)=V0,y(1)=V1,y′(0)=0,y′(1)=0.This arises from the time-dependent Navier-Strokes equations for axisymmetric flow of an incompressible fluid contained between infinite disks that occupy the planesz=-d and z=d. The disks are porous and the fluid is injected or extracted normally with velocity V0atz=-d and V1 at z=d. Here, =d/ν, where νis the kinematic viscosity.


## 6. Numerical Illustrations

To illustrate our method and to demonstrate computationally its convergence, we have solved the following linear problem using the BGS method (See [25, 32–35]), whose exact solution is known to us. We have taken [0,1] as our region of integration. The right-hand side functions and the boundary conditions are obtained using the exact solution. The iterations were stopped when the absolute error tolerance became ≤10−12. All computations were performed using double length arithmetic.Problem 1.
The problem is to solve (5.2) subject to the boundary conditions (5.4). The exact solution is u=r4sinr. The root mean square errors (RMSEs) are tabulated in Table 1. The graph of the errors for N=32 is given in Figure 1.Table 1
The RMSE.
2nd-order method4th-order methodhMAERMSEMAERMSE1/8u0.9512 (−03)0.6737 (−03)0.6018 (−03)0.4089 (−03)u′0.3162 (−02)0.2137 (−02)0.2961 (−02)0.1402 (−02)1/16u0.2212 (−03)0.1501 (−03)0.3799 (−04)0.2502 (−04)u′0.7021 (−03)0.4995 (−03)0.3022 (−03)0.9846 (−04)1/32u0.5410 (−04)0.3576 (−04)0.1989 (−05)0.1304 (−05)u′0.1729 (−03)0.1206 (−03)0.2302 (−05)0.5629 (−05)1/64u0.1339 (−04)0.8769 (−05)0.9996 (−07)0.6578 (−07)u′0.4298 (−04)0.2972 (−04)0.1562 (−05)0.2983 (−06)1/128u0.3363 (−05)0.2174 (−05)0.4198 (−08)0.3345 (−08)u′0.1072 (−04)0.7383 (−05)0.9892 (−07)0.1538 (−07)Figure 1Problem 2.
The boundary value problem is to solve (5.15). The exact solution is (1-x2)exp(x). The maximum absolute errors (MAE) and RMSE are tabulated in Table 2.Table 2
The MAE and RMSE.
2nd-order method4th-order methodhMAERMSEMAERMSEλ=1.01/8u0.8435 (−04)0.5652 (−04)0.2744 (−05)0.1880 (−05)u′0.3552 (−03)0.2322 (−03)0.1058 (−04)0.7178 (−05)1/16u0.2292 (−04)0.1499 (−04)0.1854 (−06)0.1217 (−06)u′0.7762 (−04)0.5391 (−04)0.6695 (−06)0.4351 (−06)1/32u0.5844 (−05)0.3771 (−05)0.1151 (−07)0.7446 (−08)u′0.1885 (−04)0.1322 (−04)0.4115 (−07)0.2618 (−07)λ=1001/8u0.1961 (−02)0.1401 (−02)0.5662 (−04)0.3877 (−04)u′0.6514 (−02)0.4687 (−02)0.3824 (−03)0.1910 (−03)1/16u0.4808 (−03)0.3312 (−03)0.3821 (−05)0.2558 (−05)u′0.1597 (−02)0.1110 (−02)0.1929 (−04)0.9922 (−05)1/32u0.1197 (−03)0.8103 (−04)0.2435 (−06)0.1609 (−06)u′0.4023 (−03)0.2714 (−03)0.1135 (−05)0.5843 (−06)Problem 3.
The system of nonlinear equation (5.11) is to be solved subject to the natural boundary conditions. The exact solutions are u=cos(r) and v=exp(r). The MAE and RMSE are tabulated in Table 3.Table 3
The MAE and RMSE.
2nd-order method4th-order methodhMAERMSEMAERMSE1/10u0.7141 (−05)0.4859 (−05)0.2673 (−05)0.1818 (−05)v0.3339 (−05)0.2240 (−05)0.6974 (−06)0.4745 (−06)u′0.2488 (−04)0.1683 (−04)0.1266 (−04)0.6538 (−05)v′0.1126 (−04)0.8137 (−05)0.2548 (−05)0.1656 (−05)1/20u0.1806 (−05)0.1182 (−05)0.2467 (−06)0.1639 (−06)v0.8550 (−06)0.5595 (−06)0.4846 (−07)0.3212 (−07)u′0.6486 (−05)0.4128 (−05)0.1167 (−05)0.5946 (−06)v′0.2659 (−05)0.1957 (−05)0.1795 (−06)0.1119 (−06)1/40u0.4514 (−06)0.2917 (−06)0.2048 (−07)0.1346 (−07)v0.2152 (−06)0.1391 (−06)0.3315 (−08)0.2161 (−08)u′0.1620 (−05)0.1020 (−05)0.1027 (−06)0.4911 (−07)v′0.6757 (−06)0.4828 (−06)0.1222 (−07)0.7515 (−08)1/80u0.1125 (−06)0.7216 (−07)0.1510 (−08)0.9836 (−09)v0.5307 (−07)0.3442 (−07)0.2931 (−09)0.1883 (−09)u′0.4038 (−06)0.2525 (−06)0.8190 (−08)0.3662 (−08)v′0.1686 (−06)0.1192 (−06)0.9925 (−09)0.6537 (−09)1/160u0.3050 (−07)0.1946 (−07)0.1213 (−09)0.7945 (−10)v0.1461 (−07)0.9331 (−08)0.1358 (−10)0.8769 (−11)u′0.1076 (−06)0.6812 (−07)0.6420 (−09)0.2913 (−09)v′0.4591 (−07)0.3237 (−07)0.5009 (−10)0.3048 (−10)Observe that for the linear singular problem forα=1, h1=1/64, h2=1/128, log(eh1/eh2)/log(h1/h2)=log(0.6578(-07)/0.3345(-08))/log2≈4. Thus, we obtain fourth-order convergence for u. Similarly, we get fourth-order convergence for u′.


## 7. Concluding Remarks

The numerical results confirm that the proposed finite difference methods yield second- and fourth-order convergence for the solution and its derivative for the fourth-order ordinary differential equation. Difference formulas for mesh points near the boundary are obtained without using the fictitious points. The proposed method is applicable to problems in polar coordinates and the derivative of the solution is obtained as the by-product of the method. We employ the BGS method to solve the block matrix systems of the linear singular problem and the BSOR method to solve the nonlinear singular problem. We have solved here a physical problem that arises in the axisymmetric flow of an incompressible fluid.



---

*Source: 626419-2012-03-15.xml*