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

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 equation 𝑢^𝑖𝑣(𝑥)=𝑓(𝑥,𝑢(𝑥),𝑢^′(𝑥),𝑢^′′(𝑥),𝑢^′′′(𝑥)), 𝑎<𝑥<𝑏, subject to boundary conditions 𝑢(𝑎)=𝐴₀, 𝑢^′(𝑎)=𝐴₁, 𝑢(𝑏)=𝐵₀, and 𝑢^′(𝑏)=𝐵₁, where 𝐴₀, 𝐴₁, 𝐵₀, and 𝐵₁ 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.

1. Introduction

Consider the fourth-order boundary value problem 𝑢^𝑖𝑣(𝑥)=𝑓𝑥,𝑢(𝑥),𝑢^′(𝑥),𝑢^′′(𝑥),𝑢^′′′(𝑥),𝑎<𝑥<𝑏, (1.1) subject to the prescribed natural boundary conditions𝑢(0)=𝐴₀,𝑢^′(0)=𝐴₁,𝑢(1)=𝐵₀,𝑢^′(1)=𝐵₁, (1.2) or equivalently, for 𝑢^′(𝑥)=𝑣(𝑥),𝑢^𝑖𝑣(𝑥)=𝑓𝑥,𝑢(𝑥),𝑣(𝑥),𝑢^′′(𝑥),𝑣^′′(𝑥),𝑎<𝑥<𝑏, (1.3) subject to the natural boundary conditions 𝑢(0)=𝐴₀,𝑣(0)=𝐴₁,𝑢(1)=𝐵₀,𝑣(1)=𝐵₁, (1.4) where 𝐴₀, 𝐴₁, 𝐵₀, and 𝐵₁ are real constants and −∞<𝑎≤𝑥≤𝑏<∞.

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) where 𝑢 is the lateral deflection, 𝑞 is the intensity of a distributed lateral load, 𝑃 is the axial compressive force applied to the beam, and 𝐸𝐼 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 [𝑝(𝑥)𝑦^′′]^′′+𝑞(𝑥)𝑦=𝑟(𝑥) 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 𝑁+1 subintervals each of width ℎ=(𝑏−𝑎)/(𝑁+1), where 𝑁 is a positive integer. We seek the solution of (la) or (2a) at the grid points, 𝑥_𝑘=𝑘ℎ, 𝑘=1,2,…,𝑁. Let 𝑢_𝑘 and 𝑢^′_𝑘 denote the approximate solutions, and let 𝑈_𝑘=𝑢(𝑥_𝑘) and 𝑈^′_𝑘=𝑢^′(𝑥_𝑘) be the exact solution values of 𝑢(𝑥) and 𝑢^′(𝑥) at the grid point 𝑥=𝑥_𝑘, respectively. Also, 𝑥₀=𝑎 and 𝑥_𝑁+1=𝑏.

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 𝑢(𝑥) and its derivative 𝑢^′(𝑥) 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: for 𝑘=1(1)𝑁, let𝑢_𝑘^′′=𝑢^′_𝑘+1−𝑢^′_𝑘−1,2ℎ𝑢_𝑘^′′′=𝑢^′_𝑘+1−2𝑢^′_𝑘+𝑢^′_𝑘−1ℎ²,𝑓_𝑘=𝑓𝑥_𝑘,𝑢_𝑘,𝑢^′_𝑘,𝑢_𝑘^′′,𝑢_𝑘^′′′. (2.1) Then, the difference method of order two for the given differential (1.1) is given by,−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴6𝑓_𝑘, (2.2a) and the corresponding difference method for the derivative 𝑢^′(𝑥) is given by−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1=0. (2.2b)Also, let 𝑢^′′_𝑘±1=±3𝑢^′_𝑘±1∓4𝑢^′_𝑘±𝑢^′_𝑘∓1,2ℎ𝑢_𝑘^′′′=15𝑢_𝑘+1−𝑢_𝑘−12ℎ³−3𝑢^′_𝑘+1+8𝑢^′_𝑘+𝑢^′_𝑘−12ℎ²,𝑢^′′_𝑘±1=−11𝑢_𝑘±1−16𝑢_𝑘+5𝑢_𝑘∓12ℎ²±4𝑢^′_𝑘±1−𝑢^′_𝑘∓1ℎ,𝑢_𝑘^′′=2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1ℎ²−𝑢^′_𝑘+1−𝑢^′_𝑘−12ℎ,𝑢^′′′_𝑘±1=∓27𝑢_𝑘±1−48𝑢_𝑘+21𝑢_𝑘∓12ℎ³+15𝑢^′_𝑘±1−9𝑢^′_𝑘∓12ℎ²,𝑢^′′′_𝑘±1=∓99𝑢_𝑘±1−48𝑢_𝑘−51𝑢_𝑘∓12ℎ³+39𝑢^′_𝑘±1+96𝑢^′_𝑘+15𝑢^′_𝑘∓12ℎ², (2.3)(2.4)(2.5)(2.6)(2.7)(2.8) and set𝑓_𝑘±1=𝑓𝑥_𝑘±1,𝑢_𝑘±1,𝑢^′_𝑘±1,𝑢^′′_𝑘±1,𝑢^′′′_𝑘±1,𝑓_𝑘±1=𝑓𝑥_𝑘±1,𝑢_𝑘±1,𝑢^′_𝑘±1,𝑢^′′_𝑘±1,𝑢^′′′_𝑘±1,𝑓_𝑘=𝑓𝑥_𝑘,𝑢_𝑘,𝑢^′_𝑘,𝑢_𝑘^′′,𝑢_𝑘^′′′. (2.9a)(2.9b)(2.9c)Then, the difference method of order four for the differential equation and the corresponding difference method for the derivative 𝑢^′_𝑘 are given by−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴90𝑓_𝑘+1+𝑓_𝑘−1+13𝑓_𝑘,−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1=ℎ⁴60𝑓_𝑘+1−𝑓_𝑘−1. (2.10)(2.11)

Note that 𝑢₀, 𝑢^′₀, 𝑢_𝑁+1 and 𝑢^′_𝑁+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 point 𝑥_𝑘, the given differential equation can be written as𝑢_𝑘^𝑖𝑣=𝑓𝑥_𝑘,𝑢_𝑘,𝑢^′_𝑘,𝑢_𝑘^′′,𝑢_𝑘^′′′=𝑓_𝑘,𝑘=1(1)𝑁. (3.1a) Similarly,𝑓_𝑘±1=𝑓𝑥_𝑘±1,𝑢_𝑘±1,𝑢^′_𝑘±1,𝑢^′′_𝑘±1,𝑢^′′′_𝑘±1,𝑘=1(1)𝑁. (3.1b)

Using the Taylor expansion about the grid point 𝑥_𝑘, we first obtain−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴𝑓90^𝑘+1+𝑓_𝑘−1+13𝑓_𝑘+𝑇₁, (3.2) where 𝑇₁=𝑂(ℎ⁸).

Now, we need the 𝑂(ℎ²) approximation for 𝑢^′′′_𝑘±1. Let𝑢^′′′_𝑘+1=1ℎ³𝑎₁₀𝑢_𝑘+𝑎₁₁𝑢_𝑘+1+𝑎₁₂𝑢_𝑘−1+1ℎ²𝑏₁₀𝑢^′_𝑘+𝑏₁₁𝑢^′_𝑘+1+𝑏₁₂𝑢^′_𝑘−1. (3.3) Expanding each term on the right-hand side of (3.3) in the Taylor series about the point 𝑥_𝑘 and equating the coefficients of ℎ^𝑝 (𝑝=−3,−2,−1,0, and 1)to zero, we get𝑎₁₀,𝑎₁₁,𝑎₁₂,𝑏₁₀,𝑏₁₁,𝑏₁₂=24,−272,−212,0,152,−92. (3.4)

Thus, we obtain𝑢^′′′_𝑘+11=−2ℎ³27𝑢_𝑘+1−48𝑢_𝑘+21𝑢_𝑘−1+12ℎ²15𝑢^′_𝑘+1−9𝑢^′_𝑘−1=𝑢^′′′_𝑘+1−2ℎ²5𝑢^𝑣_𝑘+𝑂ℎ³. (3.5a) Similarly, we obtain𝑢^′′′_𝑘−1=12ℎ³27𝑢_𝑘−1−48𝑢_𝑘+21𝑢_𝑘+1+12ℎ²15𝑢^′_𝑘−1−9𝑢^′_𝑘+1=𝑢^′′′_𝑘−1−2ℎ²5𝑢^𝑣_𝑘−𝑂ℎ³. (3.5b) Further, from (2.3) we obtain:𝑢^′′_𝑘±1=𝑢^′′_𝑘±1−ℎ²3𝑢_𝑘^𝑖𝑣±𝑂ℎ³. (3.6) Hence, we can verify that 𝑓_𝑘±1=𝑓_𝑘±1+𝑂(ℎ²).

Next, we obtain the 𝑂(ℎ⁴) approximation for 𝑢_𝑘^′′. Let𝑢_𝑘^′′=1ℎ²𝑎₂₀𝑢_𝑘+𝑎₂₁𝑢_𝑘+1+𝑎₂₂𝑢_𝑘−1+1ℎ𝑏₂₀𝑢^′_𝑘+𝑏₂₁𝑢^′_𝑘+1+𝑏₂₂𝑢^′_𝑘−1. (3.7) Using the Taylor series expansion and equating the coefficients of ℎ^𝑝 (𝑝=−2,−1,0,1,2, and 3) to zero, we get𝑎₂₀,𝑎₂₁,𝑎₂₂,𝑏₂₀,𝑏₂₁,𝑏₂₂=−4,2,2,0,−12,12. (3.8)

Therefore, 𝑢_𝑘^′′=2ℎ²𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1−1𝑢2ℎ^′𝑘+1−𝑢^′_𝑘−1=𝑢_𝑘^′′+𝑂ℎ⁴. (3.9) Similarly, 𝑢_𝑘^′′′=152ℎ³𝑢_𝑘+1−𝑢_𝑘−1−32ℎ²𝑢^′_𝑘+1+8𝑢^′_𝑘+𝑢^′_𝑘−1=𝑢_𝑘^′′′+𝑂ℎ⁴, (3.10)𝑢^′′_𝑘±1=−12ℎ²11𝑢_𝑘±1−16𝑢_𝑘+5𝑢_𝑘∓1±1ℎ4𝑢^′_𝑘±1−𝑢^′_𝑘∓1=𝑢^′′_𝑘±1∓ℎ³𝑢15^𝑣𝑘+𝑂ℎ⁴. (3.11) Next, we obtain the 𝑂(ℎ³) approximation for 𝑢^′′′_𝑘+1. Let𝑢^′′′_𝑘+1=1ℎ³𝑎₃₀𝑢_𝑘+𝑎₃₁𝑢_𝑘+1+𝑎₃₂𝑢_𝑘−1+1ℎ²𝑏₃₀𝑢^′_𝑘+𝑏₃₁𝑢^′_𝑘+1+𝑏₃₂𝑢^′_𝑘−1. (3.12) Equating the coefficients of ℎ^𝑝 (𝑝=−3,−2,−1,0,1, and 2) to zero, we get𝑎₃₀,𝑎₃₁,𝑎₃₂,𝑏₃₀,𝑏₃₁,𝑏₃₂=24,−992,512,48,392,152. (3.13)

Thus, we obtain𝑢^′′′_𝑘+1=−12ℎ³99𝑢_𝑘+1−48𝑢_𝑘−51𝑢_𝑘−1+12ℎ²39𝑢^′_𝑘+1+96𝑢^′_𝑘+15𝑢^′_𝑘−1=𝑢^′′′_𝑘+1−ℎ³𝑢10^𝑘𝑣𝑖+𝑂ℎ⁴. (3.14) Similarly,𝑢^′′′_𝑘−1=12ℎ³99𝑢_𝑘−1−48𝑢_𝑘−51𝑢_𝑘+1+12ℎ²39𝑢^′_𝑘−1+96𝑢^′_𝑘+15𝑢^′_𝑘+1=𝑢^′′′_𝑘−1+ℎ³𝑢10^𝑘𝑣𝑖+𝑂ℎ⁴. (3.15)

Let 𝛼_𝑘=𝜕𝑓𝜕𝑢_𝑘^′′,𝛽_𝑘=𝜕𝑓𝜕𝑢_𝑘^′′′. (3.16)

From (3.5a), (3.5b), (3.6), and (2.9a) it follows that 𝑓_𝑘±1 provides the 𝑂(ℎ²) approximation for 𝑓_𝑘±1 and𝑓_𝑘±1=𝑓_𝑘±1−ℎ²3𝑢_𝑘^𝑖𝑣𝛼_𝑘−2ℎ²5𝑢^𝑣_𝑘𝛽_𝑘±𝑂ℎ³. (3.17)

Also, from (2.9b), (3.11), (3.14), and (3.15) we obtain the 𝑂(ℎ³) approximation for 𝑓_𝑘±1:𝑓_𝑘±1=𝑓_𝑘±1∓ℎ³𝑢15^𝑣𝑘𝛼_𝑘∓ℎ³𝑢10^𝑘𝑣𝑖𝛽_𝑘+𝑂ℎ⁴. (3.18) From (2.9c), (3.9), and (3.10) it follows that 𝑓_𝑘 provides an 𝑂(ℎ⁴) approximation for 𝑓_𝑘𝑓_𝑘=𝑓_𝑘+𝑂ℎ⁴. (3.19) From (3.18) and (3.19) we get𝑓_𝑘+1+𝑓_𝑘−1+13𝑓_𝑘=𝑓_𝑘+1+𝑓_𝑘−1+13𝑓_𝑘+𝑂ℎ⁴. (3.20) 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 𝑂(ℎ⁸). Similarly, we can verify that the local truncation error associated with the difference scheme (2.11) is of 𝑂(ℎ⁷).

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 𝑢∈𝐶⁸[𝑎,𝑏], and the function 𝑓 is sufficiently differentiable with respect to its arguments. Then, the difference methods (2.10) and (2.11) with the approximations of 𝑢^′′ and 𝑢^′′′ listed in (2.3)–(2.8) are of 𝑂(ℎ⁴).

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) where 𝐀₁₁, 𝐀₁₂, 𝐀₂₁, and 𝐀₂₂ are the 𝑁th-order tri-diagonal matrices and 𝐝₁ and 𝐝₂ 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 𝐀₁₁𝐮^(𝑛+1)=𝜔−𝐀₁₂𝐮^′(𝑛)+𝐝₁+(1−𝜔)𝐀₁₁𝐮^(𝑛)𝐀,𝑛=0,1,2,…,₂₂𝐮^′(𝑛+1)=𝜔−𝐀₂₁𝐮^(𝑛+1)+𝐝₂+(1−𝜔)𝐀₂₂𝐮^′(𝑛),𝑛=0,1,2,…, (3.22) 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.

If 𝑓(𝑥,𝑢(𝑥),𝑣(𝑥),𝑢^′(𝑥),𝑣^′(𝑥))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𝚽()𝚿()𝐮,𝐯=0,𝐮,𝐯=0, (3.23) where𝑢𝐮=¹𝑢₂⋮𝑢_𝑁𝑣,𝐯=¹𝑣₂⋮𝑣_𝑁()=𝜙,𝚽𝐮,𝐯¹()𝜙𝐮,𝐯₂()⋮𝜙𝐮,𝐯_𝑁()()=𝜓𝐮,𝐯,𝚿𝐮,𝐯¹()𝜓𝐮,𝐯₂()⋮𝜓𝐮,𝐯_𝑁(𝐮,𝐯). (3.24)

Let𝐓𝐉=¹¹𝐓₁₂𝐓₂₁𝐓₂₂ (3.25)

be the Jacobian of Φ and Ψ, which is the 2𝑁th-order block tridiagonal matrix, where T₁₁=𝜕𝜙₁,𝜙₂,…,𝜙_𝑁𝜕𝑢₁,𝑢₂,…,𝑢_𝑁=𝜕𝜙₁𝜕𝑢₁𝜕𝜙₁𝜕𝑢₂𝜕𝜙₂𝜕𝑢₁𝜕𝜙₂𝜕𝑢₂𝜕𝜙₂𝜕𝑢₃𝟎0𝜕𝜙_𝑁𝜕𝑢_𝑁−1𝜕𝜙_𝑁𝜕𝑢_𝑁,T₁₂=𝜕𝜙₁,𝜙₂,…,𝜙_𝑁𝜕𝑣₁,𝑣₂,…,𝑣_𝑁,T₂₁=𝜕𝜓₁,𝜓₂,…,𝜓_𝑁𝜕𝑢₁,𝑢₂,…,𝑢_𝑁,T₂₂=𝜕𝜓₁,𝜓₂,…,𝜓_𝑁𝜕𝑣₁,𝑣₂,…,𝑣_𝑁 (3.26)

are the 𝑁th-order tridiagonal matrices.

In the NBSOR method starting with any initial approximation (u⁽⁰⁾,v⁽⁰⁾) of (𝐮^(𝑘),𝐯^(𝑘)), 𝑘=1(1)𝑁, we define𝐮^(𝑛+1)=𝐮^(𝑛)+Δ𝐮^(𝑛)𝐯,𝑛=0,1,2,…,_(𝑛+1)=𝐯^(𝑛)+Δ𝐯^(𝑛),𝑛=0,1,2,…, (3.27) where Δ𝐮 and Δ𝐯 are intermediate values obtained by solving the matrix equation for NBSOR method given by𝐓₁₁𝐓₁₂𝐓₂₁𝐓₂₂=Δ𝐮Δ𝐯−𝚽−𝚿. (3.28) The above system can be solved for Δ𝐮^(𝑛) and Δ𝐯^(𝑛) by using the block SOR method (inner iterative method) as follows:𝐓₁₁Δ𝐮^(𝑛+1)=𝜔−𝚽𝐮^(𝑛),𝐯^(𝑛)−𝐓₁₂Δ𝐯^(𝑛)+(1−𝜔)𝐓₁₁Δ𝐮^(𝑛)𝐓,𝑛=0,1,2,…,₂₂Δ𝐯^(𝑛+1)=𝜔−𝚿𝐮^(𝑛),𝐯^(𝑛)−𝐓₂₁Δ𝐮^(𝑛+1)+(1−𝜔)𝐓₂₂Δ𝐯^(𝑛),𝑛=0,1,2,…, (3.29) where 𝜔∈(0,2) is a relaxation parameter and 𝑛=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 (𝐮⁽⁰⁾,𝐯⁽⁰⁾) be close to the solution.

4. Convergence and Stability Analysis

Consider the model problem 𝑢^𝑖𝑣=𝑓(𝑥), where 𝑓 is a function of 𝑥 only. Applying the fourth-order difference methods (2.10) and (2.11) to the above equation, we get𝑢_𝑘−1−2𝑢_𝑘+𝑢_𝑘+1+ℎ2𝑣_𝑘−1−𝑣_𝑘+1=−ℎ⁴𝑓180^𝑘+1+𝑓_𝑘−1+13𝑓_𝑘,3ℎ𝑢_𝑘−1−𝑢_𝑘+1+𝑣_𝑘−1+4𝑣_𝑘+𝑣_𝑘+1=ℎ³𝑓60^𝑘+1−𝑓_𝑘−1, (4.1) where we denote 𝑢′=𝑣.

Let us denote by 𝐏=[1,0,1], 𝐋=[1,1,1], and 𝐌=[1,0,−1] the 𝑁th-order tridiagonal matrices.

The system of (4.1) can be written in the block form as ()h𝐋−3𝐈2𝐌3h𝐌()𝐮𝐯=𝐝𝐋+3𝐈𝐝¹₂, (4.2) where 𝐮 and 𝐯 are two 𝐍-dimensional solution vectors and 𝐝₁, 𝐝₂ are vectors consisting of right-hand side functions and some boundary values associated with (4.1). The BSOR method for the scheme is𝐮^(𝑛+1)=(1−𝜔)𝐮^(𝑛)−𝜔ℎ2(𝐋−3𝐈)⁻¹𝐌(𝐯)^(𝑛)+𝜔(𝐋−3𝐈)⁻¹𝐝₁,(𝐯)^(𝑛+1)=(1−𝜔)(𝐯)^(𝑛)−3𝜔ℎ(𝐋+3𝐈)⁻¹𝐌𝐮^(𝑛+1)+𝜔(𝐋+3𝐈)⁻¹𝐝₂, (4.3) The associated block SOR and block Jacobi iteration matrices are given by𝐋_𝜔=(1−𝜔)𝐈−𝜔ℎ2(𝐋−3𝐈)⁻¹𝐌−3𝜔ℎ(𝐋+3𝐈)⁻¹𝐌()𝐈,01−𝜔𝐁=−ℎ2(𝐋−3𝐈)⁻¹𝐌−3ℎ(𝐋+3𝐈)⁻¹,𝐌0 (4.4) and 𝜆 and 𝜂 are the eigenvalues associated with the corresponding matrices 𝐋_𝜔 and 𝐁, which are related by the equation(𝜆+𝜔−1)²=𝜆𝜔²𝜂². (4.5) Let 𝐯₁𝐯₂ be the eigenvector associated with the eigenvalue 𝜂 so that0−ℎ2(𝐋−3𝐈)⁻¹𝐌−3ℎ(𝐋+3𝐈)⁻¹𝐯𝐌0𝐯¹₂=𝐯₁𝐯₂𝜂. (4.6) That is,−ℎ2(𝐋−3𝐈)⁻¹𝐌𝐯₂=𝜂𝐯₁,−3ℎ(𝐋+3𝐈)⁻¹𝐌𝐯₁=𝜂𝐯₂. (4.7) On eliminating 𝐯₂, we obtain32(𝐋−3𝐈)⁻¹𝐌(𝐋+3𝐈)⁻¹M𝐯₁=𝜂²𝐯₁. (4.8) The rate of convergence of the BSOR method is given by −𝐥𝐨𝐠(𝜌(𝐋_𝜔)).

The rate of convergence of the BSOR method is dependent on the eigenvalues of B (through relation (4.5)), which are given by (3𝜏/2)=𝜂², where 𝜏 are the eigenvalues of (𝐋−3𝐈)⁻¹𝐌(𝐋+3𝐈)⁻¹𝐌.

Hence, we can determine the optimal parameter as 𝜔₀=21+1−(3/2)𝜏, (4.9) where 𝜏=𝑆((𝐋−3𝐈)⁻¹𝐌(𝐋+3𝐈)⁻¹𝐌).

Thus, we can determine the convergence factor 𝜆=𝜔₀−1=1−1−(3/2)𝜏1+1−(3/2)𝜏. (4.10) For convergence, we must have |𝜆|<1 to give the range 0<2𝜏<3 (4.11)

Hence, 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 𝑢^𝑖𝑣=𝑓(𝑥) converges if 0<𝜏<2/3, where 𝜏=𝑆((𝐋−3𝐈)⁻¹𝐌(𝐋+3𝐈)⁻¹𝐌), 𝐋=[1,1,1],and𝐌=[1,0,−1] are the 𝑁×𝑁 tridiagonal matrices.

Now, we discuss the stability analysis.

An iterative method for (4.1) can be written as 𝐮^(𝑘+1)=12𝐏𝐮^(𝑘)+ℎ4𝐌𝐯^(𝑘)𝐯+𝐑𝐇𝐔,_(𝑘+1)=−34ℎ𝐌𝐮^(𝑘)−14𝐏𝐯^(𝑘)+𝐑𝐇𝐕, (4.12) where 𝐮^(𝑘), 𝐯^(𝑘) are solution vectors at the 𝑘th 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 𝐮^(𝑘+1)𝐯^(𝑘+1)𝐮=𝐆^(𝑘)𝐯^(𝑘)+𝐑𝐇, (4.13) where1𝐆=2𝐏ℎ4𝐌−3𝐌4ℎ−14𝐏,𝐑𝐇=RHURHV. (4.14)

The eigenvalues of 𝐏 and 𝐌 are 2cos𝑘𝜋/(𝑁+1) and 2𝑖cos𝑘𝜋/(𝑁+1), respectively, where 𝑘=1,2,…,𝑁. The characteristic equation of the matrix 𝐆 is given by1det2𝐏−𝜉𝐈_𝑁ℎ4𝐌−3𝐌4ℎ−14𝐏−𝜉𝐈_𝑁=0. (4.15)

Thus the eigenvalues of 𝐆 are given by det−14𝐏−𝜉𝐈_𝑁1×det2𝐏−𝜉𝐈_𝑁+3𝐌16−14𝐏−𝜉𝐈_𝑁⁻¹𝐌=0. (4.16) 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Δ⁴𝑑𝑢≡²𝑑𝑟²+1𝑟𝑑𝑑𝑟²𝑢=𝑓(𝑟),0<𝑟<1, (5.1) or equivalently𝑢^𝑖𝑣=𝑏(𝑟)𝑢^′′′+𝑐(𝑟)𝑢^′′+𝑑(𝑟)𝑢^′+𝑓(𝑟),0<𝑟<1, (5.2) where 𝑏(𝑟)=−2/𝑟,𝑐(𝑟)=1/𝑟²,𝑑(𝑟)=−1/𝑟³. (5.3)

The above equation represents fourth-order ordinary differential equation in cylindrical polar coordinates.

The boundary conditions are given by𝑢(0)=𝐴₀,𝑢^′(0)=𝐴₁,𝑢(1)=𝐵₀,𝑢^′(1)=𝐵₁, (5.4) where 𝐴₀, 𝐴₁, 𝐵₀, and 𝐵₁ are constants.

Applying the difference scheme (2.2a), (2.2b) to the singular equation (5.2), we obtain a second-order difference method−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴6𝑏_𝑘𝑢_𝑘^′′′+𝑐_𝑘𝑢_𝑘^′′+𝑑_𝑘𝑢^′_𝑘+𝑓_𝑘,−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1=0,𝑘=1(1)𝑁. (5.5) Applying the fourth-order difference scheme (2.10) to the singular equation (5.2), we obtain−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴𝑏90^𝑘+1𝑢^′′′_𝑘+1+𝑐_𝑘+1𝑢^′′_𝑘+1+𝑑_𝑘+1𝑢^′_𝑘+1+𝑓_𝑘+1+𝑏_𝑘−1𝑢^′′′_𝑘−1+𝑐_𝑘−1𝑢^′′_𝑘−1+𝑑_𝑘−1𝑢^′_𝑘−1+𝑓_𝑘−1+13𝑏_𝑘𝑢_𝑘^′′′+13𝑐_𝑘𝑢_𝑘^′′+13𝑑_𝑘𝑢^′_𝑘+13𝑓_𝑘,𝑘=1(1)𝑁, (5.6) where 𝑏_𝑘±1=𝑏𝑟_𝑘±1,𝑐_𝑘±1=𝑐𝑟_𝑘±1,𝑑_𝑘±1=𝑑𝑟_𝑘±1,𝑓_𝑘±1=𝑓𝑟_𝑘±1. (5.7)

Note that the scheme fails when the solution is to be determined at 𝑘=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 𝑟=0 (see [29]). We consider the following approximations:𝑏_𝑘±1=𝑏_𝑘±ℎ𝑏^′_𝑘+ℎ²𝑏2!^𝑘′′+𝑂±ℎ³,𝑐_𝑘±1=𝑐_𝑘±ℎ𝑐^′_𝑘+ℎ²𝑐2!^𝑘′′+𝑂±ℎ³,𝑑_𝑘±1=𝑑_𝑘±ℎ𝑑^′_𝑘+ℎ²𝑑2!^𝑘′′+𝑂±ℎ³,𝑓_𝑘±1=𝑓_𝑘±ℎ𝑓^′_𝑘+ℎ²𝑓2!^𝑘′′+𝑂±ℎ³. (5.8)

Using the approximation (5.8) in (5.6) and neglecting higher-order terms, we can rewrite (5.6) in compact operator form as−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ²90−24𝑏^′_𝑘+18𝑐_𝑘−4ℎ²𝑐_𝑘^′′𝛿²_𝑥𝑢_𝑘+ℎ18045𝑏_𝑘−75ℎ²𝑏_𝑘^′′−6ℎ²𝑐^′_𝑘2𝜇_𝑥𝛿_𝑥𝑢_𝑘+ℎ⁴90−45𝑏_𝑘ℎ²+75𝑏_𝑘^′′+6𝑐^′_𝑘+15𝑑_𝑘+ℎ²𝑑_𝑘^′′𝑢^′_𝑘+ℎ²18015𝑏_𝑘+27ℎ²𝑏_𝑘^′′+6ℎ²𝑐^′_𝑘+2ℎ²𝑑_𝑘2𝜇_𝑥𝛿_𝑥𝑢^′_𝑘+ℎ³18024𝑏^′_𝑘−3𝑐_𝑘+5ℎ²𝑐_𝑘^′′+2ℎ²𝑑^′_𝑘2𝜇_𝑥𝛿_𝑥𝑢_𝑘+ℎ⁴9015𝑓_𝑘+ℎ²𝑓_𝑘^′′,𝑘=1(1)𝑁. (5.9)

Similarly, using the difference scheme (2.11), a fourth-order approximation for the derivative 𝑢^′ for the singular equation (5.2) in the compact form may be written as−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1=ℎ60−24𝑏_𝑘𝛿²_𝑥𝑢_𝑘+ℎ²60−3𝑏^′_𝑘2𝜇_𝑥𝛿_𝑥𝑢_𝑘+ℎ³602𝑐_𝑘+3𝑏^′_𝑘𝛿²_𝑥𝑢^′_𝑘+ℎ²6012𝑏_𝑘+ℎ²𝑐^′_𝑘+𝑑_𝑘2𝜇_𝑥𝛿_𝑥𝑢^′_𝑘+ℎ³606𝑏^′_𝑘+2ℎ²𝑑^′_𝑘𝑢^′_𝑘+ℎ⁴602ℎ𝑓^′_𝑘,𝑘=1(1)𝑁, (5.10) where 𝛿_𝑟𝑢_𝑘=𝑢_𝑘+1/2−𝑢_𝑘−1/2, 𝜇_𝑟𝑢_𝑘=(𝑢_𝑘+1/2+𝑢_𝑘−1/2)/2.

Finite difference equations (5.9) and (5.10) along with the boundary conditions (5.4) gives a 2𝑁×2𝑁 linear system of equations for the unknowns 𝑢₁,𝑢₂,…,𝑢_𝑁,𝑢^′₁,𝑢^′₂,…,𝑢^′_𝑁. 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/(𝑘±1), hence very easily solved for 𝑘=1(1)𝑁 in the region (0,1).

Consider the coupled nonlinear singular equations𝑢^𝐼𝑉=𝑎(𝑟)𝑢^′𝑣^′′+𝑣^′𝑢^′′(𝑟)𝑣+𝑓,0<𝑟<1,_𝐼𝑉=−𝑎(𝑟)𝑢^′𝑢^′′+𝑔(𝑟),0<𝑟<1, (5.11) where 𝑎(𝑟)=1/𝑟, with known boundary conditions 𝑢(0), 𝑣(0), 𝑢^′(0), 𝑣^′(0), 𝑢(1), 𝑣(1), 𝑢^′(1), and 𝑣^′(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−2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=ℎ⁴6𝑎_𝑘𝑢^′_𝑘𝑣_𝑘^′′+𝑣^′_𝑘𝑢_𝑘^′′+𝑓_𝑘,−2𝑣_𝑘+1−2𝑣_𝑘+𝑣_𝑘−1+ℎ𝑣^′_𝑘+1−𝑣^′_𝑘−1=ℎ⁴6−𝑎_𝑘𝑢^′_𝑘𝑢_𝑘^′′+𝑔_𝑘,−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1𝑣=0,−3_𝑘+1−𝑣_𝑘−1+ℎ𝑣^′_𝑘+1+4𝑣^′_𝑘+𝑣^′_𝑘−1=0, (5.12) where 𝑎_𝑘=𝑎(𝑟_𝑘), 𝑓_𝑘=𝑓(𝑟_𝑘), and 𝑔_𝑘=𝑔(𝑟_𝑘).

The fourth-order difference scheme for solving the system (5.11) for 𝑢, 𝑣, 𝑢^′, and 𝑣^′ is given by −2𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1+ℎ𝑢^′_𝑘+1−𝑢^′_𝑘−1=𝐻₁15𝑓_𝑘+ℎ²𝑓_𝑘^′′+𝐻₁𝑎₁₁𝑢^′_𝑘+1−12ℎ²11𝑣_𝑘+1−16𝑣_𝑘+5𝑣_𝑘−1+1ℎ4𝑣^′_𝑘+1−𝑣^′_𝑘−1+𝑣^′_𝑘+1−12ℎ²11𝑢_𝑘+1−16𝑢_𝑘+5𝑢_𝑘−1+1ℎ4𝑢^′_𝑘+1−𝑢^′_𝑘−1+𝐻₁𝑎₁₂𝑢^′_𝑘−1−12ℎ²11𝑣_𝑘−1−16𝑣_𝑘+5𝑣_𝑘+1−1ℎ4𝑣^′_𝑘−1−𝑣^′_𝑘+1+𝑣^′_𝑘−1−12ℎ²11𝑢_𝑘−1−16𝑢_𝑘+5𝑢_𝑘+1−1ℎ4𝑢^′_𝑘−1−𝑢^′_𝑘+1+𝐻₁𝑎₁₀𝑢^′_𝑘2ℎ²𝑣_𝑘+1−2𝑣_𝑘+𝑣_𝑘−1−1𝑣2ℎ^′𝑘+1−𝑣^′_𝑘−1+𝑣^′_𝑘2ℎ²𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1−1𝑢2ℎ^′𝑘+1−𝑢^′_𝑘−1,−2𝑣_𝑘+1−2𝑣_𝑘+𝑣_𝑘−1+ℎ𝑣^′_𝑘+1−𝑣^′_𝑘−1=𝐻₁15𝑔_𝑘+ℎ²𝑔_𝑘^′′−𝐻₁𝑎₁₁𝑢^′_𝑘+1−12ℎ²11𝑢_𝑘+1−16𝑢_𝑘+5𝑢_𝑘−1+1ℎ4𝑢^′_𝑘+1−𝑢^′_𝑘−1−𝐻₁𝑎₁₂𝑢^′_𝑘−1−12ℎ²11𝑢_𝑘−1−16𝑢_𝑘+5𝑢_𝑘+1−1ℎ4𝑢^′_𝑘−1−𝑢^′_𝑘+1−13𝐻₁𝑎₁₀𝑢^′_𝑘2ℎ²𝑢_𝑘+1−2𝑢_𝑘+𝑢_𝑘−1−1𝑢2ℎ^′𝑘+1−𝑢^′_𝑘−1,−3𝑢_𝑘+1−𝑢_𝑘−1+ℎ𝑢^′_𝑘+1+4𝑢^′_𝑘+𝑢^′_𝑘−1=𝐻₂2ℎ𝑓^′_𝑘+𝐻₂𝑎₁₁𝑢^′_𝑘+13𝑣^′_𝑘+1−4𝑣^′_𝑘+𝑣^′_𝑘−12ℎ+𝑣′_𝑘+13𝑢^′_𝑘+1−4𝑢^′_𝑘+𝑢^′_𝑘−12ℎ−𝐻₂𝑎₁₂𝑢^′_𝑘−1−3𝑣^′_𝑘−1+4𝑣^′_𝑘−𝑣^′_𝑘+12ℎ+𝑣^′_𝑘−1−3𝑢^′_𝑘−1+4𝑢^′_𝑘−𝑢^′_𝑘+1,𝑣2ℎ−3_𝑘+1−𝑣_𝑘−1+ℎ𝑣^′_𝑘+1+4𝑣^′_𝑘+𝑣^′_𝑘−1=𝐻₂2ℎ𝑔^′_𝑘−𝐻₂𝑎₁₁𝑢^′_𝑘+13𝑢^′_𝑘+1−4𝑢′_𝑘+𝑢^′_𝑘−12ℎ−𝐻₂𝑎₁₂𝑢^′_𝑘−1−3𝑢^′_𝑘−1+4𝑢′_𝑘−𝑢^′_𝑘+1,2ℎ (5.13) where𝑎_𝑘=𝑎𝑟_𝑘,𝑓_𝑘=𝑓𝑟_𝑘,𝑔_𝑘=𝑔𝑟_𝑘,𝑎₁₁=𝑎_𝑘+ℎ𝑎^′_𝑘+ℎ²2𝑎_𝑘^′′,𝑎₁₂=𝑎_𝑘−ℎ𝑎^′_𝑘+ℎ²2𝑎_𝑘^′′,𝑎₁₀=𝑎_𝑘. (5.14)

The scheme (5.13) is free from the terms 1/𝑘±1, hence very easily solved for 𝑘=1(1)𝑁 in the region (0,1). The system (5.13) can be solved using the NBSOR method.

Consider the boundary value problem [31]𝑦^𝑖𝑣−𝜆𝑦𝑦^′′(𝑥),𝑦(0)=𝑓=𝑉₀,𝑦(1)=𝑉₁,𝑦^′(0)=0,𝑦^′(1)=0. (5.15)

This arises from the time-dependent Navier-Strokes equations for axisymmetric flow of an incompressible fluid contained between infinite disks that occupy the planes 𝑧=−𝑑 and 𝑧=𝑑. The disks are porous and the fluid is injected or extracted normally with velocity 𝑉₀at𝑧=−𝑑 and 𝑉₁ at 𝑧=𝑑. Here, =𝑑/𝜈, 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⁻¹². 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 𝑢=𝑟⁴sin𝑟. The root mean square errors (RMSEs) are tabulated in Table 1. The graph of the errors for 𝑁=32 is given in Figure 1.

Table 1 

The RMSE.

Figure 1 

Problem 2. The boundary value problem is to solve (5.15). The exact solution is (1−𝑥²)exp(𝑥). The maximum absolute errors (MAE) and RMSE are tabulated in Table 2.

Table 2 

The MAE and RMSE.

Problem 3. The system of nonlinear equation (5.11) is to be solved subject to the natural boundary conditions. The exact solutions are 𝑢=cos(𝑟) and 𝑣=exp(𝑟). The MAE and RMSE are tabulated in Table 3.

Table 3 

The MAE and RMSE.

Observe that for the linear singular problem for 𝛼=1, ℎ₁=1/64, ℎ₂=1/128, log(𝑒_ℎ1/𝑒_ℎ2)/log(ℎ₁/ℎ₂)=log(0.6578(−07)/0.3345(−08))/log2≈4. Thus, we obtain fourth-order convergence for 𝑢. Similarly, we get fourth-order convergence for 𝑢^′.

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.

Acknowledgments

This research was supported by “The Council of Scientific and Industrial Research” under Research Grant no. 09/045(0836)2009-EMR-I. The authors thank the reviewers for their valuable suggestions, which substantially improved the standard of the paper.