Bernatz R. Fourier Series and Numerical Methods for Partial Differential Equations

Подождите немного. Документ загружается.

252 FINITE ELEMENT METHOD

Section 10.2 illuminates the aspects of the finite element method, introduced in

this section, for the case of a ID elliptic boundary value problem. In addition, the

equivalence (meaning identical solutions sets) of the original BVP and the reformu-

lations (10.3) and (10.7) is addressed in Section 10.2.2 for the case of the simple

operator.

10.2 1D ELLIPTICAL EXAMPLE

Finite element methods are applied to a simple ID, second-order BVP in this section.

The simple nature of the problem provides a clear mathematical and physical under-

standing of the various formulation relationships. The equivalence in reformulations

is established next. Finally, the details of the Galerkin method are given in a careful

manner so that it may be easier to understand similar steps for higher dimension

cases.

The ID BVP to be considered is

f -u"{x) = f(x) (PDE)

BVP<^ (10.9)

[ u(a) = u(b) = 0 (BC)

A BVP of this form arises in the case of

thin, long elastic medium stretched between

fixed points x = a and x = b. The dependent variable u represents the displacement

from the reference point (u

—

0) of the medium at location x. The nonhomogeneous

term / represents a load on the elastic medium acting in direction perpendicular to

the orientation of the medium. An alternate application corresponding to the BVP

is the case of a long, slender metallic rod where u represents the temperature of the

medium at x in response to a heat source given by f(x). Of course, the temperature

at the ends of the rod (x = a and x = b) is zero. The heat conductivity coefficient

for this simple case is unity.

Note: The linear operator in this example is D =

—

¿2~.

A unique solution u

from V, the set of all twice-differentiable functions on [a,b] with boundary values

u(a) = u(b) = 0, for BVP (10.9) exists if the function / is piecewise continuous

and bounded on the interval

[a,b].

The finite element process begins by developing

reformulations of BVP (10.9) as shown in Section

10.2.1.

10.2.1 Reformulations

The functional

of the variational form of BVP (10.9) is given by

I(v)=

-(Av,v)-(f,v)

where the scalar product (·, ·) is defined to be

(v,u) — \ v(x)u(x)da

1D ELLIPTICAL EXAMPLE 253

Using this definition of the scalar product, the method of integration by parts, and the

given boundary values of v G V, the functional / has an equivalent form

I(v) = \{v',v')-{v,f)

If the BVP (10.9) applies to the elastic medium with load intensity /, then the

functional / gives a measure of the total potential energy of the medium where

\(v

, v') represents the internal elastic energy and (/, v) is an expression of the load

potential.

The objective of the minimization problem is to find u eV such that

I(u)<I(v) VveV (10.10)

In the context of the elastic medium problem, it is evident that the solution to the

PDE is the stationary displacement u(x) that minimizes the total potential energy.

The weak formulation of the BVP (10.9) leading to the Galerkin approximation is

found by first multiplying both sides of the ODE by a test function v and integrating

over the interval [a,

b].

That is,

-u"(x) = f(x) => -u"(x)v(x) = f(x)v(x)

pb pb

=> / — u"{x)v(x) = / f(x)v(x)dx

Ja Ja

(-u",v)

= (f,v)

The weak solution u

V is such that

(Du,v) = (f,v) VveV (10.11)

For sufficiently continuous u"\ integration by parts and the fact that v(a) = v(b) = 0

may be used to give

pb pb

(Du,v) = / u"{x)v{x)dx = / u'{x)v'{x)dx = (u',v

)

Ja Ja

so that the weak formulation is equivalent to finding u G V such that

(

\v')

= (f,v) WveV (10.12)

The equivalence of the BVP

(10.9),

the variational problem expressed in Inequality

(10.10),

and the weak formulation in Equation (10.11) is established in Section 10.2.2

prior to providing details for the finite element solution to the BVP.

10.2.2 Equivalence in Forms

Section 10.2.1 presented a simple ID, second-order BVP and two reformulations of

the BVP: the variational form given by Inequality (10.10), and the weak form given

254 FINITE ELEMENT METHOD

by Equation (10.11). The objective in this section is to show that, under certain

conditions on u", the BVP and its two reformulations are equivalent in that they

admit the same unique solution u.

The simple nature of

the

BVP allows that, indeed, a unique solution u exists to the

BVP provided the function / is sufficiently continuous and bounded. By definition

of the weak formulation, if u solves the differential equation, it follows that u solves

the weak formulation. For the converse, suppose u solves the weak formulation so

that, for every v G V it follows:

(Du,v) = (f,v) =»

(-u",v)

= (f,v)

rb pb

=>*

/ —u"{x)v(x)dx

—

I f(x)v(x)dx

Ja Ja

rb pb

I — u" {x)v{x)dx

—

I f(x)v(x)dx = 0

Ja Ja

(-u"(x)

- f(x))v(x)dx = 0

Because v is arbitrary, if -u"(x) —f(x) is continuous on [a,

b],

the last equality above

is true only if — u"(x) = f(x). Consequently, if a solution to the weak formulation

has a continuous second derivative, then it solves the original BVP (10.9).

Equivalence between the variational problem and the weak formulation is estab-

lished next to complete the object of the current section. Suppose u is a solution to

the minimization problem so that I(u) < I(v) for all v G V. Let v be an arbitrary

element of V. The linear nature of V ensures that u + av G V for any real number

a. Define the real-valued function g as g (a) = I(u + av). It follows that g has a

minimum for a = 0. Now,

g (a) = I{u + av)

= -{(u + av)', (u + av)

) - (/, (u + av))

= ^(u

)+a(v',u

) + \a

{v',v') - (f,u)-a(f,v)

and

g'(a)

= (v',u')+a(v',v')-(f,v)

Because g has a minimum at a = 0, it follows that g'(0) = (v

) + 0(v

)

—

(/, υ) = 0 => (ν',η') = (/, v), which establishes the solution u to the variational

problem solves the weak formulation of the PDE.

Let u satisfy the weak formulation so that (?/, v!) = (/, v) is true for all v G V.

To show that u minimizes the function 7, let v be an arbitrary element of V, and

let w = v

—

u. The linearity of V ensures tí; is a member so that it follows by u

satisfying the weak formulation that

<«',

w') = {f, w) =>

<«',

υ') - («', u'> = (/, υ) - (/, u)

1D ELLIPTICAL EXAMPLE 255

=> <u

)-</,w> = (u

)-(/,v>

=> I(u) = -\(u',u

) + <u',t/> - ±<t/,t;'> +!(*)

=> J(u) = -±(ti;>'>+J(t;)

=» J(u) < I{y)

where the last implication follows because (it/, w') > 0.

10.2.3 Finite Element Solution

The details a finite element solution for a ID BVP will be described in this section.

The specific BVP is

( -u"(x) = 10x(l - x) (DE)

BVP^

(10.13)

[ u(0) = u(l) = 0 (BC)

The solution u is the unique function from the set V of functions with continuous

second derivatives have boundary values u(0) = w(l) = 0. An approximation u

u from a subspace V

of V will be determined by solving the weak formulation on

the set V

. This constitutes the Galerkin method for finite elements.

To establish a finite-dimensional subspace of V, the interval [0,1] is partitioned

into subintervals of uniform length. Each subinterval is an "element." For this

example, 10 uniform elements are used with the endpoints given by 0, x\ =0.1,

X2 = 0.2, ..., £9 = 0.9, and xio = 10. The values x¿ are the "nodes" at which the

values for the approximation u

is calculated. The nodal locations are depicted as

open circles on the x-axis in Figure 10.1. The subspace V

is the set of piecewise

linear functions on [0,1] with linear segments corresponding to the finite elements of

[0,1].

Any piecewise linear function on [0,1] is uniquely determined by its values

at the nodes x¿. Therefore, a basis for V

is given by the set of piecewise linear

functions φι{χ) where φι{χι) - 1 and

4>Í{XJ)

= 0 for j φ i. Plots of the basis

functions φ\ and φ^ are shown in Figure 10.1. Consequently, any v

G V

is such

that

= ^2^i(x)

i=l

The Galerkin solution u

is also expressed as a linear combination of the basis

functions with

3 = 1

Using the weak formulation given in Equation (10.12), the Galerkin solution is

the function u

G V

such that

(u'

) = (/, v

)

e V

256 FINITE ELEMENT METHOD

0.1

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure

10.1

Nodes and basis functions φ\{χ) and

<t>2(x).

Because functions φι form a basis of the space

this is equivalent to

(0;.(x),^(x))

(/(x),^(x)) ¿

1,2,...,

which constitutes nine equations

the nine unknowns Uj. The resulting system

equations in matrix-vector form is

Au =

where

and

all &12 · · · Û19 '

«21

&22 · · · «29

agi α

... α

<*ij

= (Φ'ϋΦ'ο) =

Ûg

Φ)

fti

(10.14)

Ax)dx

= (/,&) = /

¡(χ)φί(χ)άα

With the given definitions

the basis functions

φ^ it

follows that the diagonal

entries an

20, the off-diagonal entries

a^j =

—10 for \i

—

1, and

α^ =

0 for

all others (see Exercise 10.3). For the entries in

= I 10x(l

—

χ)φι(χ)άχ

10x(l

χ)φ

(χ)άχ

=—-

(10.15)

600

and so on for the other eight values.

Once all the entries in matrix

and vector

have been determined, the system of

nine equations is easily solved because of the tridiagonal structure of

In this case,

solving for

gives

327

58 847 31 25 31 847 58

327

4000'

375' 4000' Ϊ25' 96' 125' 4000' 375' 4000

)

2D ELLIPTICAL EXAMPLE 257

The approximated values for u(x) are plotted as open circles in Figure 10.2 and the

exact analytic solution is graphed as a solid line. The true solution to BVP (10.13)

0.25

Figure 10.2 True solution u(x) (solid line) and the finite element approximation (open

circles).

is easily determined to be u{x) = |x

—

+ |x. Calculating u at the nine finite

element nodes gives

327 58 847 31 25 31 847 58 327 \

4000'

375' 4000' Ï25' 95' 125' 4000' 375' 4000 J

the same values given by the finite element method.

10.3 2D ELLIPTICAL EXAMPLE

The finite element method based on the weak formulation for approximating the

solution to a 2D Poisson problem is presented in this section. The BVP is

BVP

(10.16)

I u(x,y)\dn=0 (BC)

The weak formulation of BVP

( 10.16)

is presented next.

10.3.1 Weak Formulation

The weak formulation to the BVP is found by multiplying the PDE in BVP (10.16)

by a test function v(x,y) and the integrating over the domain Ω to give

where

(Du,v) = (f,v)

(Du,v) = / / —V

u(x,y)v(x,y)dxdy

(10.17)

258 FINITE ELEMENT METHOD

and

(fi

)= / / f{x,y)v(x,y)dxdy

/Ω

The function v is a member of the set V of all functions with continuous second

partíais on the interior of

and having zero value of the boundary 9Ω of V. Recall

that it is necessary for the solution to the BVP (10.16) to be continuous up to at least

its second

partíais.

Using Green's Theorem and the fact that v is zero on the boundary

of Ω, it follows that

(Du,v) — I I Vu(x,y) -\7v(x,y) dxdy

J Jn

so that the weak formulation of the BVP

( 10.16)

is to find u G V such that

/ / Vu{x,y) -\7v(x,y) dxdy = / / f(x,y)v(x,y)dxdy (10.18)

for all v eV.

10.3.2 Finite Element Approximation

For the sake of developing the finite element approximation to u, the domain Ω is the

unit square. That is,

Ω = {(x, y)\0 < x < 1 and 0 < y < 1}

The approximation

will be determined at 16 evenly spaced locations with Ω, as

shown in Figure 10.3. The domain Ω is partitioned into triangular elements with the

vertices of the triangles coinciding with the nodes. The nodes are numbered 1-16

and the triangular elements are numbered 1-50. The finite-dimensional subspace Vh,

from which the approximate solution

will be found, is the set of all continuous

functions that are linear over each of the triangular elements shown in Figure 10.3.

A set of basis φι(χ, y) functions for Vh are function having value one at node i and

zero at all other nodes. The functions are constructed by defining functions ψ^(χ, y)

on the triangles having node i as a vertex. Each function ipj(x, y) is linear and of the

form ax + by + c. Figure 10.4 shows the graph of basis function φ\(χ, y) and the

functions

i/)j(x,

y). From the figure, one can see that

Φι

(x,

V>i(x,

ψ2(χ,

ψζ(χ,

^14

(z,

^13

(^,

^12(2,

Consequently, the set {φι, φ^,

··.,

φ\§) is a basis for the set Vh of piecewise-linear

functions on Ω, with value zero on the boundary <9Ω.

The Galerkin solution to Equation (10.17) from the set Vh is the function

J2j

ßj4>j

e V

, such that

(Ç^i) · ^(Y^a^dxdy = J J fY^a^idxdy (10.19)

2D ELLIPTICAL EXAMPLE 259

Figure 10.3 Nodes and triangular elements for the 2D

finite

element approximation.

Figure 10.4 Basis function φ\(χ, y) construction using ^j(x, y), j

14,

13,

and 12.

for all v = Σί

ιΦί ^ ^· Because functions φι form a basis for V^, this is equivalent

to finding ßj such that

ν(

βόΦό)

' ^Φΐ dxdy = / / /Φί dxdy i =

1,2,3,...,

from which follows

/ / ΣΡόΦΦΐ

^) dxdy = [ [

f</>i

dxdy

J Jn ■ J Jn

i =

1,2,3,...,

16 (10.20)

260 FINITE ELEMENT METHOD

which generates a system of 16 equations in the 16 unknowns ßj, Mß = b. The

matrix M is called the stiffness matrix with entries given by

Mij = ^Φο · V0i dxdy

and the values for constant vector

given by

= f Í

ΪΦι

Because the basis functions φι are zero on much of the domain Ω, the matrix M is

sparse. In fact, the nature of the basis functions is such that there is only nonzero

"interaction" with the basis function for a given node and the four nearest-neighbor

nodes.

Using the details for determining the entries in M established in Exercise

10.11,

the "interaction" pattern and corresponding "weights" that result are shown

in Figure 10.5. Additionally, the commutativity of the inner product results in a

symmetric stiffness matrix as w(

i 0

i-10

j-2

i-2

Figure 10.5 Stiffness matrix pattern based on near nodes

0 ... 0

1 ... 0

0 ... 0

The

"-1"

term in the first row, fifth column (identified by the square) corresponds to

the case with i

—

1, where the near node in the j +

location is node number

Also,

note the "0" entry in row 4 and column 5 (identified by the square) results because

there is no common support (no nonzero "interaction") for nodes number 4 and 5.

©©00

0 0 0 0

i- 1 i i + l i + 2

The stiffness matrix for this case is

M :

-1

ERROR ANALYSIS 261

To complete the example, suppose the function / in BVP (10.16) is given as

f(x, y) = 2x

+ 2y

—

2y, and the domain Ω is the unit square {(#, y)\0 <

x < 1 and 0 < y < 1}. For the sake of determining the error in the finite element

approximation, the true solution to BVP (10.16) is u(x, y) — (x

—

x)(y

—

y).

The stiffness matrix M has already been determined above. What must be deter-

mined is the right-hand-side vector b with entry i defined as bi = J J

(2x

+ 2y

—

— 2y)4>i

dxdy. Once b is calculated, the system of equations is solved using

Gauss-Jordan elimination, which is a reasonable linear system technique given the

small size of the system.

The finite element solution and the true solution to this BVP are plotted in Figure

10.6.

The true solution surface is plotted with the contours and the finite element

solution values are plotted in the circular symbols at the computational nodes. The

plot indicates the finite element solution provides a "good" approximation to the true

solution. Calculation of the error at each of the nodes indicates the maximum error

of magnitude -^ ~ 0.0017777778 at the nodes (0.4,0.4), (0.6,0.4), (0.4,0.6), and

(0.6,0.6).

Figure 10.6

True

solution (solid

contours) and the

finite

element solution (circular

symbols).

The sparse nature of the stiffness matrix M makes Gauss-Jordan elimination

inefficient for large systems. For larger systems, direct methods, such as Cholesky's,

where the matrix M is factored into a product of lower and upper triangular matrices,

are better means for solving the linear system. An alternative to direct methods for

symmetric matrices involves solving and equivalent quadratic minimization problem

given by \

· My

—

b · y\. Here, various iterative methods for finding the unique

minimizing vector include the gradient method, as well as the conjugate gradient

method.

10.4 ERROR ANALYSIS

A discussion of

the

error associated with the finite element is provided in this section.

The stated results are not proved. The interested reader should consult books such