A Note on the Finite Element Method with Singular

Basis Functions

Hideaki Kaneko ^{∗} Department of Mathematics and Statistics Old Dominion University Norfolk, Virginia 23529-0077

Peter A. Padilla Crew Vehicle Integration Branch Flight Dynamics and Control Division NASA-Langley Research Center MS-152 Hampton, VA 23681

Abstract

In this note, we make a few comments concerning the paper of Hughes and Akin [?]. Our primary goal is to demonstrate that the rate of convergence of numerical solutions of the finite element method with singular basis functions depends upon the location of additional collocation points associated with the singular elements.

Key Words: Finite element method with singular basis functions, the collocation method, Singular interpolation elemnts.

∗This author was partially supported by NASA under grant NCC1-213 when much of the research presented in this paper was done.

1 Finite Element Method With Singular Basis Functions

In the paper [?], Hughes and Akin made an interesting observation concerning the finite element analysis that incorporates singular element functions. A need for introducing some singular elements as part of basis functions in certain finite element analysis arises out of the following considerations. The solution of certain problems, such as a field problem [?], exhibits highly singular behavior due to geometric features of the spatial domain. On the other hand, in other circumstances, the solution is overwhelmingly affected by the nature of loading and the problem of singularity can be ignored. To satisfy both situations just described, it is thought that an incorporation of singular elements that emulate the solution with the standard polynomial elements may perhaps be desirable. This is the point that was exploited in [?] by Hughes and Akin. In order to make the computations of the finite element method with singular elements more efficient, they consider the following algorithm for constructing interpolation functions. A construction of such algorithm was motivated by the idea that “it is of practical interest to develop techniques for systematically defining shape functions for singularity modeling (and for developing special elements in general), which circumvent the interpolation problem” ([?] p.176). The algorithm that they developed goes as follows:

ALGORITHM Suppose that there are n shape functions N_{a}, a = 1, 2, …, n which satisfy the interpolation property on the first m nodes r_{b}, viz., N_{a}(r_{b}) = δ_{ab}, a, b = 1, 2, …, m (m < n). Their idea here is to reshape N_{a}’s so that the interpolation property is satisfied on all n nodes. The algorithm is given by

Step 1 N_{m+1} ← \frac{N_{m+1}(r)−∑ a=1^{m}N_{m+1}(r_{a})N_{a}(r)}{N_{m+1}(r_{m+1})−∑_{a=1}^{m} N_{m+1}(r_{a})N_{a}(r_{m+1})}
Step 2 N_{a}(r) ← N_{a}(r) − N_{a}(r_{m+1})N_{m+1}(r),  a = 1, 2, …, m
Step 3 If m + 1 < n, replace m by m + 1 and repeat Steps 1 to 3.
If m + 1 = n, stop.

Before we comment on the algorithm, we would like to draw the reader’s attention to the recent paper by Dydo and Busby [?]. They observed that “using the algorithm as a subroutine during the finite element shape function generation requires re-evaluation of all terms of the shape functions at each nodal point and eventually leads to excessive round-off error”. To prevent this difficulty, they develop an equivalent algorithm in matrix form to generate finite elements with special properties. The results reported in [?] demonstrate a considerable increase in computational efficiency.

We now return to our discussion of the algorithm. To demonstrate this algorithm, we borrow one of the examples from [?]. Let r_{1} = 0, r_{2} = \frac{1}{2} and r_{3} = 1. The shape functions that we reconstruct are N_{1}(r) = 1 − 2r, N_{2}(r) = 2r and N_{3}(r) = r^{α} where α representing some real number. Note that N_{a}(r_{b}) = δ_{ab}, 1 ≤ a, b ≤ 2. An application of the above algorithm gives

\begin{aligned}N_{1}(r) & ← 1 − 2r + [\frac{r^{α}−2( \frac{1}{2} )^{α}r}{1−2( \frac{1}{2} )^{α}} ] \\ N_{2}(r) & ← 2r − 2 [\frac{r^{α}−2( \frac{1}{2} )^{α}r}{1−2( \frac{1}{2} )^{α}} ] \\ N_{3}(r) & ← \frac{r^{α}−2( \frac{1}{2} )^{α}r}{1−2( \frac{1}{2} )^{α}} . \end{aligned}

Of course, the newly defined shape functions satisfy

N_{a}(r_{b}) = δ_{ab}  1 ≤ a, b ≤ 3.  (1)

It is important to note at this point that the algorithm is subject to the location of the interpolation points r_{b}, m + 1 ≤ b ≤ n. Clearly, step 1 of algorithm is sensitive to the location. Namely, if these interpolation points are such that

N_{m+1}(r_{m+1}) − ∑_{a=1}^{m} N_{m+1}(r_{a})N_{a}(r_{m+1}) = 0  (2)

then the algorithm does not work. Out of this observation, there seems to arise a profound and difficult problem in the area of approximation theory. The problem is important in that the success of the finite element method using the collocation scheme with singular basis functions hinges on a resolution of this problem. To describe it, let W_{p}^{k} denote the Sobolev space,

W_{p}^{k} = {f|f^{(k)} ∈ L_{p}(Ω)}

where f^{(k)} denote the kth generalized derivative and Ω is a bounded region in R. The theory extends easily to higher dimensions. Now let U ≡ span[N_{a}]_{a=m+1}^{n} where {N_{a}}_{a=m+1}^{n} represent the singular elements. Also denote by S_{h}^{k} the approximation space spanned by {N_{a}}_{a=1}^{m}. Here S_{h}^{k} is usually taken as the space of piecewise polynomials of degree k − 1 with length of each subinterval h. Our goal is to approximate each element of U ⊕ W_{p}^{k} by an element from U ⊕ S_{h}^{k} by interpolation. That is, if {r_{b}}_{b=1}^{m} is such that N_{a}(r_{b}) = δ_{ab} for 1 ≤ a, b ≤ m, then for each f ∈ U ⊕ W_{p}^{k}, we must find v ∈ U ⊕ S_{h}^{k} that satisfies

v(r_{b}) = f(r_{b})  b = 1, …, n,  (3)

where {r_{b}}_{b=m+1}^{n} are specified interpolation points for {N_{a}}_{a=m+1}^{n}. Denote the interpolation projector of U ⊕ W_{p}^{k} to U ⊕ S_{h}^{k} by P_{h}. Namely P_{h} is defined so that

P_{h}f(s) = v(s)  s ∈ Ω.  (4)

Notice that P_{h}^{2} = P_{h}. In order to achieve a convergence by the collocation scheme in the finite element method, we must examine the following inequality. Here we assume that the number of singular basis functions, n − m, is fixed.

|det[N_{a}(r_{b})]|_{a,b=1}^{n} ≥ ϵ > 0  for all n.  (5)

This inequality is a necessary and sufficient condition for the algorithm of Hughes and Akin to work. It is important to remark at this point that the success of algorithm depends upon the existence of a solution to the interpolation problem (3) which in turn is equivalent to the condition (5). For each fixed index n (hence for fixed m) , it is not difficult to find n − m interpolation points r_{b}, m + 1 ≤ b ≤ n, that correspond to the singular basis elements N_{a}, m + 1 ≤ a ≤ n for which the inequality in (5) is satisfied. What is difficult here is the question of locating n − m points for as many singular elements for which condition (5) is satisfied for all n. The problem of finding n − m interpolation points for singular basis functions that work for all n is currently under investigations. In the following section, we proceed our discussion of the interpolation problem which leads to a finite element analysis with singular elements. The discussion will provide information concerning the rate of approximation of interpoalation. Interpolation examples at the end of the section show that rates of convergence are quite sensitive to locations of interpolation points for singular elements.

2 Convergence Analysis

When condition (5) is satistied, one can deduce the rate of convergence of the projector P_{h} to the identity operator I. As is well known -e.g., [?], the convergence rate of such interpolation projectors determine the rate of convergence of the finite element method that uses collocation scheme. The following theorem of Cao and Xu [?] is useful. We sketch a proof for completeness.

Lemma 2.1 Let X be a Banach space. Assume that U_{1} and U_{2} are two subspaces of X with U_{1} ⊆ U_{2}. Moreover assume that P_{1}: X → U_{1} and P_{2}: X → U_{2} are linear operators. If P_{2} is a projection, then

∥x − P_{2}x∥_{X} ≤ (1 + ∥P_{2}∥_{X})∥x − P_{1}x∥_{X} for all x ∈ X.

Proof: Let x ∈ X. We write

x − P_{2}x = (x − P_{1}x) + (P_{1}x − P_{2}x).

Since P_{1}x ∈ U_{1} and U_{1} ⊆ U_{2}, we have P_{2}P_{1}x = P_{1}x. Hence,

x − P_{2}x = x − P_{1}x + P_{2}P_{1}x − P_{2}x
= (I − P_{2})(x − P_{1}x).

It follows that

∥x − P_{2}x∥_{X} ≤ (1 + ∥P_{2}∥_{X})∥x − P_{1}x∥_{X} for all x ∈ X.

□

This lemma seems to suggest a reason for which the method of interpolation using a set of singular elements is sensitive to the location of interpolation points corresponding to the singular basis functions. We let P_{1} be the interpolatory projection of U ⊕ W_{p}^{k} onto S_{h}^{k}, namely for each f ∈ U ⊕ W_{p}^{k}, P_{1}f ∈ S_{h}^{k} is defined by

P_{1}f(r_{a}) = f(r_{a}), a = 1, …, m.

Theorem 2.2 Assume that y = u + v with u ∈ U and v ∈ W_{p}^{k}. Let P_{h} be the projection defined by (4). Then

∥P_{h}y − y∥_{p} ≤ (1 + ∥P_{h}∥_{p})∥y − P_{1}y∥_{p}

where 1 ≤ p ≤ ∞.

We make two observations here. First, condition (5) plays the essential role for the pointwise convergence of P_{h} to the identity operator I, which, in turn, by the uniform bounded principle, guarantees the uniform boundedness of ∥P_{h}∥ for h > 0. Second, because of a singular component present in y, it is expected that ∥y − P_{1}y∥_{p} does not converge at the optimal rate. Two observations seem to provide some explanations to the phenomenon revealed in the example below. A comprehensive study of the power of approximations by splines can be found in [?].

A similar study of incorporating the singular elements into an approximating basis was given in the numerical analysis for the weakly singular Fredholm equations. The results from [?] were recently generalized to a class of nonlinear Hammerstein equations in [?], which in turn extends the work in [?]. These papers are concerned with the singularity preserving Galerkin methods and the analysis associated with these methods appear more straightforward. We close this note by demonstrating the sensitivity of the location of interpolation points for singular elements.

EXAMPLE: Let f(x) = \sqrt{x} + \sqrt{1 − x} + x^{2}. We wish to approximate f over [0, 1] by an element from U ⊕S_{h}^{2}, where U = span[\sqrt{x}, \sqrt{1 − x}]. Let \{x_{i}\}_{i=0}^{n} be the uniform partition of [0, 1] defined by x_{i} = \frac{i}{n}, i = 0, 1, …, n and h = \frac{1}{n}. The interpolation points used to define an element from S_{h}^{2} are taken to be the zeros of the second degree Legendre polynomials transformed into [x_{i−1}, x_{i}] for i = 1, 2, …, n. The following data shows that (a) when the interpolation points for the singular elements are taken to be t_{1} = \frac{1}{5} and t_{2} = \frac{4}{5} for each n, the convergence is O(h^{1/2}), whereas (b) when t_{1} = \frac{1}{2^{n}} and t_{2} = 1 − t_{1}, then the convergence is of the order O(h^{2}).

interpolation point t_{1} =  1/2^{n}  1/5
n = 4  0.0513168  0.0513168
n = 6  0.0008437  0.0235540
convergence exponent =  1.99  0.56

Table 1. Error and convergence rate data for the example

References

[1] R. E. Barnhill and J. R. Whiteman, Error Analysis of Finite Element Methods with Triangles for Elliptic Boundary Value Problems, Mathematics of Finite Element Methods and Applications, ed. J. R. Whiteman, London: Academic Press, (1973), 83-112.

[2] Y. Cao and Y. Xu, Singularity Preserving Galerkin Methods for Weakly Singular Fredholm Integral Equations, Journal of Integral Equations and Applications, 6(1994), 303-334.

[3] J. R. Dydo and H. R. Busby, An Improved Algorithm for Developing ‘Special’ Finite Element Shape Functions, Communications in Numerical Methods in Engineering, 13(1997), 687-693.

[4] T. J. R. Hughes and J. E. Akin, Techniques for developing ‘special’ finite element shape functions with particular reference to singularities, International Journal for Numerical Methods in Engineering, 15(1980), 733-751.

[5] T. J. R. Hughes, The Finite Element Method, Linear Static and Dynamical Finite Element Analysis , Prentice-Hall Inc. Englewood Cliffs, N. J. (1987).

[6] H. Kaneko, R. Noren and P. Padilla, Singularity Preserving Galerkin Method for Hammerstein Equations with Logarithmic Kernel, Advances in Computational Mathematics(in press).

[7] H. Kaneko and Y. Xu, Superconvergence of the Iterated Galerkin Method for Hammerstein Equations, SIAM Journal of Numerical Analysis,33(1996), 1048-1064.

[8] L. Schumaker, Spline Functions Basic Theory , John Wiley and Sons, N.Y., 1981.