Adaptive Finite Element Method for Fredholm Integro-Differential Equation

In this work, we apply an adaptive finite element method for solving Fredholm integro-differential equations. The adaptive refinement is based on a recovery-type error estimator. The results of some numerical examples are shown to illustrate that prescribed accuracy is achieved with minimum degrees of freedom.


Introduction
An integro-differential equation (IDE) is an equation in which differential and integral operators appear at the same time.The mathematical formulations in many applied areas of sciences and engineering contain integro-differential equations.These include the study of fluid, physics, chemistry, biology, mechanics, astronomy, potential theory, electrostatics and chemical kinetics [1]- [3].Also some partial differential equations are converted into IDEs as in [4], where a steady state diffusion equation which contains two unknowns is replaced by a single integro-differential equation containing only one unknown.IDEs with constant limits of integration are called Fredholm IDEs.There are varieties of numerical methods that have been developed to solve Fredholm IDEs.Some of these methods are direct method based on Fourier and block-pulse functions [5], compact finite difference method [6], Taylor and Chebyshev matrix methods [7]- [10], an extrapolation method [11], method of regularization [12], Adomian decomposition method [13], and series solution method [14].One of the most powerful methods utilized to solve IDEs is the finite element method (FEM).FEMs are numerical methods that successfully obtain good approximate solutions to different types of equations.FEMs have many subclasses where the C0-Galerkin FEM is considered the simplest of them.One of the main features that differ FEMs from other numerical methods is the combination of variational methods and piece-wise polynomial approximation.Also, FEMs deal with complicated domains with stability and flexibility.FEMs, as they are known today, were introduced in 1960 [15].Different forms of IDEs were solved using FEMs (see for example [16]- [23]).Adaptive procedures for the FEM approximate solutions are important techniques that enhance the accuracy and the efficiency of finite element solutions.They perform refinement only in parts of the domain where the error in the solution is large by using error estimators.A priori and a posteriori error estimates for the optimal control problems governed by integral or IDEs of Fredholm type were established in [23].In this article, an adaptive FEM procedure is proposed to solve Fredholm IDE of the second kind of the following form subject to boundary conditions: u(a) = d 1 and u(b) = d 2 .In this adaptive procedure, a recovery is constructed for the FEM solution using superconvergent patch recovery (SPR) technique.The mesh is refined to obtain higher accuracy solution till the error calculated between the FEM solution and the recovered solution reaches a prescribed tolerance.

Finite element solution
At first, the weighted residual is formed as where R denotes the residual which is expressed as where u is the approximate solution and w is a weighting function to be chosen.Substituting In this step, the derivative inside the integral is reduced from second order to first order by integrating by parts.Thus (2.4) becomes Now, we divide the domain into finite number of elements and let the FEM approximate solution u h be given in the form where φ j (x) are the basis functions, c j are the coefficients to be determined and h denotes the maximum subinterval length.
Choosing w i (x) = φ i (x) (The Galerkin assumption), so that the residual is kept orthogonal to the space of functions used in the approximation of u.By substituting (2.6) in (2.5) and rearranging, we get We obtain a system of equations which has the following matrix form where ) ) ) Applying boundary conditions gives directly the first and last nodal values: i.e., u 1 = d 1 and u m = d 2 .These Dirichlet boundary conditions replace the first and last rows in the global system of equations (2.8).Solving these equations for C, the approximate solution u h (x) is obtained.

Solution recovery
We are concerned with the matter of enhancing the accuracy of the approximate solution by a posteriori treatment of the data.Such process is referred to as recovery technique and the error estimator is classified as a recovery-based error estimator.Here, the SPR technique is utilized to produce superconvergent displacements.We use the recovery procedure presented in [24].This procedure differs from the general gradient recovery technique proposed in [25] and [26] in two things.Firstly, recovery is made for the whole element instead of each node and an element patch is used instead of a node patch.Secondly, the local error is calculated directly over the element from the recovered solution for this element.First, the element patch is constructed which consists of the current element e and the adjacent elements sharing common nodes with it.For each element, the recovered solution R h u h is expressed as this is a single and continuous polynomial representation of the function describing the solution over the element where a denotes the unknown parameters and Q(x) is the basis of the assumed polynomial for R h u h .Let the superconvergent finite element solutions in this patch be S h and the number of S h be ns.To determine the parameters a, the expansion of R h u h is made to fit the superconvergent points on the element patch in a least squares manner.Thus we have to minimize the functional This leads to a = A −1 b, (3.18)where After calculating the parameters a, the local error (over the element e) is calculated directly from the recovered solution.
One of the major applications of the recovery methods is its use in the calculation of a posteriori error estimators.
Having the recovered solutions, the error can now be evaluated simply by replacing the exact value of u, which is in general unknown, by the recovered value R h u h which is better than the direct finite element solution.Here we choose the L 2 norm so the error estimator ε h can be written as follows To ensure the accuracy and the quality of the error estimator, the effectivity index θ is calculated.The effectivity index is the ration between the recovery error estimator and the exact error ε h , i.e.
The error estimator is said to be asymptotically exact if θ approaches unity as the mesh partition size approaches zero.

Adaptive procedure
In adaptive techniques, iterative procedures are performed starting with an initial mesh until a condition that involves a predefined tolerance ε is satisfied.In this work, the SPR technique is applied at every iteration and the local error ε h,e between the recovered solution and the finite element solution is calculated for every element.The maximum of all local errors is checked to see if it is less than the predefined tolerance ε.If this is the case, the adaptive procedure is put to an end.Otherwise, we mark some elements to be refined, and repeat the adaptive procedure.The marking strategy also depends on the local error (ε h,e ) at each element as follows.A ratio τ is to chosen such that 0 < τ < 1.If (ε h,e k > τ(max i (ε h,e i ))), the k th element is marked for refinement.The algorithm used is summarized in the following steps: 1. Define the initial mesh and the desired tolerance ε.
2. Calculate the finite element solution u h .
3. Apply solution recovery using SPR and determine the local error ε h,e for each element.

If (max
) then mark element k, refine marked elements and go back to step 2.

Numerical examples
The following examples are solved by FEM using quadratic element approximation then we use fourth order polynomial for the recovered solution.The tolerance is chosen as ε = 10 −6 .Figure 1 shows the good agreement between the exact solution and the recovered solution.It also indicates that this accuracy is achieved with minimum degrees of freedom (DOFs) as the refinement is carried out only at the parts of the domain where the solution exhibits fast changes.This is confirmed by the graph of final mesh in Figure 2. Figure 3 shows that the effectivity index approaches unity after several iterations.Figure 4 and Figure 5 confirm that good agreement between the exact solution and the recovered solution is obtained with minimum DOFs and that the refinement is performed only at parts of the domain where the solution exhibits fast changes.Figure 6 shows that the effectivity index approaches unity after several iterations.

Conclusion
An adaptive finite element scheme is presented to solve Fredholm IDEs.The scheme is based on using a recovery based error estimator to choose the parts of the domain to be refined.The recovery technique utilized is a function recovery type which is based on SPR technique.The numerical examples illustrate the advantage of this scheme as the it is successful in marking the elements where the solution exhibits fast changes and refine them only, hence achieving high accuracy with minimum DOFs.

Figure 4 :
Figure 4: Exact solution versus recovered solution for Example 5.2.