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

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


Introduction
Consider the fourth-order boundary value problem u iv x f x, u x , u x , u x , u x , a < x < b, 1 where A 0 , A 1 , B 0 , and B 1 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 : 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][18][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 fourthorder 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, x k kh, k 1, 2, . . . , N. Let u k and u k denote the approximate solutions, and let U k u x k and U k u x k be the exact solution values of u x and u x at the grid point x x k , respectively. Also, x 0 a and x N 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.

The Finite Difference Method
The method is described as follows: for k 1 1 N, let

2.1
Then, the difference method of order two for the given differential 1.1 is given by, and the corresponding difference method for the derivative u x is given by Also, let Then, the difference method of order four for the differential equation and the corresponding difference method for the derivative u k are given by Note that u 0 , u 0 , u N 1 and u N 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 .

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 x k , the given differential equation can be written as Similarly, Using the Taylor expansion about the grid point x k , we first obtain

5
Now, we need the O h 2 approximation for u k±1 . Let Expanding each term on the right-hand side of 3.3 in the Taylor series about the point x k and equating the coefficients of h p p −3, −2, −1, 0, and 1 to zero, we get Thus, we obtain

3.5a
Similarly, we obtain

3.5b
Further, from 2.3 we obtain: Hence, we can verify that f k±1 f k±1 O h 2 . Next, we obtain the O h 4 approximation for u k . Let Using the Taylor series expansion and equating the coefficients of h p p −2, −1, 0, 1, 2, and 3 to zero, we get Advances in Numerical Analysis Therefore, Similarly, Next, we obtain the O h 3 approximation for u k 1 . Let Equating the coefficients of h p p −3, −2, −1, 0, 1, and 2 to zero, we get Thus, we obtain

3.15
Let Advances in Numerical Analysis 7 From 3.5a , 3.5b , 3.6 , and 2.9a it follows that f k±1 provides the O h 2 approximation for f k±1 and Also, from 2.9b , 3.11 , 3.14 , and 3.15 we obtain the O h 3 approximation for f k±1 : From 2.9c , 3.9 , and 3.10 it follows that f k provides an O h 4 approximation for f k From 3.18 and 3.19 we get 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. If the differential 1.1 is linear, then the difference method 2.10 and 2.11 in the matrix form can be written as where A 11 , A 12 , A 21 , and A 22 are the Nth-order tri-diagonal matrices and d 1 and d 2 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 − ω A 11 u n , n 0, 1, 2, . . . ,

Advances in Numerical Analysis
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 f 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.24
Let J ⎡ ⎣ T 11 T 12 be the Jacobian of Φ and Ψ, which is the 2Nth-order block tridiagonal matrix, where are the Nth-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 u n 1 u n Δu n , n 0, 1, 2, . . . , v n 1 v n Δv n , n 0, 1, 2, . . . , 3.27 where Δu and Δv are intermediate values obtained by solving the matrix equation for NBSOR method given by 3.28 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
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.

Convergence and Stability Analysis
Consider the model problem u iv f x , where f is a function of x only. Applying the fourthorder difference methods 2.10 and 2.11 to the above equation, we get where we denote u v.
The system of 4.1 can be written in the block form as where u and v are two N-dimensional solution vectors and d 1 , d 2 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
The associated block SOR and block Jacobi iteration matrices are given by and λ and η are the eigenvalues associated with the corresponding matrices L ω and B, which are related by the equation Let v 1 v 2 be the eigenvector associated with the eigenvalue η so that That is,

4.7
Advances in Numerical Analysis

11
On eliminating v 2 , we obtain 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 of B through relation 4.5 , which are given by 3τ/2 η 2 , where τ are the eigenvalues of Hence, we can determine the optimal parameter as ω 0 2 Thus, we can determine the convergence factor 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. Now, we discuss the stability analysis. An iterative method for 4.1 can be written as

4.12
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 Advances in Numerical Analysis The eigenvalues of P and M are 2 cos kπ / N 1 and 2i cos kπ/ N 1 , respectively,  where k 1, 2, . . . , N. The characteristic equation of the matrix G is given by Thus the eigenvalues of G are given by 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.

Application to Singular Equation
Consider a singular fourth-order linear ordinary differential equation of the form The above equation represents fourth-order ordinary differential equation in cylindrical polar coordinates.
The boundary conditions are given by where A 0 , A 1 , B 0 , and B 1 are constants.

13
Applying the difference scheme 2.2a , 2.2b to the singular equation 5.2 , we obtain a second-order difference method

5.5
Applying the fourth-order difference scheme 2.10 to the singular equation 5.2 , we obtain Note that the scheme fails when the solution is to be determined at k 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
Using the approximation 5.8 in 5.6 and neglecting higher-order terms, we can rewrite 5.6 in compact operator form as

Advances in Numerical Analysis
15f k h 2 f k , k 1 1 N.

5.9
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
where δ r u k u k 1/2 − u k−1/2 , μ r u k u k 1/2 u k−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 u 1 , u 2 , . . . , u N , u 1 , u 2 , . . . , u N . 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 u IV a r u v v u f r , 0 < r < 1, v IV −a r u u g r , 0 < r < 1,

5.11
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 where a k a r k , f k f r k , and g k g r k .
The fourth-order difference scheme for solving the system 5.11 for u, v, u , and v is given by where a k a r k , f k f r k , g k g r k , a 11 a k ha k h 2 2 a k , a 12 a k − ha k h 2 2 a k , a 10 a k .

5.14
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
This arises from the time-dependent Navier-Strokes equations for axisymmetric flow of an incompressible fluid contained between infinite disks that occupy the planes z −d and z d. The disks are porous and the fluid is injected or extracted normally with velocity V 0 at z −d and V 1 at z d. Here, d/ν, where ν is the kinematic viscosity.

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.  Table 1. The graph of the errors for N 32 is given in Figure 1.  Table 2.
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.

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.