#### Abstract

We introduce improved element-free Galerkin method based on block pulse wavelet integration for numerical approximations to the solution of a system of fourth-order boundary-value problems associated with obstacle, unilateral, and contact problems. Moving least squares (MLS) approach is used to construct shape functions with optimized weight functions and basis. Numerical results for test problems are presented in this article to elaborate the pertinent features for the proposed technique. Comparison with existing techniques shows that our proposed method based on integration technique provides better approximation at reduced computational cost.

#### 1. Introduction

Variational inequality theory has turned out to be an effective and potential tool for solution of mathematical models for obstacle, unilateral and contact problems. These problems arise in study of fluid flow through porous media, elasticity, transportation, optimal control, nonlinear optimization, and operations research [14]. The area of obstacle problems arising in fluid flow through porous media and elasticity forms an significant basis for the applications of variational inequalities. It has been revealed by Kikuchi and Oden [5] that the problem of equilibrium of elastic bodies in contact with a rigid frictionless foundation can be studied in the structure of variational inequalities. In a variational inequalities formulation, the location of the free boundary becomes an intrinsic part of the solution, and no special techniques are needed to locate it. Various numerical methods were employed to solve such type of problems [14]. If the obstacle is known, then the variational inequalities can be characterized by a system of differential equations by using the penalty function method of Lewy and Stampacchia [6]. In this paper, we apply element-free Galerkin method based on block pulse functions to develop a numerical method for obtaining smooth approximations to the solution of a system of fourth-order obstacle boundary-value problems of the type as follows: With boundary conditions, consider where and are continuous functions on and , respectively. The parameters , , and () are real constants and the continuity conditions of at and .

This problem was solved by many authors using different order of polynomial and nonpolynomial spline. Papamichael and Worsey [4] and Khalifa and Noor [3], have developed first- and second-order methods for solving such type of boundary value problems. Siraj-ul-Islam et al. [7] have established and analyzed smooth approximations for second, third, and fourth order linear/nonlinear boundary-value problems.

Galerkin-finite element method is one of the most popular well-developed methods, has been extensively used in the computational mechanics, due to its robustness, versatility, and convenience. However, the finite element method has inherent problems such as locking and poor derivative solutions [8]. FEM requires considerable amount of labor and time in generation of a predefined mesh/remesh of the problem domain. In structural analysis, the stresses calculated employing FEM algorithms are discontinuous and less accurate. When large deformations are analyzed, distortions in elements cause loss in solution accuracy in the case of large deformations [9]. Other numerical techniques such as the BEM (boundary element method) is limited by the availability of the infinite space fundamental solution for (at least the highest linear) differential operator of the problem.

The physical phenomena modeled by fourth-order obstacle boundary value problems are complicated, particularly, when they involve complex materials (e.g., viscoelastic materials with poisson ratio nearly 0.5) and where load is highly oscillatory. In conventional displacement-based approaches like finite element method (FEM), the interpolation of displacement requires continuity to ensure the convergence of FEM procedure for 4th order boundary value problems. This interpolation involves very intricate shape functions. These shape functions occupy large number of degrees of freedom in every element, including nodal rotations, displacement and so forth. For detailed discussion for solving such problems, the reader is referred to works [10] and references therein.

Improved element-free Galerkin (EFG) method provides an alternative choice of numerical tools which presents an attractive option in solution of such problems due to its flexibility in the absence of nodal interconnectivity. A meshfree (or meshless) method is a method used to establish system algebraic equations for the whole problem domain without the use of a predefined mesh for the domain discretization [9]. Several meshless techniques have been developed, for example, the “diffuse element method” [11] the “element-free Galerkin (EFG)” method [12] “meshless local Petrov-Galerkin (MLPG)” method [8], the “smooth particle hydrodynamics (SPH)” method, the “reproducing Kernel particle” method [13] “boundary node” method [14, 15] “boundary point interpolation” methods [16] and so forth [8].

#### 2. Numerical Technique

Meshfree method has been applied to physical problems using Gaussian weight function. Moving least square method is used to generate the shape functions. Galerkin weak form of the boundary value problem is formulated to give the system equations. Penalty/Lagrange method is used to impose the boundary conditions. Numerical integration of the system equation has been performed using block-pulse function.

##### 2.1. Meshless Shape Function

Shape functions for meshless techniques need to satisfy certain conditions such as adherence to partition of unity (), compact domain of influence, and adapting to randomness of nodes, to name a few. Computational efficiency is significantly affected by choice of shape function.

MLS. Moving least square (MLS) scheme was developed by mathematicians for interpolation, data fitting, and surface construction. Meshless techniques maintain the local character of the numerical implementation, by using a local interpolation to represent the trial function with the fictitious values of the unknown variable at some randomly located nodes. The local interpolation techniques generally used are MLS, partition of unity method, reproducing kernel particle method, hp clouds, Shepard function, and so forth. The moving least squares (MLSs) approximation has reasonably high accuracy and can be generalized to work with -dimensional problems. It provides continuous approximation for the field function over the entire problem domain [9, 17].

Here, moving least square (MLS) scheme is considered to interpolate the field variable. MLS approximation of a field variable function is defined as where = number of terms in the basis, = polynomial basis function,   = nonconstant coefficients, and .

is a complete monomial basis vector of order ; is a vector containing coefficients: , . The unknown coefficients of approximation are computed by minimizing the difference between the local approximation at a point and the nodal parameter for the node : .

Here the sample point may be a nodal point under consideration or a quadrature point. The support of the nodal point is usually taken to be a circle of radius , centered at . The weighted residual functional is evaluated as where is weight function associated with the node calculated at point . for all points in the support domain of node , and Minimization of weighted residual functional results in where is a non-singular matrix. The weight moment matrix is where is is the weight function matrix: The following Gaussian weight function is used in the present study: Shape perimeter () substantially affects the accuracy of the solution. Weight function for different values of shape parameters and its derivatives are illustrated in Figures 1 and 2 respectively.

Obtain value from (7), and substituting in (5) is the nodal shape function For any node , Transform the strong form of obstacle boundary value problems into symmetric variational or weak form as follows: where is the trial function and is the test function. We get the following linear equation system after employing “Integration by parts” technique and by using the shape functions and boundary conditions Solving the system of equation, we get that nodal fictitious values are obtained. Required numerical solution can be determined by multiplication of nodal fictitious values and shape functions :

##### 2.2. EFG

Element-free Galerkin (EFG) technique was introduced by Belytschko et al. [12]. EFG results in generation of nodes with variable interconnectivity, thus providing huge flexibility in discretization. Solution is expanded on a basis of MLS shape functions instead of piecewise polynomial basis as in case of the Finite element method. System equations are formulated using Galerkin weak form whereas integration is performed using background cells. EFG is conforming due to the use of consistent and compatible shape functions [9].

##### 2.3. Block-Pulse Function-Based Integration

The accuracy of numerical approximation is dependent on the integration technique. Quadrature rule is commonly used for numerical integration. In this method, interpolating polynomials are employed to find weights for nodes in the domain. Conventional Gaussian quadrature scheme involves large number of points for accurate numerical solution. In other numerical integration techniques such as Newton-Cotes quadrature rule, fine discretization leads to erroneous results due to high degree of polynomial interpolation [18]. Accuracy of Gaussian quadrature rule depends on the selection of nodes and corresponding weights. Method of undetermined coefficients may be used for quadrature rule which requires tabulation of nodes and weights before numerical integration.

In order to avoid these difficulties, Block-Pulse functions (BPF) have been used in this study to find numerical solution of system integral equations. Wavelet transforms were primarily developed for signal analysis, but they have also been used in applications like image compression, data compression, de-noising data and numerical approximation. Various types of wavelets include Battle-Lemarie, B-spline, Chebyshev, Daubechies, Legendre and Block-Pulse functions [18].

Numerical Integration based on wavelets is becoming popular in recent years [18]. In this study, Block-Pulse functions have been used due to their suitable properties and exceptional performance in numerical integration [18]. Their performance particularly, in approximation of highly oscillatory and improper integral is more accurate as compared to existing method in literature [1823].

Block-Pulse wavelet is one of the simplest wavelet. BPF with number of sets is defined as [18]: Here and and .

Let us consider the integrand , for , we get:

Theorem 1. The integral can be approximated as [18]

Proof. Consider To calculate the coefficients we consider the nodal points, The equation can be discretized to give Thus the integral can be approximated based on -set block-pulse function as or Optimal weights are calculated using built-in procedure in terms of wavelets. There is no need for tabulation of weights, and no intermediate technique is required for integration. BPF-based integration eliminates the need for quadrature nodes. Integrands are treated explicitly without solving nonlinear system which results from unknown nodes and weights.

##### 2.4. Consistency

A numerical method is called convergent if it approaches the exact solution as the discretization is refined. Consistency and stability guarantee convergence. The choice of monomial basis vector influences the consistency of the EFG method. MLS approximation reproduces all components that appear in . EFG method is consistent with order if all monomials up to order are included in [24]. If , this implies that is linear complete and if , is quadratic complete, and so forth. For consistency test we consider MLS approximation which reproduces all components of polynomials. Consider

#### 3. Application and Test Problems

We consider the following example fourth-order obstacle boundary value problem (BVP) for finding such that where , a given force acting on the beam and , is the elastic obstacle. Now the problem will be discussed in the framework of variational inequality approach. Define set , as is a closed convex set in , (Sobolev space), which is in fact a Hilbert space. It can be easily shown that energy functional associated with the obstacle problem (27) is it can easily proved that the form , defined by (30), is bilinear, symmetric, and positive (coercive) and the functional , is a linear continuous functional [13]. Minimum of the functional defined by (29) on the closed convex set in can be characterized by the variational inequality Thus we conclude that the obstacle problem (27) is equivalent to solving variational inequality problem (32). This equivalence has been used to study the existence of a unique of (27); see [13], and using the idea of Lewy and Stampacchia, [6] can be written as where is the obstacle function and is a penalty function defined by We assume that the obstacle function, , is defined by From (33)–(35), we obtain the following system of equations With the boundary conditions, consider

Problem 1. In this example, we consider system of differential equation (38), when The analytical solution for this boundary value problem is where ,  ,  . Numerical and analytical solution of Problem 1 is shown in Figure 3.

Problem 2. For , With the same boundary conditions as in previous example, consider Numerical and analytical solution of Problem 2 is shown in Figure 5.

Test problems were solved for different number of nodes and constant basis with the help of improved EFG based on numerical integration using Block pulse functions/wavelets (BPF). Figures 4 and 6 show the errors in solution of test Problems 1 and 2, respectively. Maximum absolute errors in solution are tabulated and compared with the results of exact solution, which show better accuracy.

##### 3.1. Results and Discussion

Maximum errors obtained using our technique based on block pulse function (BPF) for our test problems are shown in Tables 1 and 2. Comparison of maximum absolute errors with exact solution demonstrates the accuracy of our method. Results using quadrature rule-based integration are also calculated to show the performance of BPF-based integration technique. Maximum absolute errors for test problems are shown in Tables 1 and 2. It is clearly seen that integration of system equations using BPF has produced better results as compared with quadrature rule method for numerical integration.

Effect of Integration Points. It has been observed, that the accuracy of meshless computational technique depends upon the size of influence of domain, penalty factor, and Integration points. In this paper, we have studied that the accuracy of numerical method generally improves as we increase the number of integration points; however, if the number of gauss points exceeds some limiting value, it may have some adverse effects on the accuracy. It has been determined that arrangement of background cells and ratio of total number of Gauss points to the total number of nodes, plays very important roles in accuracy. To avoid all these complications to improve the accuracy of the meshless method, we purpose new numerical integration technique as mentioned above.

#### 4. Conclusion

In this study, improved numerical technique “element free Galerkin (EFG) technique using Block-Pulse function-based integration” has been presented. The results obtained by the suggested method for given problems as mentioned ((38) and (40)) exhibit its ability to provide improved solutions. The improved EFG technique showed fast convergence and provided better results at reduced number of nodes, as mentioned in Tables 1 and 2. It has capability of handling of improper, highly oscillatory function. Details of the results are shown the figures and tables.