Finite Difference Method for Solving a System of Third-Order Boundary Value Problems

Abstract

We develop a new-two-stage finite difference method for computing approximate solutions of a system of third-order boundary value problems associated with odd-order obstacle problems. Such problems arise in physical oceanography (Dunbar (1993) and Noor (1994), draining and coating flow problems (E. O. Tuck (1990) and L. W. Schwartz (1990)), and can be studied in the framework of variational inequalities. We show that the present method is of order three and give numerical results that are better than the other available results. Numerical example is presented to illustrate the applicability and efficiency of the new method.

1. Introduction

Variational inequalities have had a great impact and influence in the development of almost all branches of pure and applied sciences. It has been shown that the variational inequalities provide a novel and general framework to study a wide class of problems arising in various branches of pure and applied sciences. The ideas and techniques of variational inequalities are being used in a variety of diverse fields and proved to be innovative and productive, see [1–15] and the references therein. In recent years, variational inequalities have been extended and generalized in several directions. A useful and important generalization of variational inequalities is called the general variational inequalities involving two continuous operators, which was introduced by Noor [10] in 1988. It has been shown that a wide class of nonsymmetric and odd-order obstacle problems arising in industry, economics, optimization, mathematical and engineering sciences can be studied in the unified and general framework of the general variational inequalities, see [1–5, 11–19] and the references therein. Despite of their importance, little attention has been given to develop some efficient numerical technique of solving such type of problems. In principle, the finite difference techniques and other related methods cannot be applied directly to solve the obstacle-type problems. Using the technique of the penalty method, one rewrite the general variational inequalities as general variational equations. We note that, if the obstacle function is known, then one can use the idea and technique of Lewy and Stampacchia [9] to express the general variational equations as a system of third-order boundary value problems. This resultant system of equations can be solved, which is the main advantage of this approach. The computational aspect of this method is its simple applicability for solving obstacle problems. Such type of penalty function method in conjunction with spline and finite difference techniques has been used quite effectively as a basis for solving system of third-order boundary value problems, see [1–5, 7, 13, 15–17, 19]. Our approach to these problems is to consider them in a general manner and specialize them later on. To convey an idea of the technique involved, we first introduce and develop a new two stage finite difference scheme for solving a third-order boundary value problems. An example involving the order-order obstacle problem is given to illustrate the efficiency and its comparison with other methods.

For simplicity and to convey an idea of the obstacle problems, we consider a system of third-order boundary value problem of the type, which was first considered by Noor [11] 𝑢^′′′=𝑓(𝑥)𝑝(𝑥)𝑢(𝑥)(𝑥)𝑓(𝑥),𝑎≤𝑥≤𝑐,+𝑓+𝑟,𝑐≤𝑥≤𝑑,,𝑑≤𝑥≤𝑏, (1.1a) with the boundary conditions𝑢(𝑎)=𝛼,𝑢^′(𝑎)=𝛽₁,𝑢^′(𝑏)=𝛽₂, (1.1b)and the continuity conditions of 𝑢, 𝑢^′, and 𝑢^′′ at 𝑐 and 𝑑. Here, 𝑓 and 𝑝 are continuous functions on [𝑎,𝑏] and [𝑐,𝑑], respectively. The parameters 𝑟, 𝛼, 𝛽₁, and 𝛽₂ are real finite constants. Using the penalty method technique, one can easily show that a wide class of unrelated obstacle, unilateral, moving and free boundary value problems arising in various branches of pure and applied sciences can be characterized by the system of third-order boundary value problems of type (1.1a) and (1.1b), see, for example, [1–17] and the references therein. In general, it is not possible to obtain the analytical solution of (1.1a) and (1.1b) for arbitrary choices of 𝑓(𝑥) and 𝑝(𝑥). We usually resort to some numerical methods for obtaining an approximate solution of (1.1a) and (1.1b).

The available finite difference and collocation methods are not suitable for solving system of boundary value problems of the form defined by (1.1a) and (1.1b). Such methods have a serious drawback in the accuracy regardless of the order of the convergent of the method being used, see [2, 4, 7, 13, 16, 17, 19]. On the other hand, Al-Said [1], Al-Said and Noor [3], Al-Said et al. [5] and Noor and Al-Said [16] have developed first- and second-order two-stage difference methods for solving (1.1a) and (1.1b), which gives numerical results that better than those produced by the first-, second- and, third-order methods considered in [2, 4, 7, 13, 16, 17, 19].

Motivated by the above works, we suggest a new two-stage numerical algorithm for solving the system of third-order boundary value problem (1.1a) and (1.1b). We prove that the present method is of order two, and it outperforms other collocation and finite difference methods when solving (1.1a) and (1.1b). In Section 2, we derive the numerical method for solving (1.1a) and (1.1b). Section 3 is devoted for the convergence analysis of the method. The numerical experiments and comparison with other methods are given in Section 4.

2. Numerical Method

For simplicity, we take 𝑐=(3𝑎+𝑏)/4 and 𝑑=(𝑎+3𝑏)/4 in order to develop the numerical method for solving the system of differential equations (1.1a) and (1.1b). For this purpose we divide the interval [𝑎,𝑏] into 𝑛 equal subintervals using the grid points 𝑥_𝑖=𝑎+𝑖ℎ,𝑖=0,1,2,…,𝑛, 𝑥₀=𝑎, 𝑥_𝑛=𝑏 andℎ=𝑏−𝑎𝑛, (2.1) where 𝑛 is a positive integer chosen such that both 𝑛/4 and 3𝑛/4 are also positive integers.

Using Taylor series expansions along with the method of undetermined coefficients, boundary and continuity conditions to develop the following finite difference scheme:9𝑢_1/2−𝑢_3/2=8𝑢₀+3ℎ𝑢^′₀−1ℎ160³6𝑢₀^′′′+51𝑢^′′′_1/2+3𝑢^′′′_3/2+𝑡₁,for𝑖=1,−2𝑢_1/2+3𝑢_3/2−𝑢_5/2=+ℎ𝑢^′₀−1ℎ1920³809𝑢^′′′_1/2+1062𝑢^′′′_3/2−31𝑢^′′′_5/2+𝑡₂𝑢,for𝑖=2,_𝑖−5/2−3𝑢_𝑖−3/2+3𝑢_𝑖−1/2−𝑢_𝑖+1/2=16ℎ³−𝑢^′′′_𝑖−5/2+6𝑢^′′′_𝑖−3/2+𝑢^′′′_𝑖+1/2+𝑡_𝑖𝑢,for3≤𝑖≤𝑛−1,_𝑛−5/2−3𝑢_𝑛−3/2+2𝑢_𝑛−1/2=ℎ𝑢^′_𝑛−1ℎ1920³−31𝑢^′′′_𝑛−5/2+1062𝑢^′′′_𝑛−3/2+809𝑢^′′′_𝑛−1/2+𝑡_𝑛,for𝑖=𝑛, (2.2) where𝑢^′′′_𝑖+1/2=𝑓_𝑖+1/2𝑛,for0≤𝑖≤4−1,3𝑛4𝑝≤𝑖≤𝑛−1_𝑖+1/2𝑢_𝑖+1/2+𝑓_𝑖+1/2𝑛+𝑟,for4≤𝑖≤3𝑛4−1, (2.3)𝑓_𝑖+1/2=𝑓(𝑥_𝑖+1/2),𝑖=0,1,2,…,𝑛−1, the relation (2.2) forms a system of 𝑛 linear equations in the unknowns 𝑢_𝑖−1/2,𝑖=1,2,…,𝑛. The local truncation errors associate with (2.2) are given by𝑡_𝑖=49ℎ320⁶𝑢₀^(𝑣𝑖)+𝑂ℎ⁷−,for𝑖=1,53ℎ3840⁶𝑢₀^(𝑣𝑖)+𝑂ℎ⁷,for𝑖=2,1523ℎ7680⁶𝑢^(𝑣𝑖)_i+𝑂ℎ⁷−,for3≤𝑖≤𝑛−1,53ℎ3840⁶𝑢_𝑖^(𝑣𝑖)+𝑂ℎ⁷,for𝑖=𝑛, (2.4) which suggest that the scheme (2.2) is a third-order accurate.

Remark 2.1. The new method can be considered as an improvement of the previous finite difference methods at the midknots developed in [1, 3, 4, 16] for solving the third-order obstacle problem. Thus, the matrix remains the same for the sake of comparison with other methods. We would like to point out that the same goes for the finite difference methods at the knots for solving third-order boundary value problems, see the references.

3. Convergence Analysis

In this section, we investigate the convergence analysis of the method developed in Section 2. For this purpose, we first let 𝐮=(𝑢_𝑖+1/2), 𝐰=(𝑤_𝑖+1/2), 𝐜=(𝑐_𝑖), 𝐭=(𝑡_𝑖), and 𝐞=(𝑒_𝑖+1/2) be 𝑛-dimensional column vectors. Here 𝑒_𝑖+1/2=𝑢_𝑖+1/2−𝑤_𝑖+1/2 is the discretization error. Thus, we can write our method as follow: 𝐀𝐮=𝐜+𝐭, (3.1a)𝐀𝐰=𝐜, (3.1b)𝐀𝐞=𝐭, (3.1c)

where 𝐀=𝐀₀+1ℎ1920³𝐁𝐏, (3.2)𝐏=diag(𝑝_𝑖−1/2),𝑖=1,2,…,𝑛, with 𝑝_𝑖−1/2≠0 for 𝑛/4<𝑖≤3𝑛/4, 𝐀₀=9−10⋯⋯⋯0−23−10⋯⋯01−33−10⋯00⋱⋱⋱⋱⋯0⋮⋱⋱⋱⋱⋱⋮0⋯01−33−10⋯⋯01−32, (3.3) and the lower triangular matrix𝐁=612360⋯⋯⋯01062−310⋯⋯003200⋯00⋱…⋱⋱⋱0⋮⋱⋱⋱⋱⋱⋮0⋯003200⋯⋯0809. (3.4) For the vector 𝐜, we have𝑐_𝑖=8𝛼+3ℎ𝛽₁−1ℎ160³𝐹₁,𝑖=1,ℎ𝛽₁−1ℎ1920³𝐹₂−1,𝑖=2,6ℎ³𝐹_𝑖𝑛,3≤𝑖≤4−1,3𝑛4−1+3≤𝑖≤𝑛−1,6ℎ³𝐹_𝑖𝑛+𝑟,𝑖=4𝑛,𝑖=4−1+1,6ℎ³𝐹_𝑖𝑛+7𝑟,𝑖=4−1+2,6ℎ³𝐹_𝑖,𝑛+6𝑟4+3≤𝑖≤3𝑛4−1−1,6ℎ³𝐹_𝑖+5𝑟,𝑖=3𝑛4,𝑖=3𝑛4−1+1,ℎ12³𝐹_𝑖−𝑟,𝑖=3𝑛4+2,ℎ𝛽₂−1ℎ1920³𝐹_𝑛,𝑖=𝑛, (3.5) where𝐹_𝑖=51𝑓_1/2+3𝑓_3/2,𝑖=1,809𝑓_1/2+1062𝑓_3/2−3125𝑓_5/2,𝑖=2,−𝑓_𝑖−5/2+8𝑓_𝑖−3/2+𝑓_𝑖+1/2,3≤𝑖≤𝑛−1,−31𝑓_𝑛−5/2+1062𝑓_𝑛−3/2+908𝑓_𝑛−1/2,𝑖=𝑛. (3.6)

Our main purpose now is to derive a bound on ‖𝐞‖, where ‖⋅‖ represents the ∞-norm. It has been shown in [1] that 𝐀₀⁻¹ exits and satisfies the equation𝐀₀⁻¹=4𝑛³−𝑛+3,48 (3.7) which upon using (2.1) we obtain𝐀₀⁻¹=3ℎ³−(𝑏−𝑎)ℎ²+4(𝑏−𝑎)³48ℎ³. (3.8) Thus, using (3.1a)–(3.1c) and (3.2) and the fact that ‖𝐁‖=2560 and ‖𝐏‖≤|𝑝(𝑥)|, we get‖𝐞‖≤1523𝜆𝑀₆ℎ³𝑝(𝑥ℎ25603−4𝜆)≅𝑂³, (3.9) where 𝜆=(1/48)[ℎ³−(𝑏−𝑎)ℎ²+(𝑏−𝑎)³] and 𝑀₆=max|𝑦^(𝑣𝑖)(𝑥)|, see [1] for more details. Relation (3.9) indicates that (3.1b) is a third-order convergent method.

Remark 3.1. The matrix 𝐴₀ and its inverse were first introduced in [1]. The derivation of (3.7) (the elements of 𝐴₀⁻¹) given in [1] was derived using the definition of the matrix inverse. The derivation involves a long and complicated algebraic manipulations. The interested reader may try to derive it.

4. Applications and Computational Results

To illustrate the application of the numerical method developed in the previous sections, we consider the third-order obstacle boundary value problem of finding 𝑢 such that−𝑢^′′′[],[],≥𝑓,onΩ=0,1𝑢≥𝜓,onΩ=0,1−𝑢^′′′[],𝑢(0)−𝑓𝑢−𝜓=0,onΩ=0,1=0,𝑢^′(0)=0,𝑢^′(1)=0, (4.1) where 𝑓(𝑥) is a continuous function and 𝜓(𝑥) is the obstacle function. We study the problem (4.1) in the framework of variational inequality approach. To do so, we first define the set 𝐾 as𝐾=𝑣∶𝑣∈𝐻²₀(Ω)∶𝑣≥𝜓onΩ, (4.2) which is a closed convex set in 𝐻²₀(Ω), where 𝐻²₀(Ω) is a Sobolev (Hilbert) space, see [8]. One can easily show that the energy functional associated with the problem (4.1) is𝐼[𝑣]=−¹₀𝑑³𝑣𝑑𝑥³𝑑𝑣𝑑𝑥𝑑𝑥−2¹₀𝑓(𝑥)𝑑𝑣𝑑𝑥𝑑𝑥,∀𝑑𝑣=𝑑𝑥∈𝐾₁₀𝑑²𝑣𝑑𝑥²²𝑑𝑥−2¹₀𝑓(𝑥)𝑑𝑣=(𝑣(𝑣,(𝑣=𝑑𝑥𝑑𝑥𝑇𝑣,𝑔)−2𝑓,𝑔)𝑇𝑢,𝑔)¹₀𝑑²𝑢𝑑𝑥²𝑑²𝑣𝑑𝑥²(𝑣=𝑑𝑥,𝑓,𝑔)¹₀𝑓(𝑥)𝑑𝑣𝑑𝑥𝑑𝑥, (4.3)(4.4) and 𝑔=𝑑/𝑑𝑥 is the linear operator.

It is clear that the operator 𝑇 defined by (4.4) is linear, 𝑔-symmetric, and 𝑔-positive. Using the technique of Noor [13, 15], one can easily show that the minimum 𝑢∈𝐻 of the functional 𝐼[𝑣] defined by (4.3) associated with the problem (4.1) on the closed convex set 𝐾 can be characterized by the inequality of the type𝑇𝑢,𝑔(𝑣)−𝑔(𝑢)≥𝑓,𝑔(𝑣)−𝑔(𝑢),∀𝑔(𝑣)∈𝐾, (4.5) which is exactly the general variational inequality, considered by Noor [10] in 1988. It is worth mentioning that a wide class of unrelated odd-order and nonsymmetric equilibrium problems arising in regional, physical, mathematical, engineering, and applied sciences can be studied in the unified and general framework of the general variational inequalities, see [1–20].

Using the penalty function technique of Lewy and Stampacchia [9], we can characterize the problem (4.1) as−𝑢^′′′𝑢(0)+𝜈𝑢−𝜓𝑢−𝜓=𝑓,0<𝑥<1,=𝑢^′(0)=𝑢^′(1)=0, (4.6) where𝜈{𝑡}=1,for𝑡≥00,for𝑡<0 (4.7) is a discontinuous function and is known as the penalty function and 𝜓 is the given obstacle function defined by𝜓(𝑥)=1−1,for0≤𝑥≤4,341≤𝑥≤1,1,for43≤𝑥≤4. (4.8) From equations (4.3)–(4.8), we obtain the following system of differential equations:𝑢^′′′=1𝑓,for0≤𝑥≤4,341≤𝑥≤1,𝑢+𝑓−1,for43≤𝑥≤4, (4.9) with the boundary conditions𝑢(0)=𝑢^′(0)=𝑢^′(1)=0 (4.10) and the condition of continuity of 𝑢,𝑢^′ and 𝑢^′′ at 𝑥=1/4 and 3/4. Note that the system of differential equations (4.7) is a special form of the system (1.1a) with 𝑝(𝑥)=1 and 𝑟=−1.

Example 4.1. For 𝑓=0, the system of differential equations (4.9) reduces to 𝑢^′′′=10,for0≤𝑥≤4,341≤𝑥≤1,𝑢−1,for43≤𝑥≤4, (4.11) with the boundary conditions (4.10). The analytical solution for this problem is 𝑢(𝑥)=12𝑎₁𝑥²1,0≤𝑥≤4,1+𝑎₂𝑒^𝑥+𝑒^−𝑥/2𝑎₃√cos32𝑥+𝑎₄√sin32𝑥,143≤𝑥≤4,12𝑎₅𝑥(𝑥−2)+𝑎₆,34≤𝑥≤1. (4.12) We can find the constants 𝑎_𝑖,𝑖=1,2,…,6, by solving a system of linear equations constructed by applying the continuity conditions of 𝑢,𝑢^′, and 𝑢^′′ at 𝑥=1/4 and 3/4, see [3] for more details.

For different values of ℎ, the boundary value problem defined by (4.10) and (4.11) was solved using the numerical methods developed in the previous sections and the observed maximum errors ‖𝐞‖ are listed in Table 1. Also, we give in Table 1 the numerical results for the finite difference methods at midkonts introduced in [1, 3, 5, 16]. It is clear from this table that our present method produced better results than the other ones. However, the numerical results may indicate that we have second-order approximations. This is due to the fact that the third derivative is not continuous across the interfaces.

Table 1 

Observed maximum errors ‖𝑒‖.

Now, let 𝑤_𝑖−1/2 be an approximate value of 𝑢_𝑖−1/2, for 𝑖=1,2,…,𝑛, computed by the numerical method developed in the previous sections. Then having the values of 𝑤_𝑖−1/2, for 𝑖=1,2,…,𝑛, we can compute 𝑤_𝑖≈𝑢_𝑖 using the second-order interpolating𝑢_𝑖=12𝑢_𝑖+1/2+𝑢_𝑖−1/2+𝑂ℎ², (4.13) for 𝑖=1,2,…,𝑛. Note that we make use of the boundary condition 𝑢^′(1)=0 to compute the value of 𝑤_𝑛. The computations of 𝑤_𝑖, 𝑖=1,2,…,𝑛, give us the opportunity for a fare comparisons with the other methods discussed in [2, 7, 17, 19], which approximate the solution of problems (4.11) at the nodes. In Table 2, we list the maximum value for the errors max_𝑖|𝑢_𝑖−𝑤_𝑖| for different ℎ values for our present method and the others. From Table 2, it can be noted that our present method gives the best results. We mentioned here in passing that the numerical results for the method developed in [4] are worse than those given in Table 2 and are not presented here.

Table 2 

Observed maximum errors.

5. Conclusion

As mentioned in [1, 3, 5, 16] where two-stage first- and second-order method was developed, we have noticed from our experiments that the maximum value of the error occurs near the center of the interval and not around 𝑥=1/4 nor 3/4, where the solution satisfies the extraconditions. On the other hand, it was noticed from the experiments done by the authors in [2, 4, 17] that the maximum error occurs at near the indicated values of 𝑥. Thus, we can conclude that the extra conditions at 𝑥=1/4 and 3/4 have little effect on the performance of the methods that first approximate the solution at the midknots, whereas for the other two methods these added conditions introduce a serious drawback in the accuracy.

Acknowledgments

This paper is supported by the Visiting Professor Program of King Saud University, Riyadh, Saudi Arabia. The authors would like to thank the referees for their valuable input and comments.