Solution of Boundary Value Obstacle Problems Using MQ-RBF and IMQ-RBF

Abstract

A kind of numerical method which is based on multiquadric RBF, inverse multiquadric RBF, and Wu-Schaback operators is presented for solving second-order and third-order boundary value problems associated with obstacle, unilateral, and contact problems. The algorithms are proved to be highly accurate and easy to implement. Some numerical tests are also presented to show the efficiency of the algorithms.

1. Introduction

We consider the numerical solutions of the following second-order boundary value problems: 𝑢^′′(𝑥)=𝑓(𝑥)𝑔(𝑥)𝑢(𝑥)(𝑥)𝑓(𝑥),𝑎≤𝑥≤𝑐,+𝑓+𝑟,𝑐≤𝑥≤𝑑,,𝑑≤𝑥≤𝑏, (1) with the boundary conditions 𝑢(𝑎)=𝛼,𝑢(𝑏)=𝛽 (2) and the continuity conditions of 𝑢 and 𝑢^′ at 𝑐 and 𝑑. We also consider the following third-order boundary value problems: 𝑢^′′′(𝑥)=𝑓(𝑥)𝑔(𝑥)𝑢(𝑥)(𝑥)𝑓(𝑥),𝑎≤𝑥≤𝑐,+𝑓+𝑟,𝑐≤𝑥≤𝑑,,𝑑≤𝑥≤𝑏, (3) with the boundary conditions 𝑢(𝑎)=𝛼₁,𝑢^′(𝑎)=𝛼₂,𝑢^′(𝑏)=𝛽 (4) and the continuity conditions of 𝑢, 𝑢^′, and 𝑢^′′ at 𝑐 and 𝑑. In the two problems, 𝑓(𝑥) and 𝑔(𝑥) are continuous functions and 𝛼, 𝛽, 𝛼₁, 𝛼₂ are real constants. The two problems arise usually in obstacle, contact, unilateral, and equilibrium problems associated with economics, transportation, oceanography, fluid flow through porous media, and many other fields of pure and applied sciences; see [1–12]. In the literature, some techniques were used to solve the previous two systems of second-order and third-order boundary value problems associated with obstacles. For example, for problem (1), Noor and Khalifa [13] used first-order accuracy collocation method with cubic splines as basis functions. Al-Said et al. [14] developed second- and fourth-order finite difference and spline methods. Al-Said [5] developed and analyzed quadratic and cubic splines methods. He proved that both quadratic and cubic spline methods can be used to produce second-order smooth approximations for the solution of (1) and its first derivative over the whole interval [𝑎,𝑏]. Al-Said [7] also used cubic spline polynomial functions for solving problem (1).

For the problem (3), Khalifa and Noor [15] applied quintic spline method for contact problems. Noor and Khalifa [16] applied quartic spline method for odd-order obstacle problems. Gao and Chi [17] applied quartic B-spline method for third-order obstacle problems and Siraj-ul-Islam et al. [12] proposed nonpolynomial spline methods. Some of these methods require huge computational work for finding the numerical solutions.

In this paper, we use multiquadric RBF (MQs) and inverse multiquadric RBF (IMQs) to solve systems of second-order and third-order boundary value problems associated with obstacle. Recently, radial base functions (RBF) have been widely applied to the numerical solution of many types of PDE; see [18, 19]. The advantage of techniques of using RBF is that mesh grid can be avoided and high accuracy can be achieved. However, to the authors’ knowledge, little work has been done involving the usage of RBF in the solution of boundary value obstacle problems. Among Radial Basis Functions, MQs and IMQs are especially suitable for solving second- and third-order boundary value problems associated with obstacles because we can integrate IMQs to get MQs. Our technique is that we first use IMQs or MQs quasi-interpolation to approximate the derivatives of the solution of problems of (1) and (3); then we get their integrals to approximate the solutions. The following part of this paper is organized as follows. In Section 2, we review some basics of RBFs. In Section 3, we present the numerical schemes for problem (1) and (3). In Section 4, we do some numerical experiments. Section 5 is our conclusion.

2. RBF and MQ Quasi-Interpolation

RBF was proposed in 1951 [20] to deal with geological problem.

For a given region Ω∈𝑅^𝑚 and a set distinct interpolation points, 𝑋=𝑥₁,𝑥₂,…𝑥_𝑁⊂Ω. (5) Given a function 𝑓∶𝑋→𝑅, construct the interpolation to 𝑓𝑆_𝑓,𝑋(𝑥)=^𝑁_𝑗=1𝛼_𝑗𝜑𝑥−𝑥_𝑗,for𝑥∈Ω, (6) where ‖⋅‖ denote Euclidean norm and 𝜑 is a certain RBF. The undetermined coefficients {𝛼₁,𝛼₂,…,𝛼_𝑁} can be obtained by the linear system ^𝑁_𝑗=1𝛼_𝑗𝜑𝑥_𝑖−𝑥_𝑗=𝑓𝑥_𝑖,𝑖=1,2,…,𝑁. (7) We adopt multiquadric (MQs) RBF 𝜙(𝑟)=√𝑠²+𝑟² (8) and inverse multiquadric (IMQs) RBF 𝜑(𝑟)=𝑠²𝑠²+𝑟^23/2 (9) to solve second-order and third-order boundary value obstacle problems, where 𝑠∈𝑅⁺ is a shape parameter. From [21], we know that the linear system (7) has unique solution. It is easy to see that √𝑠²+𝑟^2′′=𝑠²𝑠²+𝑟^23/2. (10)

We also adopt a kind of Quasi-interpolation to numerically solve problem (1) and (3). The Quasi-interpolation of a univariate function 𝑓∶[𝑎,𝑏]→𝑅 with MQs is constructed as follows. On the scattered points 𝑎=𝑥₀<𝑥₁<⋯<𝑥_𝑛=𝑏,ℎ∶=max_{1≤𝑖≤𝑛}𝑥_𝑖−𝑥_𝑖−1. (11) Quasi-interpolant of 𝑓 usually is as follows: 𝑀𝑓(𝑥)=^𝑛_𝑖=0𝑓𝑥_𝑖𝜓_𝑖(𝑥), (12) where 𝜓_𝑖(𝑥) is linear combinations of the MQs. Consider 𝜙_𝑖𝑥−𝑥_𝑖=𝑠²+𝑥−𝑥_𝑖²,𝑠∈𝑅⁺, (13) where 𝑠 is a shape parameter.

Given data {(𝑥_𝑖,𝑓_𝑖)}^𝑛_𝑖=1, Wu–Schabacks 𝐿_𝐷 operator is 𝐿_𝐷𝑓(𝑥)=^𝑛_𝑗=0𝑓𝑥_𝑗𝜓_𝑗(𝑥), (14) where 𝜓₀(𝑥)=12+𝜙₁𝑥−𝑥₁−𝑥−𝑥₀2𝑥₁−𝑥₀,𝜓₁(𝑥)=𝜙₂𝑥−𝑥₂−𝜙₁𝑥−𝑥₁2𝑥₂−𝑥₁−𝜙₁𝑥−𝑥₁−𝑥−𝑥₀2𝑥₁−𝑥₀,𝜓_𝑛−1(𝑥)=𝑥_𝑛−𝑥−𝜙_𝑛−1𝑥−𝑥_𝑛−12𝑥_𝑛−𝑥_𝑛−1−𝜙_𝑛−1𝑥−𝑥_𝑛−1−𝜙_𝑛−2𝑥−𝑥_𝑛−22𝑥_𝑛−1−𝑥_𝑛−2,𝜓_𝑛(𝑥)=12+𝜙_𝑛−1𝑥−𝑥_𝑛−1−𝑥_𝑛−𝑥2𝑥_𝑛−𝑥_𝑛−1,𝜓_𝑗(𝑥)=𝜙_𝑗−1𝑥−𝑥_𝑗−1−𝜙_𝑗𝑥−𝑥_𝑗2𝑥_𝑗−1−𝑥_𝑗−𝜙_𝑗𝑥−𝑥_𝑗−𝜙_𝑗−1𝑥−𝑥_𝑗−12𝑥_𝑗−𝑥_𝑗−1,𝑗=2,…,𝑛−2. (15) The advantage of 𝐿_𝐷 operator is that high accuracy can be achieved without solving any linear system.

For ∫𝜓_𝑖(𝑥)𝑑𝑥, we need only to integrate 𝜙_𝑖(|𝑥−𝑥_𝑖|), 𝜙_𝑖𝑥−𝑥_𝑖=𝑑𝑥𝑠²+𝑥−𝑥_𝑗²=𝑑𝑥𝑥−𝑥_𝑗𝑥−𝑥_𝑗²+𝑠²2+𝑠²2ln𝑥−𝑥_𝑗+𝑥−𝑥_𝑗²+𝑠². (16) For ∫𝑑𝑥∫𝜓_𝑖(𝑥)𝑑𝑥, which is denoted with Φ_𝑖(𝑥), we have 𝑑𝑥𝜙_𝑖𝑥−𝑥_𝑖=𝑑𝑥𝑑𝑥𝑠²+𝑥−𝑥_𝑗²𝑑𝑥=𝑥−𝑥_𝑗²+𝑠^21.56+𝑠²2𝑥−𝑥_𝑗ln𝑥−𝑥_𝑗+𝑥−𝑥_𝑗²+𝑠²−𝑠²2𝑥−𝑥_𝑗²+𝑠². (17)

3. Numerical Algorithms

To solve problem (1) and (3), we set 𝑎=𝑥₁<𝑥₂<⋯<𝑥_𝑚=𝑐=𝑥_𝑚+1<𝑥_𝑚+2<⋯<𝑥_𝑚+𝑛=𝑑=𝑥_𝑚+𝑛+1<𝑥_𝑚+𝑛+2<⋯<𝑥_{𝑚+𝑛+𝑘}=𝑏. (18) Construct (𝐿_𝐷𝑓)(𝑥)=∑^𝑚_𝑗=1𝑓(𝑥_𝑗)𝜓_𝑗(𝑥) on [𝑎,𝑐]  to approximate  𝑓(𝑥), 𝑥∈[𝑎,𝑐]; and (𝐿_𝐷𝑓)(𝑥)=∑^𝑘_𝑗=1𝑓(𝑥_{𝑚+𝑛+𝑗})𝜓_{𝑚+𝑛+𝑗}(𝑥) on [𝑑,𝑏]  to approximate  𝑓(𝑥), 𝑥∈[𝑑,𝑏].

On the interval [𝑐,𝑑], let 𝑢^′′(𝑥)=∑^𝑛_𝑗=1𝛼_𝑗𝜑_𝑚+𝑗(|𝑥−𝑥_𝑚+𝑗|), where 𝜑_𝑚+𝑗𝑥−𝑥_𝑚+𝑗=𝑠²𝑠²+𝑥−𝑥_𝑚+𝑗^23/2, (19) and 𝑠 is a shape parameter and 𝛼_𝑗, 𝑗=1,…,𝑛 are coefficients to be identified. We can easily verify that 𝜑_𝑗𝑥−𝑥_𝑗𝑑𝑥=𝑥−𝑥_𝑗𝑠²+𝑥−𝑥_𝑗²+𝜆,𝑥−𝑥_𝑗𝑠²+𝑥−𝑥_𝑗²+𝜆𝑑𝑥=𝜙_𝑗𝑥−𝑥_𝑗+𝜆𝑥+𝜇, (20)(21) where 𝜆, 𝜇 are constants. Therefore, 𝑢(𝑥)=𝑑𝑥^𝑛_𝑗=1𝛼_𝑗𝜑_𝑚+𝑗𝑥−𝑥_𝑚+𝑗=𝑑𝑥_𝑛𝑗=1𝛼_𝑗𝜙_𝑚+𝑗𝑥−𝑥_𝑚+𝑗[],+𝜆𝑥+𝜇,𝑥∈𝑐,𝑑 (22) where 𝜆, 𝜇 are constants and 𝜙_𝑗(|𝑥−𝑥_𝑗|) are defined in (15). Therefore, we have the approximation to the solution to problem (1) as follows: 𝑢=𝑑𝑥^𝑚_𝑗=1𝑓𝑥_𝑗𝜓_𝑗(𝑥)=𝑑𝑥_𝑚𝑗=1𝑓𝑥_𝑗Φ_𝑗(𝑥)+𝑐₁𝑥+𝑐₂,𝑥∈[𝑎,𝑐_𝑛𝑗=1],𝛼_𝑗𝜙_𝑚+𝑗𝑥−𝑥_𝑚+𝑗+𝜆𝑥+𝜇,𝑥∈[𝑐,𝑑],𝑑𝑥^𝑘_𝑗=1𝑓𝑥_{𝑚+𝑛+𝑗}𝜓_{𝑚+𝑛+𝑗}(𝑥)=𝑑𝑥_𝑘𝑗=1𝑓𝑥_{𝑚+𝑛+𝑗}Φ_{𝑚+𝑛+𝑗}(𝑥)+𝑐₃𝑥+𝑐₄,𝑥∈[𝑑,𝑏], (23) where Φ_𝑗(𝑥)=∫𝑑𝑥∫𝜓_𝑖(𝑥)𝑑𝑥 and is computed in (17), and 𝑐₁, 𝑐₂, 𝜆, 𝜇, 𝑐₃, 𝑐₄ are constants. We need to determine the coefficients 𝑐₁, 𝑐₂, 𝛼₁,…,𝛼_𝑛, 𝜆, 𝜇, 𝑐₃, 𝑐₄ to identify the numerical solution (23) to problem (1).

Let 𝐻₁(𝑥)=^𝑚_𝑗=1𝑓𝑥_𝑗Φ_𝑗(𝑥)+𝑐₁𝑥+𝑐₂[],𝐻,𝑥∈𝑎,𝑐₂(𝑥)=^𝑛_𝑗=1𝛼_𝑗𝜙_𝑚+𝑗𝑥−𝑥_𝑚+𝑗,[],𝐻+𝜆𝑥+𝜇𝑥∈𝑐,𝑑₃(𝑥)=^𝑘_𝑗=1𝑓𝑥_{𝑚+𝑛+𝑗}Φ_{𝑚+𝑛+𝑗}(𝑥)+𝑐₃𝑥+𝑐₄,𝑥∈[𝑑,𝑏]. (24) Substitute (24) into problem (1) at points (18) to get the following linear system: 𝐻₁(𝑎)=𝛼,𝐻₁(𝑐)=𝐻₂(𝑐),𝐻^′₁(𝑐)=𝐻^′₂(𝑐),_𝑛𝑗=1𝛼_𝑗𝜑_𝑚+𝑗𝑥_𝑚+𝑖−𝑥_𝑚+𝑗=𝑔𝑥_𝑚+𝑖𝐻₂𝑥_𝑚+𝑖+𝑓𝑥_𝑚+𝑖𝐻+𝑟,𝑖=1,…,𝑛,₂(𝑑)=𝐻₃(𝑑),𝐻^′₂(𝑑)=𝐻^′₃(𝑑),𝐻₃(𝑏)=𝛽, (25) where 𝐻^′₁(𝑥), 𝐻^′₃(𝑥) can be computed in (16) and 𝐻^′₂(𝑥) can be computed in (20).

Solving the linear system (24), we get 𝑐₁, 𝑐₂, 𝛼₁,…,𝛼_𝑛, 𝜆, 𝜇, 𝑐₃, 𝑐₄ to identify the numerical solution (23) to problem (1).

To solve problem (3), we use ^𝑚_𝑗=1𝛼_𝑗𝜑_𝑗𝑥−𝑥_𝑗,_𝑛𝑗=1𝛼_𝑚+𝑗𝜑_𝑚+𝑗𝑥−𝑥_𝑚+𝑗,_𝑘𝑗=1𝛼_{𝑚+𝑛+𝑗}𝜑_{𝑚+𝑛+𝑗}𝑥−𝑥_{𝑚+𝑛+𝑗} (26) to approximate 𝑢^′′′ on intervals [𝑎,𝑐], [𝑐,𝑑], and [𝑑,𝑏], respectively, where 𝜑_𝑚+𝑗𝑥−𝑥_𝑚+𝑗=𝑠²𝑠²+𝑥−𝑥_𝑚+𝑗^23/2 (27) and 𝑠 is a shape parameter. To get 𝑢^′′, 𝑢^′, and 𝑢, we need only to compute 𝜑_𝑗𝑥−𝑥_𝑗𝑑𝑥=𝑥−𝑥_𝑗𝑠²+𝑥−𝑥_𝑗^21/2𝜑+𝜆,𝑑𝑥_𝑗𝑥−𝑥_𝑗𝑑𝑥=𝜙_𝑗𝑥−𝑥_𝑗+𝜆𝑥+𝜇, (28)(29) where 𝜆,  𝜇 are constants and 𝜙_𝑗 is defined in (13). And we also have 𝜙_𝑗𝑥−𝑥_𝑗+𝜆𝑥+𝜇𝑑𝑥=𝜁_𝑗(𝑥)+𝜆²2𝑥²+𝜇+𝜈, (30) where 𝜈 is a constant and 𝜁_𝑗(𝑥)=𝑥−𝑥_𝑗𝑥−𝑥_𝑗²+𝑠²2+𝑠²2ln𝑥−𝑥_𝑗+𝑥−𝑥_𝑗²+𝑠². (31) Therefore, we need to identify the following numerical solution to (3): 𝑢=^𝑚_𝑗∑₌₁𝛼_𝑗𝜁_𝑗(𝑥)+𝑐₁𝑥²+𝑐₂𝑥+𝑐₃,_𝑛𝑗∑₌₁𝛼_𝑚+𝑗𝜁_𝑚+𝑗(𝑥)+𝜆𝑥²+𝜇𝑥+𝜈,_𝑘𝑗∑₌₁𝛼_{𝑚+𝑛+𝑗}𝜁_{𝑚+𝑛+𝑗}(𝑥)+𝑐₄𝑥²+𝑐₅𝑥+𝑐₆. (32) To determine the coefficients 𝑐₁, 𝑐₂, 𝑐₃, 𝛼₁,…,𝛼_{𝑚+𝑛+𝑘}, 𝜆, 𝜇, 𝜈, 𝑐₃, 𝑐₄, 𝑐₅, denote 𝐻₁(𝑥)=^𝑚_𝑗=1𝛼_𝑗𝜁_𝑗(𝑥)+𝑐₁𝑥²+𝑐₂𝑥+𝑐₃,𝐻₂(𝑥)=^𝑛_𝑗=1𝛼_𝑚+𝑗𝜁_𝑚+𝑗(𝑥)+𝜆𝑥²𝐻+𝜇𝑥+𝜈,₃(𝑥)=^𝑘_𝑗=1𝛼_{𝑚+𝑛+𝑗}𝜁_{𝑚+𝑛+𝑗}(𝑥)+𝑐₄𝑥²+𝑐₅𝑥+𝑐₆, (33) and 𝐻^′₁(𝑥), 𝐻^′₂(𝑥), 𝐻^′₃(𝑥) can be computed in (29) and 𝐻^′′₁(𝑥), 𝐻^′′₂(𝑥), 𝐻^′′₃(𝑥) can be computed in (28). Substitute (25) into problem (3) at points (18) to get the following linear system: 𝐻₁(𝑎)=𝛼₁,𝐻^′₁(𝑎)=𝛼₂,𝐻₁(𝑐)=𝐻₂(𝑐),𝐻^′₁(𝑐)=𝐻^′₂(𝑐),𝐻^′′₁(𝑐)=𝐻^′′₂(𝑐),_𝑛𝑗=1𝛼_𝑚+𝑗𝜑_𝑚+𝑗𝑥_𝑚+𝑖−𝑥_𝑚+𝑗=𝑔𝑥_𝑖𝐻₂𝑥_𝑚+𝑖+𝑓𝑥_𝑚+𝑖𝐻+𝑟,𝑖=2,…,𝑛−1,₂(𝑑)=𝐻₃(𝑑),𝐻^′₂(𝑑)=𝐻^′₃(𝑑),𝐻^′′₂(𝑑)=𝐻^′′₃(𝑑),𝐻^′₃(𝑏)=𝛽. (34) Then, we can solve (34) to identify the numerical solution (32) to problem (3).

4. Numerical Experiment

We do numerical experiments to test the efficiency of our numerical methods.

Example 1. Consider the following second-order boundary obstacle problem as given in [22]: 𝑢^′′=𝜋0,0≤𝑥≤4,𝜋𝑢−1,43≤𝑥≤43𝜋,0,4𝜋≤𝑥≤𝜋, (35) with boundary condition 𝑢(0)=0, 𝑢(𝜋)=0. The exact solution of this problem is given as 4𝑢=𝑟₁𝜋𝑥,0≤𝑥≤4,41−𝑟₂𝜋cosh4,𝜋−𝑥43≤𝑥≤44𝜋,𝑟₁(),3𝜋−𝑥4𝜋≤𝑥≤𝜋, (36) where 𝑟₁=𝜋+4coth(𝜋/4), 𝑟₂=𝜋sinh(𝜋/4)+4cosh(𝜋/4).
The results are reported in Table 1. The errors between numerical and exact solution in terms of 𝐿_∞ and 𝐿₂ are listed in Table 1. The accuracy of the proposed method is measured using the 𝐿₂ and 𝐿_∞ error norms defined as 𝐿₂=1𝑛+1^𝑛_𝑖=0𝑢^exact_𝑖−𝑢^num_𝑖²,𝐿_∞=max_{0≤𝑖≤𝑛}𝑢^exact_𝑖−𝑢^num_𝑖. (37) In Table 1, we set 𝑥_𝑖=(𝑖−1)ℎ₁,𝑖=1,2,…,6,ℎ₁=𝜋,𝑥20_𝑖=𝜋4+(𝑖−6)ℎ₂,𝑖=6,7,…,11,ℎ₂=𝜋,𝑥10_𝑖=34𝜋+(𝑖−11)ℎ₁,𝑖=12,13,…,17,ℎ₁=𝜋.20 (38)

Table 1 

Comparison of exact and numerical solutions for Example 1.

Example 2. Consider the following third-order boundary obstacle problem as given in [17]: 𝑢^′′′(𝑥)=10,0≤𝑥≤4,𝑢(𝑥)1−1,43≤𝑥≤4,30,4≤𝑥≤1, (39) with boundary conditions 𝑢(0)=0, 𝑢^′(0)=0, 𝑢^′(1)=0.
The exact solution of this problem is 𝑢(𝑥)=12𝑎₁𝑥²1,0≤𝑥≤4,1+𝑎₂𝑒^𝑥+𝑒^−𝑥/2𝑎₃√cos32𝑥+𝑎₄√sin32𝑥,143≤𝑥≤4,𝑎₅𝑥12𝑥−1+𝑎₆,34≤𝑥≤1, (40) where 𝑎₁𝑎=0.14520742362098,₂𝑎=−0.21130240827197,₃𝑎=−0.78610085318732,₄𝑎=−0.24585768969643,₅𝑎=0.05860440434801,₆=0.04768241777632. (41)

The results are reported in Table 2. The errors between numerical and exact solution in terms of 𝐿_∞ and 𝐿₂ are listed in Table 2.

Table 2 

Comparison of exact and numerical solutions for Example 2.

From the results of Tables 1 and 2, we can see that the shape parameter 𝑠 plays no significant role in the computation.

5. Conclusions

In this paper, based on the multiquadric RBF and inverse multiquadric RBF, numerical schemes are presented for solving system of second-order and third-order boundary value problems associated with obstacle. The schemes are shown to be highly accurate and easy to implement. We point out three notes here. First, to solve second-order problem, we use 𝐿_𝐷 operator to approximate the derivatives on the interval [𝑎,𝑐] and [𝑑,𝑏] because by using 𝐿_𝐷 operator, we can on one hand avoid solving linear systems and on the other hand achieve high accuracy. Second, to solve third-order problem, we do not employ 𝐿_𝐷 operators to approximate the derivatives on the interval [𝑎,𝑐] and [𝑑,𝑏] because it is no easy to integrate 𝐿_𝐷𝑓 for three times and no explicit integral expressions can be found. Third, why we choose MQs and IMQs rather than other RBF? That is because MQs and IMQs can achieve high accuracy and IMQs are the second derivatives of MQs and this property bring much convenience to our calculation. Further study can be carried out for the applications of RBF in the numerical solution of higher-order boundary value obstacle problems.