# Finite Difference Method for Solving a System of

Authors: Muhammad Aslam Noor; Eisa Al-Said; Khalida Inayat Noor
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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)u′′′={f(x),a≤x≤c,p(x)u(x)+f(x)+r,c≤x≤d,f(x),d≤x≤b,
with the boundary conditions(1.1b)u(a)=α,u′(a)=β1,u′(b)=β2,and the continuity conditions of u, u′, and u′′ at c and d. Here, f and p are continuous functions on [a,b] and [c,d], respectively. The parameters r, α, β1, and β2 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 f(x) and p(x). 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 takec=(3a+b)/4 and d=(a+3b)/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 [a,b] into n equal subintervals using the grid points xi=a+ih,i=0,1,2,…,n, x0=a, xn=b and(2.1)h=b-an,
where n is a positive integer chosen such that both n/4 and 3n/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:(2.2)9u1/2-u3/2=8u0+3hu0′-1160h3[6u0′′′+51u1/2′′′+3u3/2′′′]+t1,fori=1,-2u1/2+3u3/2-u5/2=+hu0′-11920h3[809u1/2′′′+1062u3/2′′′-31u5/2′′′]+t2,fori=2,ui-5/2-3ui-3/2+3ui-1/2-ui+1/2=16h3[-ui-5/2′′′+6ui-3/2′′′+ui+1/2′′′]+ti,for3≤i≤n-1,un-5/2-3un-3/2+2un-1/2=hun′-11920h3[-31un-5/2′′′+1062un-3/2′′′+809un-1/2′′′]+tn,fori=n,
where(2.3)ui+1/2′′′={fi+1/2,for0≤i≤n4-1,3n4≤i≤n-1pi+1/2ui+1/2+fi+1/2+r,forn4≤i≤3n4-1,fi+1/2=f(xi+1/2),i=0,1,2,…,n-1, the relation (2.2) forms a system of n linear equations in the unknowns ui-1/2,i=1,2,…,n. The local truncation errors associate with (2.2) are given by(2.4)ti={49320h6u0(vi)+O(h7),fori=1,-533840h6u0(vi)+O(h7),fori=2,15237680h6ui(vi)+O(h7),for3≤i≤n-1,-533840h6ui(vi)+O(h7),fori=n,
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 Section2. For this purpose, we first let u=(ui+1/2), w=(wi+1/2), c=(ci), t=(ti), and e=(ei+1/2) be n-dimensional column vectors. Here ei+1/2=ui+1/2-wi+1/2 is the discretization error. Thus, we can write our method as follow: (3.1a)Au=c+t,(3.1b)Aw=c,(3.1c)Ae=t,where(3.2)A=A0+11920h3BP,P=diag⁡(pi-1/2),i=1,2,…,n, with pi-1/2≠0 for n/4<i≤3n/4, (3.3)A0=[9-10⋯⋯⋯0-23-10⋯⋯01-33-10⋯00⋱⋱⋱⋱⋯0⋮⋱⋱⋱⋱⋱⋮0⋯01-33-10⋯⋯01-32],
and the lower triangular matrix(3.4)B=[612360⋯⋯⋯01062-310⋯⋯003200⋯00⋱…⋱⋱⋱0⋮⋱⋱⋱⋱⋱⋮0⋯003200⋯⋯0809].
For the vector c, we have(3.5)ci={8α+3hβ1-1160h3F1,i=1,hβ1-11920h3F2,i=2,-16h3Fi,3≤i≤n4-1,3n4+3≤i≤n-1,-16h3[Fi+r],i=n4,i=n4+1,-16h3[Fi+7r],i=n4+2,-16h3[Fi+6r],n4+3≤i≤3n4-1,-16h3[Fi+5r],i=3n4,i=3n4+1,-112h3[Fi-r],i=3n4+2,hβ2-11920h3Fn,i=n,
where(3.6)Fi={51f1/2+3f3/2,i=1,809f1/2+1062f3/2-3125f5/2,i=2,-fi-5/2+8fi-3/2+fi+1/2,3≤i≤n-1,-31fn-5/2+1062fn-3/2+908fn-1/2,i=n.Our main purpose now is to derive a bound on∥e∥, where ∥·∥ represents the ∞-norm. It has been shown in [1] that A0-1 exits and satisfies the equation(3.7)‖A0-1‖=4n3-n+348,
which upon using (2.1) we obtain(3.8)‖A0-1‖=3h3-(b-a)h2+4(b-a)348h3.
Thus, using (3.1a)–(3.1c) and (3.2) and the fact that ∥B∥=2560 and ∥P∥≤|p(x)|, we get(3.9)‖e‖≤1523λM6h32560[3-4λ|p(x)|]≅O(h3),
where λ=(1/48)[h3-(b-a)h2+(b-a)3] and M6=max⁡|y(vi)(x)|, see [1] for more details. Relation (3.9) indicates that (3.1b) is a third-order convergent method.Remark 3.1.
The matrixA0 and its inverse were first introduced in [1]. The derivation of (3.7) (the elements of A0-1) 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 findingu such that(4.1)-u′′′≥f,on⁡Ω=[0,1],u≥ψ,on⁡Ω=[0,1],[-u′′′-f][u-ψ]=0,on⁡Ω=[0,1],u(0)=0,u′(0)=0,u′(1)=0,
where f(x) is a continuous function and ψ(x) 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 K as(4.2)K={v:v∈H02(Ω):v≥ψon⁡Ω},
which is a closed convex set in H02(Ω), where H02(Ω) is a Sobolev (Hilbert) space, see [8]. One can easily show that the energy functional associated with the problem (4.1) is(4.3)I[v]=-∫01(d3vdx3)(dvdx)dx-2∫01f(x)(dvdx)dx,∀dvdx∈K=∫01(d2vdx2)2dx-2∫01f(x)(dvdx)dx=〈Tv,g(v)〉-2〈f,g(v)〉,(4.4)〈Tu,g(v)〉=∫01(d2udx2)(d2vdx2)dx,〈f,g(v)〉=∫01f(x)dvdxdx,
and g=d/dx is the linear operator.It is clear that the operatorT defined by (4.4) is linear, g-symmetric, and g-positive. Using the technique of Noor [13, 15], one can easily show that the minimum u∈H of the functional I[v] defined by (4.3) associated with the problem (4.1) on the closed convex set K can be characterized by the inequality of the type(4.5)〈Tu,g(v)-g(u)〉≥〈f,g(v)-g(u)〉,∀g(v)∈K,
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(4.6)-u′′′+ν{(u-ψ)}(u-ψ)=f,0<x<1,u(0)=u′(0)=u′(1)=0,
where(4.7)ν{t}={1,fort≥00,fort<0
is a discontinuous function and is known as the penalty function and ψ is the given obstacle function defined by(4.8)ψ(x)={-1,for0≤x≤14,34≤x≤1,1,for14≤x≤34.
From equations (4.3)–(4.8), we obtain the following system of differential equations:(4.9)u′′′={f,for0≤x≤14,34≤x≤1,u+f-1,for14≤x≤34,
with the boundary conditions(4.10)u(0)=u′(0)=u′(1)=0
and the condition of continuity of u,u′ and u′′ at x=1/4 and 3/4. Note that the system of differential equations (4.7) is a special form of the system (1.1a) with p(x)=1 and r=-1.Example 4.1.
Forf=0, the system of differential equations (4.9) reduces to
(4.11)u′′′={0,for0≤x≤14,34≤x≤1,u-1,for14≤x≤34,
with the boundary conditions (4.10). The analytical solution for this problem is
(4.12)u(x)={12a1x2,0≤x≤14,1+a2ex+e-x/2[a3cos⁡32x+a4sin⁡32x],14≤x≤34,12a5x(x-2)+a6,34≤x≤1.
We can find the constants ai,i=1,2,…,6, by solving a system of linear equations constructed by applying the continuity conditions of u,u′, and u′′ at x=1/4 and 3/4, see [3] for more details.For different values ofh, 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 ∥e∥ 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∥e∥.
hNew method[1][3][5][16]1/204.05×10-54.96×10-51.48×10-46.74×10-41.04×10-31/409.78×10-61.24×10-53.70×10-53.04×10-42.60×10-41/802.31×10-63.10×10-69.24×10-61.37×10-46.49×10-5Now, letwi-1/2 be an approximate value of ui-1/2, for i=1,2,…,n, computed by the numerical method developed in the previous sections. Then having the values of wi-1/2, for i=1,2,…,n, we can compute wi≈ui using the second-order interpolating(4.13)ui=12[ui+1/2+ui-1/2]+O(h2),
for i=1,2,…,n. Note that we make use of the boundary condition u′(1)=0 to compute the value of wn. The computations of wi, i=1,2,…,n, 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⁡i|ui-wi| for different h 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.
hNew method[2][7][17][19]1/322.57×10-55.53×10-45.30×10-45.32×10-44.05×10-41/646.38×10-62.61×10-42.52×10-42.56×10-42.24×10-41/1281.47×10-61.27×10-41.23×10-41.26×10-41.15×10-4


## 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 x=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 x. Thus, we can conclude that the extra conditions at x=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.



---

*Source: 351764-2012-04-17.xml*