A least squares finite element method applied to b-splines

A Least Squares Finite Element Method Applied to B-Splines by

N.

G.

ZAh4ANI

Department of Mathematics Ontario NlG 2W1, Canada

and Sratistics,

University

of Guelph,

Guelph,

In this article we consiruct a class of functions on a bounded irregular region Rc RZ which are of compact support, smooth and locally polynomials. The basic tool is the use of ordinary B-splines associated with the square containing Cl. This construction is then used for approximating the solution of Poisson’s boundary value problem. 7’he approximation is carried out through a least squares finite element method applied to the above class. Aside from some computational experiments, the objective is to emphasize the ease of generating the basis elements and the role of Kernel function in a convergence proof. ABSTRACT:

I. Introduction Given the interval [a, b], a uniform mesh is imposed functions are constructed, (h = (b - a)/N, N is the number i=-1,0 ,..., N,N+l; ‘(x -xiJ3 h3+3h2(x~~-~)+3h(x-x~_~)2-3(x-~_~)3 Bi(x)=~~h3+3hZ(~+I-x)+3h(~+I-x)2-3(~+I-x)3

and the following of subintervals) for x E [xi--29 %-II X ~[q_lp

Xi]

x E [xi, %+I1

(xi+2-x)3

x E [%+I, xi+21

0

otherwise.

I? -splines These functions are known as the cubic B-splines B-splines are defined as the tensor products of one splines i.e. E&(x, y) = &(x)B,(y) see [S]. Restriction

[l],[S].The bicubic dimensional

cubic

B-

of B-splines.

We consider embedded in a single indexed sociated with S R. The schematic all the B-splines

fit, a bounded region in R* with piecewise smooth boundary, square S with a two dimensional grid system. We construct a family 9 = {Bi}E”=, where $&‘i’s are the bicubic B-splines aswhich have a nonzero restriction on the interior of the region construction of 9 is given in Fig. 1. In step 3 of this diagram which are associated with dotted mesh points are retained as

0 Tlw Franklin Institute 0016-0032/81/030195-14SO2.00/0

195

N. G. Zamani

step I h

h

step 3 FIG. 1.

elements of 9. One can prove that functions in 9 are linearly independent on a. The proof of independence is straightforward and is quite similar to that of independence on rectangular regions. For this reason we shall skip the proof. II. Application The class of functions constructed in the previous section will be used for a least squares approximation to the solution of certain boundary value problems. We deal with the simplest case, Poisson’s equation with one of the following well known boundary conditions: i-k, A=-$+$

Au=f

8R (Dirichlet)

u=g

au

Xl (Neumann)

z=g

E+f_ru=gaR

(Mixed)

o>O.

Here we assume that the functions f, g and u are smooth enough to guarantee the existence of a classical solution belonging to C”(n). Another major assumption that will be made throughout is that the solution u can be extended to a square S containing CR, and that the resulting extension belongs to C4(S). Such extension results are available, see e.g. A. Friedmann, Partial Differential Equations. In the case of Neumann’s problem, the compatibility condition

II I g f=

n

196

Jn

A Least Squares

Finite Element

must hold and in such a case the solution following two conditions,

I an

u = 0

or

can be made

particular

Applied unique

to B-Splines by one of the

(x,, y,) E 0’.

u(x,, yO) = 0

The reason we have chosen these kernel functions for such problems.

Method

constraints

is the availability

of

Algorithm Let us concentrate

on the Dirichlet

problem,

Au=f

fi

u=g On C’(0)

we define

an inner

product

Next we look for an approximation

aR of the form

of the form UM = f

@ai

i=l

and determine the coefficients so that the quantity 11~- u~\&~, is minimized. (For other treatments of least squares applied to boundary value problems, see (3), (s), (6), (9)). This leads to the following

set of normal

equations:

,El q((%, sjsi))= ((sjp~1) But the right hand

side is simply

an

n

((Bjai, U>> = (&si, f)bcn, + (4, Hence

the final form of the normal

E ai((%,

i = 1, . . ., MS

Se,)) = (&,

an

n g)b
equations

is

f)&(n) + (Bjai,g)bcaa,

j = 192, . . .y M.

i=l

Pointwise convergence We first show convergence in ().I)((,))norm. Let ii be the spline of interpolation (8) to u (after u has been extended to S) which corresponds to the square Vol. 311, No. 3, March 1981 printed in Northern Ireland

S

197

N. G. Zamani with its grid system. We then have the following string of inequalities,

=

JJ J R (Au-Azi)*+

(u-ii)*=

0(h4)+O(h8)=

0(h4).

Jll

Here we have used the fact that lID% -D%lj,= O(IZ~-“~), see (8). Hence convergence in 11&,, is established. With the aid of kernel functions from potential theory, we will be able to show pointwise convergence of u,. More precisely, let G(x, y; & n) be the Green’s function of the Dirichlet problem. We then have the following integral representations for u and u,: ub, Y> =

JJ JJ Gb,

n

k&,

Y) =

Gb,

Jil

Y; 6 q) An&,

aG an,.,

J J

Y; 5, rl)f(& rl) +

v)+

n

Taking the absolute value of the difference we get

g(4,11) (x, Y)E no,

aG 45, an,., M

s)

x, YEfiO.

and using the Schwartz inequality

JI$+MI 611

IU(X,Y)-U~((X,Y)I~JJIG(~,Y;~,~))IIAU-AL(M~+

n

JCl

or alternatively, lu(x, Y)- u&,

Y)I5 (K,(x, Y)+ K2(x, Y)) II” - %&o, = K(x, y)O(h2)

(x7 Y)E OR0

where ~,(x~)-Illl~(x,~:f~~l~~

K2k

Kb,

198

t

(*)

Y) = Y I=

K,k

Y I+

K2k

Y). Journs

of ‘l-he

Franklin lmtitute

P~RCE3J.M.

A Least Squares Finite Element Method Applied to B-Splines From the nature of Green’s function one can conclude that as long as (x, y) E rR”, the quantity K(x, y) remains finite. It is known (2) that both Neumann and Mixed equations possess kernel functions which admit similar integrability conditions to (*) in the Dirichlet problem. Hence we can apply the same algorithm to these boundary value problems and obtain convergence results. However the inner products have to be modified, the suitable inner products being given below; Mixed problem:

((u, u))= IlAu R

Au+ j (g+,,)(g+,) an

((u, u)) = IjAu R

Au + j g E+

Neumann 1

Neumann 2

We have established Theorem

14x,, y,)u(x,, yo).

an

the following result:

I

If the true solution u to the Poisson’s equation Au = f under any of the mentioned boundary conditions belongs to C”(Q) and has a C4 extension to a square S containing a, the least square algorithm gives an approximation u, with the following asymptotic error estimate: It&, Y)-4x,

y)I=G,

yP(h2)

(~3 Y)ERO.

Convergence

is uniform on compact subsets of R”. Again we first specialize to the Dirichlet problem. The following lemma (10)will be used to establish mean square convergence. Lemma. Suppose u E C’(n), then there exist constants y and S independent of u such that L,

convergence:

If we apply this lemma to u-u,, we can conclude that \Iu - u,(jL2cn, = O(h’). Similar types of a priori inequalities are available for Neumann and mixed problems. Hence we have established the following result. II Under the same assumptions as in Theorem I, the approximation u, converges to u in the L, norm and the following asymptotic error estimate holds: Theorem

11~Vol. 311, No. 3, March 1981 printed in Northern Ireland

~,ll~,cn, = 0(h2h 199

N. G. Zamani III. Numerical

Computation

This section is concerned with gorithm described for the Dirichlet originally defined as

the numerical implementation problem. Although the inner

((u, v)) = jjA~

of the alproduct was

Av + 1 uv,

n

JSI

it has been observed that introducing a weight function p(h) for the second integral accelerates the convergence dramatically, so the modified inner product is

((u, v)) = jI Au Av +p(h) n

I

UV.

JR

In all examples to follow it is assumed that p(h) = l/h4. We would like to say a few words about solving the Grammian system. In our computations a straight-forward Gaussian elimination is used. Cholesky’s method was also tried but did not prove to be superior to our naive scheme. In the worst case, a square region, the Grammian is (N + 3)’ X (N + 3)’ where N is the grid size. On a triangular region the Grammian is (N+3)(N+2)X(N+3)(N+2) 2

2

.

The dimensions considerably reduce in the air foil problem. These indicate that efficient algorithms should be employed to handle the Grammian matrix. Finally we briefly discuss the quadrature problem. Region R is approximated by a polygonal region fi, which in turn is subdivided into smaller subregions (Fig. 2). Hence we can write,

A simple trapozidal rule is used for computing the boundary different types of quadratures, depending on the shape of fii, approximate the double integrals.. In the case where Ai is a square, points {Pi},P_, are chosen as normalized Legendre polynomial of degree two in each direction

200

integrals but are used to the zeros (Fig. 3).

Jo,nnal of The

of

FrmJdin Institute PergammRessLtd.

A Least Squares Finite Element Method Applied to B-Splines Discretizotion the region

of

s

FIG. 2.

Hence the quadrature

is given by,

If fii is a polygon of general shape (Fig. 3), we let {Pi}:‘: (I = 3,4,5) corner nodal points together with the centroid of Ai (Fig. 3). The quadrature rule is given by, j jf=Area

(A,,[?

f(P,))/(I+

be the

1).

j=l

6,

Although this is a very crude quadrature, the results are satisfactory. Example 1. The equation under consideration is Au=-1 { u=o where 1R= [0, l]x[O, 11. This equation

IR aCI arises from the torsion of a bar with p2

PI .

b .

p4 .

P3 .

El p2

Nodal selection

in region 8,

FIG.3. Vol. 311, No. 3, March 1981 Printed in Northern Ireland

201

N. G. Zamani TAFJLE 1 h

l/3 l/4 l/5 l/6 l/7 l/8 119 l/10

LSFE

TAJJLE~

FE

FD 0.33E-2 O.l5E-2 0.88E-3 0.57E-3

0.51E-3 0.60E-4 0.60E-4 0.20E-4 0.30E-4 0.20E-4 O.lOE-4 0.30E-4

LSFE

h

0.70E-2 0.40E-2 0.26E-2 O.l6E-2 O.l3E-2 0.93E-3 0.81E-3 0.60E-3

l/3 l/4 l/5 l/6 l/7 l/8 l/9 l/10

O.l4E-1 0.80E-3 0.22E-2 0.20E-3 0.40E-3 0.20E-3 0.30E-3 O.lOE-3

square cross section. There is no closed form solution, however an accurate approximation is ~(0.5, 0.5) - 0.7366 which is obtained by taking a partial sum of the series solution. In Table 1 we have compared the above number with the approximations derived from three different numerical methods. These methods are: (1) LSFE-Least Squares Finite Element, (2) FD-Finite Difference (4), (3) FE-Classical Finite Element (7). In our classical finite element computations we have utilized the linear and bilinear interpolation functions, this is a handicap when compared to LSFE in which bicubic polynomials are used. In Fig. 4a, the logarithm of error is plotted against N. We like to mention that in the FD column, the point (0.5, 0.5) does not coincide with a mesh point when N is odd, so an approximation is not available. Another quantity of interest (related to the stress in a bar under torsion) is

AlIuuuL 3 4

5

6 7

6

3ILuu-u 3

4

5

6

7

6

9

IO

9

logarithm

3

5

6

IO of error

FIG. 4.

202

4

against

N

7

6

9

IO

A Least Squares Finite Element Method Applied to B-Splines TABLE 3

TAESLE~

h

LSFE

FD

FE

h

l/3 l/4 l/S l/6 l/7 l/8 l/9 l/10

0.67E-6 0.48E-6 0.20E-5 O.l9E-5 0.21E-5 0.55E-5 O.l5E-4 0.29E-4

0.37E-7 0.45E-7 0.86E-7 O.l4E-6 0.23E-6 O.OOE-0 0.96E-7 0.28E-6

0.49E-2 0.32E-2 O.l9E-2 O.l4E-2 0.98E-3 0.78E-3 0.60E-3 0.50E-3

l/3 l/4 l/5 l/6 l/7 l/8 l/9 l/10

LWE O.l3E-15 0.97E-16 0.43E-15 0.26E-15 0.74E-15 0.91E-15 0.40E-14 0.80E-14

@ulax)l~.o,o.s. For comparison purposes we will take the figure -0.3378 known in the engineering literature. Table 2 gives the error ]&&3x +0.3378] where u, is the approximate obtained from LSFE. Example 2. The equation to solve is Au =2(x*+y*)-2(x+y) { u=o

fi

a0

where 0 = [0, l] x [0, 11. The exact solution to this equation is u(x, y) = xy(1 -x)(1 - y) which is a bicubic spline. Table 3 gives the distribution of error for the three methods. In Table 4 we have recorded the error when all computations are carried out in double precision. The enormous gain of accuracy can be attributed to the fact that the solution is bicubic spline and that round off is the dominant source of error. For graphical representation we refer to Fig. 4b. Notice that since the exact solution is known, the error is defined by error = Example

max

6.y)mesh

points

lu(x, Y)-%Ax,

Yl.

3. We will next look at Au = 2 sin my - m*(x’- x) sin my fl u=(x*-x)sin7ru aci TABLE

Vol. 311, No. 3. March 1981 Printed in Northern Ireland

5

h

LSFE

Fm

FE

l/3 l/4 l/5 116 l/7 l/8 l/9 l/10

0.26E-3 0.85E-4 0.42E-4 0.23E-4 O.l8E-4 0.27E-4 0.54E-4 O.llE-3

0.93E-2 0.68E-2 0.39E-2 0.30E-2 0.21E-2 O.l6E-2 O.l2E-2 O.lOE-2

O.l9E-1 O.l3E-1 0.76E-2 0.57E-2 0.40E-2 0.32E-2 0.24E-2 0.20E-2

203

N. G. Zamani 6

TABLE

h

LSFE

FD

FE

l/3 l/4 l/5 l/6 l/7 l/8 l/9 l/10

0.68E-2 0.42E-3 O.l7E-3 O.l2E-3 0.74E-4 0.66E-4 0.98E-3 0.65E-3

0.46E-2 O.l9E-2 O.l3E-2 O.lOE-2 0.71E-3 0.55E-3 0.45E-3 0.35E-3

0.33E-2 0.29E-2 0.27E-2 0.20E-2 0.22E-2 0.23E-2 0.21E-2 0.20E-2

where 0 is first taken to be the square [0, l]x[O, l] and then taken to be a right isosceles triangle with vertices at (0, 0), (1,0) and (0,l). The results for square region are tabulated in Table 5 and Fig. 4c, whereas the results for triangular region are recorded in Table 6 and Fig. 4d. Notice that the exact solution is U(X,y) = (x*-x) sin my. Example 4. Finally we consider the same differential equation as in the previous example but now defined on an irregular region. fl is now bounded by the x-axis and the parametric polar curve r=cos38

O~f3<1r/6

(see Fig. 6).

The exact solution is as in Example 3. The results for LSFE are recorded in Table 7. The exotic behavior of error can partially be blamed on the ill-conditioning of the normal equations. This topic will be discussed in Section IV. EfFciency and accuracy There is no doubt in our mind that from efficiency point of view, the least squares finite element method is not comparable to the classical finite element technique. A quick look at the size of the matrices involved indicates that the amount of storage location for ISFE is larger than FD and FE. Another disadvantage of least squares is the necessity of using basis functions with high degree of smoothness. Of course by modifying the definition of inner product in Section II one can use nonconforming elements with lower degree of smoothness. TABLE 7

204

h

LSFE

lb l/7 l/8 l/9 l/10

0.27E-4 0.36E-4 O.l3E-3 0.23E-4 O.l8E-4

Journal of The Franklin Instihlte Pergamon Press Ltd.

A Least Squares Finite Element Method Applied to B-Splines TABLE

h

8

llu- %4llk,,

l/6 l/7 l/8 l/9 l/10

0.46E-7 O.l4E-7 O.llE-5 O.l7E-5 O.l2E-5

We are convinced that as far as accuracy is concerned, the least squares method is comparable to finite differences and classical finite elements. This can be observed in Fig. 4. From this diagram it can also be deduced that as the mesh spacing is refined, the error decreases until the round off takes over and the error deteriorates. The error is particularly good for large mesh spacing. Smoothness of the solution In the statement of Theorem I, it is assumed that the solution belongs to C”(n). This severe smoothness assumption can be ignored for computational considerations. It is worth remembering that the solution to Example 1, has logarithmic singularities at the corners of the unit square. Some practical remarks Obviously if the true solution to our problem is not known, we cannot evaluate the error explicitly. However it is possible to use the approximate solution in deriving an a postetioti upper bound for the error. To be more precise, we refer to the inequality:

14x, Y)-u&,

~)l~Wx,

Y>

which holds for (x, y) E 0’. Here u is the true solution to Au=f u=g

Q

ai2

and u, is our least squares finite element approximation. The size of the residual can be computed from the quantity in brackets provided we have some knowledge of K(x, y). Unfortunately, in most cases the value of K(x, y) cannot be computed explicitly. The quantity j\u - uM(J&, which is closely related to the terms in the above brackets is calculated for Example 4 and results are tabulated in Table 8. IV. Ill-conditionality When implementing LSFE numerically we occasionally run into an illconditioned set of normal equations. The source of this ill-conditioning is explained below. Vol. printed

3, in Northern

Ireland

205

N. G. Zamani

FIG. 5.

Consider Fig. 5, pi is the B-spline associated with the dotted node. The intersection of support of Bi with the region R has a very small area which will result in small numbers in the ith column of Grammian matrix. r-_________((ai,~,))__________ ---_______((yei, $$$J)__________ G = [((o&, By))]=

---------- ((ai, s))---------__________((98,48,))_________

Tith column Such a phenomenon is reflected in Example 4. More specifically, we have tabulated the condition numbers of the Grammian for different h, recorded in Table 9. Reducing

the condition number

There are certain measures one can take to reduce the condition number. To explain the idea, we return to Example 4. Au = 2 sin ry - m*(x* - x) sin fry u = (x*-x)

sin 7ry,

fi an

where R is half the air foil. We take p(h) = l/h4 and h =$. The grid system which appears in Fig. 6 is part of the grid which was constructed on the embedding square. The rest of the grid points have no contribution to our algorithm. TABLE 9

206

h

Condition number

l/6 l/7 l/8 l/9 l/10

0.28E16 0.35E19 0.81E15 0.31E14 0.37E13

A Least Squares

Finite Element

Method Applied to B-Splines

FIG. 6.

All

0.26E-6

0.32E-6

0.32E-6

functions ore wesent:

045E-7

0.23E-5

O.S8E-6

0.20E-5

0.2SE-5

0.45E-7

0.23E-5

0.87E6

Q2OE-5

0.26E-5

0.69E-6

0.20E-5

0.29E-5

0.73E-7

Q94E-7

O.l5E-3

O.l3E-3

024E-3

0.32E-5

0.39E-6

Cl70E-7

0.23E-5

CXSE-7

0.23E-6

k 0.35E-6

O.l7E-5

0.4 I E-4

FIG. 7. Vol. 311, No. 3. March 1981 Printed in Northern Ireland

207

N. G. Zamani Now we observe that the B-spline associated with grid point 8 has a small contribution to the approximation. The Grammian has a condition number of 0.35E19. If we delete gBs from the approximating basis i.e. d* i=l,i#8

the resulting Grammian has a condition number 0.37E13 which is considerably smaller than the one above. Figure 7 gives the distribution of error when different B-splines are being deleted from the approximating basis. V. Conclusion The proper choice of B-splines as elements of 9 enables us to solve a variety of boundary value problems. The existence of a kernel function with proper integrability conditions immediately gives pointwise convergence results. An important problem which meets these requirements is the Biharmonic equation. In the absence of a Green’s function one can utilize a ptioti inequalities to prove L, convergence. This is the situation that arises in solving complicated elliptic-parabolic partial differential equations (10). Finally, this method with slight modification of 9 can be applied to Navier’s equation, a coupled system of elliptic equations which arises in three dimensional elasticity (2). Acknowledgement The author is grateful for the support of Professors J. H. Ahlberg and G. Leibbrandt. This research was supported in part by the Natural Sciences and Engineering Research Council of Canada under grant No. A8063. References (1)J. H. Ahlberg,, E. N. Nilson and J. L. Walsh, “The Theory of Splines and Their Applications”, Academic Press, New York, 1967. (2) S. Bergman and M. Schiffer, “Kernel Functions and Elliptic Differential Equations in Mathematical Physics”, Academic Press, New York, 1953. (3) J. H. Bramble and A. H. Schatz, “Least squares method for 2mth order elliptic boundary value problems”, Maths Cornput., Vol. 25, pp. l-32, 1971. (4) E. T. Copson, “Partial Differential Equations”, Cambridge University Press, Cambridge 1975. (5) P. J. Davis and P. Rabinowitz, “Advances in Orthonormalizing Computation”, Advances in Computers, Vol. II, (Edited by F. J. Ah) New York, 1961. (6) S. G. Mikhlin, “Varational Methods in Mathematical Physics”, Macmillan, New York, 1964. (7) J. T. Oden and J. N. Reddy, “Mathematical Theory of Finite Elements”, John Wiley and Sons, New York, 1976. (8) P. M. Prenter, “Splines and Variational Methods”, John Wiley and Sons, New York, 1975. (9) S. M. Serbin. “Computational investigations of least squares type methods for the approximate solution of boundary value”, Maths Comput., Vol. 29, pp. 777793, 1975. (10)V. G. Sigillito, “Explicit a ptiori Inequalities with Applications to Boundary Value Problems”, Pitman Publishing, London, 1977.

208

Journal of ‘Ihe Franklin Institute Pergamon Press Ltd.