Introduction

Back

___ Historical Perspective

The _nite element method is a computational technique for obtaining approximate solu_

tions to the partial di_erential equations that arise in scienti_c and engineering applica_

tions_ Rather than approximating the partial di_erential equation directly as with_ e_g__

_nite di_erence methods_ the _nite element method utilizes a variational problem that

involves an integral of the di_erential equation over the problem domain_ This domain

is divided into a number of subdomains called _nite elements and the solution of the

partial di_erential equation is approximated by a simpler polynomial function on each

element_ These polynomials have to be pieced together so that the approximate solution

has an appropriate degree of smoothness over the entire domain_ Once this has been

done_ the variational integral is evaluated as a sum of contributions from each _nite el_

ement_ The result is an algebraic system for the approximate solution having a _nite

size rather than the original in_nite_dimensional partial di_erential equation_ Thus_ like

_nite di_erence methods_ the _nite element process has discretized the partial di_eren_

tial equation but_ unlike _nite di_erence methods_ the approximate solution is known

throughout the domain as a pieceise polynomial function and not just at a set of points_

Logan ___       attributes the discovery of the _nite element method to Hrennikof _     and

McHenry ___  who decomposed a two_dimensional problem domain into an assembly of

one_dimensional bars and beams_ In a paper that was not recognized for several years_

Courant __      used a variational formulation to describe a partial di_erential equation with

a piecewise linear polynomial approximation of the solution relative to a decomposition of

the problem domain into triangular elements to solve equilibrium and vibration problems_

This is essentially the modern _nite element method and represents the _rst application

where the elements were pieces of a continuum rather than structural members_

Turner et al_ ___         wrote a seminal paper on the subject that is widely regarded

_

 Introduction

as the beginning of the _nite element era_ They showed how to solve one_ and two_

dimensional problems using actual structural elements and triangular_ and rectangular_

element decompositions of a continuum_ Their timing was better than Courant_s __    _

since success of the _nite element method is dependent on digital computation which

was emerging in the late ____s_ The concept was extended to more complex problems

such as plate and shell deformation _cf_ the historical discussion in Logan ___ _ Chapter

__ and it has now become one of the most important numerical techniques for solving

partial di_erential equations_ It has a number of advantages relative to other methods_

including

_ the treatment of problems on complex irregular regions_

_ the use of nonuniform meshes to re_ect solution gradations_

_ the treatment of boundary conditions involving _uxes_ and

_ the construction of high_order approximations_

Originally used for steady _elliptic_ problems_ the _nite element method is now used

to solve transient parabolic and hyperbolic problems_ Estimates of discretization errors

may be obtained for reasonable costs_ These are being used to verify the accuracy of the

computation_ and also to control an adaptive process whereby meshes are automatically

re_ned and coarsened and_or the degrees of polynomial approximations are varied so as

to compute solutions to desired accuracies in an optimal fashion ___ _ __ __ __ __ __           _

___ Weighted Residual Methods

Our goal_ in this introductory chapter_ is to introduce the basic principles and tools of

the _nite element method using a linear two_point boundary value problem of the form

L_u      __ _

d

dx

_p_x_

du

dx

_ _ q_x_u _ f_x__ _ _ x _ __ _____a_

u___ _ u___ _ __ _____b_

The _nite element method is primarily used to address partial di_erential equations and is

hardly used for two_point boundary value problems_ By focusing on this problem_ we hope

to introduce the fundamental concepts without the geometric complexities encountered

in two and three dimensions_

Problems like ______ arise in many situations including the longitudinal deformation

of an elastic rod_ steady heat conduction_ and the transverse de_ection of a supported

____ Weighted Residual Methods _

cable_ In the latter case_ for example_ u_x_ represents the lateral de_ection at position

x of a cable having _scaled_ unit length that is subjected to a tensile force p_ loaded by

a transverse force per unit length f_x__ and supported by a series of springs with elastic

modulus q _Figure ______ The situation resembles the cable of a suspension bridge_ The

tensile force p is independent of x for the assumed small deformations of this model_ but

the applied loading and spring moduli could vary with position_

__

__

__

__

__

__

__

__

q(x) u(x)

p p x

f(x)

Figure _____ De_ection u of a cable under tension p_ loaded by a force f per unit length_

and supported by springs having elastic modulus q_

Mathematically_ we will assume that p_x_ is positive and continuously di_erentiable

for x _ ___ _    _ q_x_ is non_negative and continuous on ___ _         _ and f_x_ is continuous on

___ _   _

Even problems of this simplicity cannot generally be solved in terms of known func_

tions_ thus_ the _rst topic on our agenda will be the development of a means of calculating

approximate solutions of _______ With _nite di_erence techniques_ derivatives in _____a_

are approximated by _nite di_erences with respect to a mesh introduced on ___ _      __       _

With the _nite element method_ the method of weighted residuals _MWR_ is used to

construct an integral formulation of ______ called a variational problem_ To this end_ let

us multiply _____a_ by a test or weight function v and integrate over ___ __ to obtain

_v_L_u            _ f_ _ __ ____a_

We have introduced the L_ inner product

_v_ u_ __ Z _

_

vudx ____b_

to represent the integral of a product of two functions_

The solution of ______ is also a solution of ____a_ for all functions v for which the

inner product exists_ We_ll express this requirement by writing v _ L____ ___ All functions

of class L____ __ are _square integrable_ on ___ ___ thus_ _v_ v_ exists_ With this viewpoint

and notation_ we write ____a_ more precisely as

_v_L_u            _ f_ _ __ _v _ L____ ___ ____c_

_ Introduction

Equation ____c_ is referred to as a variational form of problem _______ The reason for

this terminology will become clearer as we develop the topic_

Using the method of weighted residuals_ we construct approximate solutions by re_

placing u and v by simpler functions U and V and solving ____c_ relative to these

choices_ Speci_cally_ we_ll consider approximations of the form

u_x_ _ U_x_ _

N X

j__

cj_j_x__ _____a_

v_x_ _ V _x_ _

N X

j__

dj_j_x__ _____b_

The functions _j_x_ and _j_x__ j _ __ _ _ _ _ _N_ are preselected and our goal is to

determine the coe_cients cj _ j _ __ _ _ _ _ _N_ so that U is a good approximation of u_

For example_ we might select

_j_x_ _ _j_x_ _ sin j_x_ j _ __ _ _ _ _ _N_

to obtain approximations in the form of discrete Fourier series_ In this case_ every function

satis_es the boundary conditions _____b__ which seems like a good idea_

The approximation U is called a trial function and_ as noted_ V is called a test func_

tion_ Since the di_erential operator L_u           is second order_ we might expect u _ C____ ___

_Actually_ u can be slightly less smooth_ but C_ will su_ce for the present discussion__

Thus_ it_s natural to expect U to also be an element of C____ ___ Mathematically_ we re_

gard U as belonging to a _nite_dimensional function space that is a subspace of C____ ___

We express this condition by writing U _ SN___ __ _ C____ ___ _The restriction of these

functions to the interval _ _ x _ _ will_ henceforth_ be understood and we will no longer

write the ___ ____ With this interpretation_ we_ll call SN the trial space and regard the

preselected functions _j_x__ j _ __ _ _ _ _ _N_ as forming a basis for SN_

Likewise_ since v _ L__ we_ll regard V as belonging to another _nite_dimensional

function space _ SN called the test space_ Thus_ V _ _ SN _ L_ and _j_x__ j _ __ _ _ _ _ _N_

provide a basis for _ SN_

Now_ replacing v and u in ____c_ by their approximations V and U_ we have

_V_L_U          _ f_ _ __ _V _ _ SN_ _____a_

The residual

r_x_ __ L_U    _ f_x_ _____b_

____ Weighted Residual Methods _

is apparent and clari_es the name _method of weighted residuals__ The vanishing of the

inner product _____a_ implies that the residual is orthogonal in L_ to all functions V in

the test space _ SN_

Substituting ______ into _____a_ and interchanging the sum and integral yields

N X

j__

dj__j _L_U      _ f_ _ __ _dj _ j _ __ _ _ _ _ _N_ ______

Having selected the basis _j _ j _ __ _ _ _ _ _N_ the requirement that _____a_ be satis_ed for

all V _ _ SN implies that ______ be satis_ed for all possible choices of dk_ k _ __ _ _ _ _ _N_

This_ in turn_ implies that

__j _L_U         _ f_ _ __ j _ __ _ _ _ _ _N_ ______

Shortly_ by example_ we shall see that ______ represents a linear algebraic system for the

unknown coe_cients ck_ k _ __ _ _ _ _ _N_

One obvious choice is to select the test space _ SN to be the same as the trial space

and use the same basis for each_ thus_ _k_x_ _ _k_x__ k _ __ _ _ _ _ _N_ This choice leads

to Galerkin_s method

__j_L_u           _ f_ _ __ j _ __ _ _ _ _ _N_ ______

which_ in a slightly di_erent form_ will be our _work horse__ With _j _ C__ j _

__ _ _ _ _ _N_ the test space clearly has more continuity than necessary_ Integrals like

______ or ______ exist for some pretty _wild_ choices of V _ Valid methods exist when V

is a Dirac delta function _although such functions are not elements of L__ and when V

is a piecewise constant function _cf_ Problems _ and  at the end of this section__

There are many reasons to prefer a more symmetric variational form of ______ than

______ e_g__ problem ______ is symmetric _self_adjoint_ and the variational form should

re_ect this_ Additionally_ we might want to choose the same trial and test spaces_ as with

Galerkin_s method_ but ask for less continuity on the trial space SN_ This is typically

the case_ As we shall see_ it will be di_cult to construct continuously di_erentiable

approximations of _nite element type in two and three dimensions_ We can construct

the symmetric variational form that we need by integrating the second derivative terms

in ____a_ by parts_ thus_ using _____a_

Z _

_

v___pu___ _ qu _ f     dx _ Z _

_

_v_pu_ _ vqu _ vf_dx _ vpu_j_

_ _ _ _____

where _ __ _ d_ __dx_ The treatment of the last _boundary_ term will need greater

attention_ For the moment_ let v satisfy the same trivial boundary conditions _____b_ as

_ Introduction

u_ In this case_ the boundary term vanishes and _____ becomes

A_v_ u_ _ _v_ f_ _ _ _____a_

where

A_v_ u_ _ Z _

_

_v_pu_ _ vqu_dx_ _____b_

The integration by parts has eliminated second derivative terms from the formulation_

Thus_ solutions of ______ might have less continuity than those satisfying either ______ or

______ For this reason_ they are called weak solutions in contrast to the strong solutions

of ______ or ______ Weak solutions may lack the continuity to be strong solutions_ but

strong solutions are always weak solutions_ In situations where weak and strong solutions

di_er_ the weak solution is often the one of physical interest_

Since we_ve added a derivative to v by the integration by parts_ v must be restricted

to a space where functions have more continuity than those in L__ Having symmetry in

mind_ we will select functions u and v that produce bounded values of

A_u_ u_ _ Z _

_

_p_u___ _ qu_            dx_

Actually_ since p and q are smooth functions_ it su_ces for u and v to have bounded

values of

Z _

_

__u___ _ u_    dx_ _______

Functions where _______ exists are said to be elements of the Sobolev space H__ We_ve

also required that u and v satisfy the boundary conditions _____b__ We identify those

functions in H_ that also satisfy _____b_ as being elements of H_

_ _ Thus_ in summary_

the variational problem consists of determining u _ H_

_ such that

A_v_ u_ _ _v_ f__ _v _ H_

_ _ _______

The bilinear form A_v_ u_ is called the strain energy_ In mechanical systems it frequently

corresponds to the stored or internal energy in the system_

We obtain approximate solutions of _______ in the manner described earlier for the

more general method of weighted residuals_ Thus_ we replace u and v by their approxi_

mations U and V according to _______ Both U and V are regarded as belonging to the

same _nite_dimensional subspace SN

_ of H_

_ and _j _ j _ __ _ _ _ _ _N_ forms a basis for

SN

_ _ Thus_ U is determined as the solution of

A_V_U_ _ _V_ f__ _V _ SN

_ _ _____a_

____ Weighted Residual Methods _

The substitution of _____b_ with _j replaced by _j in _____a_ again reveals the more

explicit form

A__j_ U_ _ __j_ f__ j _ __ _ _ _ _ _N_ _____b_

Finally_ to make _____b_ totally explicit_ we eliminate U using _____a_ and interchange

a sum and integral to obtain

N X

k__

ckA__j_ _k_ _ __j_ f__ j _ __ _ _ _ _ _N_ _____c_

Thus_ the coe_cients ck_ k _ __ _ _ _ _ _N_ of the approximate solution _____a_ are deter_

mined as the solution of the linear algebraic equation _____c__ Di_erent choices of the

basis _j_ j _ __ _ _ _ _ _N_ will make the integrals involved in the strain energy _____b_

and L_ inner product ____b_ easy or di_cult to evaluate_ They also a_ect the accuracy

of the approximate solution_ An example using a _nite element basis is presented in the

next section_

Problems

__ Consider the variational form ______ and select

_j_x_ _ __x _ xj__ j _ __ _ _ _ _ _N_

where __x_ is the Dirac delta function satisfying

__x_ _ __ x __ __ Z _

__

__x_dx _ __

and

_ _ x_ _ x_ _ __ _ _ xN _ __

Show that this choice of test function leads to the collocation method

L_U     _ f_x_jx_xj _ __ j _ __ _ _ _ _ _N_

Thus_ the di_erential equation ______ is satis_ed exactly at N distinct points on

___ ___

_ The subdomain method uses piecewise continuous test functions having the basis

_j_x_ __ _ __ if x _ _xj_____ xj_____

__ otherwise

_

where xj____ _ _xj _ xj_____ Using _______ show that the approximate solution

U_x_ satis_es the di_erential equation _____a_ on the average on each subinterval

_xj_____ xj______ j _ __ _ _ _ _ _N_

 Introduction

__ Consider the two_point boundary value problem

_u__ _ u _ x_ _ _ x _ __ u___ _ u___ _ __

which has the exact solution

u_x_ _ x _

sinh x

sinh _

_

Solve this problem using Galerkin_s method _____c_ using the trial function

U_x_ _ c_ sin _x_

Thus_ N _ __ ___x_ _ ___x_ _ sin _x in _______ Calculate the error in strain

energy as A_u_ u_ _ A_U_ U__ where A_u_ v_ is given by _____b__

___ A Simple Finite Element Problem

Finite element methods are weighted residuals methods that use bases of piecewise poly_

nomials having small support_ Thus_ the functions __x_ and __x_ of ______ _____ are

nonzero only on a small portion of problem domain_ Since continuity may be di_cult to

impose_ bases will typically use the minimum continuity necessary to ensure the existence

of integrals and solution accuracy_ The use of piecewise polynomial functions simplify

the evaluation of integrals involved in the L_ inner product and strain energy ____b_

____b_ and help automate the solution process_ Choosing bases with small support leads

to a sparse_ well_conditioned linear algebraic system _____c__ for the solution_

Let us illustrate the _nite element method by solving the two_point boundary value

problem ______ with constant coe_cients_ i_e__

_pu__ _ qu _ f_x__ _ _ x _ __ u___ _ u___ _ __ _______

where p            _ and q _ __ As described in Section ___ we construct a variational form of

______ using Galerkin_s method ________ For this constant_coe_cient problem_ we seek

to determine u _ H_

_ satisfying

A_v_ u_ _ _v_ f__ _v _ H_

_ _ _____a_

where

_v_ u_ _ Z _

_

vudx_ _____b_

A_v_ u_ _ Z _

_

_v_pu_ _ vqu_dx_ _____c_

____ A Simple Finite Element Problem _

With u and v belonging to H_

_ _ we are sure that the integrals _____b_c_ exist and that

the trivial boundary conditions are satis_ed_

We will subsequently show that functions _of one variable_ belonging to H_ must

necessarily be continuous_ Accepting this for the moment_ let us establish the goal of

_nding the simplest continuous piecewise polynomial approximations of u and v_ This

would be a piecewise linear polynomial with respect to a mesh

_ _ x_ _ x_ _ _ _ _ _ xN _ _ _______

introduced on ___ _     _ Each subinterval _xj___ xj__ j _ __ _ _ _ _ _N_ is called a _nite element_

The basis is created from the _hat function_

_j_x_ ___

__

_

x_xj__

xj_xj__

_ if xj__            x _ xj

xj___x

xj___xj

_ ifxj     x _ xj__

__ otherwise

_ ______a_

x x x x

1

x

0 j-1 j j+1

j (x)

N x

Figure ______ One_dimensional _nite element mesh and piecewise linear hat function

_j_x__

As shown in Figure ______ _j_x_ is nonzero only on the two elements containing the

node xj _ It rises and descends linearly on these two elements and has a maximal unit

value at x _ xj _ Indeed_ it vanishes at all nodes but xj _ i_e__

_j_xk_ _ _jk __ _ __ if xk _ xj

__ otherwise

_ ______b_

Using this basis with _______ we consider approximations of the form

U_x_ _

N__

Xj__

cj_j_x__ _______

Let_s examine this result more closely_

__ Introduction

x x x x x

x

0 j-1 j j+1 N

j (x) j-1 (x)

c

c

j

j-1

j+1

c

1

U(x)

Figure _____ Piecewise linear _nite element solution U_x__

__ Since each _j_x_ is a continuous piecewise linear function of x_ their summation

U is also continuous and piecewise linear_ Evaluating U at a node xk of the mesh

using ______b_ yields

U_xk_ _

N__

Xj__

cj_j_xk_ _ ck_

Thus_ the coe_cients ck_ k _ __ _ _ _ _ _N _ __ are the values of U at the interior

nodes of the mesh _Figure ______

_ By selecting the lower and upper summation indices as _ and N__ we have ensured

that _______ satis_es the prescribed boundary conditions

U___ _ U___ _ __

As an alternative_ we could have added basis elements ___x_ and _N_x_ to the

approximation and written the _nite element solution as

U_x_ _

N X

j__

cj_j_x__ _______

Since_ using ______b__ U_x__ _ c_ and U_xN_ _ cN_ the boundary conditions are

satis_ed by requiring c_ _ cN _ __ Thus_ the representations _______ or _______ are

identical_ however_ _______ would be useful with non_trivial boundary data_

__ The restriction of the _nite element solution _______ or _______ to the element

_xj___ xj         is the linear function

U_x_ _ cj___j___x_ _ cj_j_x__ x _ _xj___ xj            _ _______

____ A Simple Finite Element Problem __

since _j__ and _j are the only nonzero basis elements on _xj___ xj       _Figure ______

Using Galerkin_s method in the form _____c__ we have to solve

N__

Xk__

ckA__j_ _k_ _ __j_ f__ j _ __ _ _ _ _ _N _ __ ______

Equation ______ can be evaluated in a straightforward manner by substituting replacing

_k and _j using _______ and evaluating the strain energy and L_ inner product according

to _____b_c__ This development is illustrated in several texts _e_g__ __       _ Section ____

We_ll take a slightly more complex path to the solution in order to focus on the computer

implementation of the _nite element method_ Thus_ write _____a_ as the summation of

contributions from each element

N X

j__

_Aj_V_U_ _ _V_ f_j  _ __ _V _ SN

_ _ ______a_

where

Aj_V_U_ _ ASj

_V_U_ _ AMj

_V_U__ ______b_

ASj

_V_U_ _ Z xj

xj__

pV _U_dx_ ______c_

AMj

_V_U_ _ Z xj

xj__

qV Udx_ ______d_

_V_ f_j _ Z xj

xj__

V fdx_ ______e_

It is customary to divide the strain energy into two parts with ASj

arising from internal

energies and AMj

arising from inertial e_ects or sources of energy_

Matrices are simple data structures to manipulate on a computer_ so let us write the

restriction of U_x_ to _xj___ xj           according to _______ as

U_x_ _ _cj___ cj        _ _j___x_

_j_x_ __ __j___x__ _j_x_      _ cj__

cj __ x _ _xj___ xj      _ _______a_

We can_ likewise_ use _____b_ to write the restriction of the test function V _x_ to _xj___ xj          

in the same form

V _x_ _ _dj___ dj       _ _j___x_

_j_x_ __ __j___x__ _j_x_      _ dj__

dj __ x _ _xj___ xj      _ _______b_

_ Introduction

Our task is to substitute ________ into ______c_e_ and evaluate the integrals_ Let us begin

by di_erentiating _______a_ while using ______a_ to obtain

U__x_ _ _cj___ cj      _ ___hj

__hj _ _ ____hj _ __hj            _ cj__

cj __ x _ _xj___ xj      _ _______a_

where

hj _ xj _ xj___ j _ __ _ _ _ _ _N_ _______b_

Thus_ U__x_ is constant on _xj___ xj              and is given by the _rst divided di_erence

U__x_ _

cj _ cj__

hj

_ x _ _xj___ xj            _

Substituting ________ and a similar expression for V __x_ into ______b_ yields

ASj

_V_U_ _ Z xj

xj__

p_dj___ dj       _ ___hj

__hj _____hj _ __hj    _ cj__

cj _dx

or

ASj

_V_U_ _ _dj___ dj     _Z xj

xj__

p _ __h_

j ___h_

j

___h_

j __h_

j _dx    _ cj__

cj __

The integrand is constant and can be evaluated to yield

ASj

_V_U_ _ _dj___ dj     Kj _ cj__

cj __ Kj _

p

hj _ _ __

__ _ __ _______

The    matrix Kj is called the element sti_ness matrix_ It depends on j through hj _

but would also have such dependence if p varied with x_ The key observation is that

Kj can be evaluated without knowing cj___ cj _ dj___ or dj and this greatly simpli_es the

automation of the _nite element method_

The evaluation of AMj

proceeds similarly by substituting ________ into ______d_ to

obtain

AMj

_V_U_ _ Z xj

xj__

q_dj___ dj       _ _j__

_j ___j___ _j   _ cj__

cj _dx_

With q a constant_ the integrand is a quadratic polynomial in x that may be integrated

exactly _cf_ Problem _ at the end of this section_ to yield

AMj

_V_U_ _ _dj___ dj     Mj  cj__cj __ Mj _

qhj

_ _  _

_  __ ________

whereMj is called the element mass matrix because_ as noted_ it often arises from inertial

loading_

____ A Simple Finite Element Problem __

The _nal integral ______e_ cannot be evaluated exactly for arbitrary functions f_x__

Without examining this matter carefully_ let us approximate it by its linear interpolant

f_x_ _ fj___j___x_ _ fj_j_x__ x _ _xj___ xj   _ ________

where fj __ f_xj__ Substituting ________ and _______b_ into ______e_ and evaluating the

integral yields

_V_ f_j _ Z xj

xj__

_dj___ dj         _ _j__

_j ___j___ _j   _ fj__

fj _dx _ _dj___ dj        lj _______a_

where

lj _

hj

_ _ fj__ _ fj

fj__ _ fj __ _______b_

The vector lj is called the element load vector and is due to the applied loading f_x__

The next step in the process is the substitution of ________ _________ and ________ into

______a_ and the summation over the elements_ Since this our _rst example_ we_ll simplify

matters by making the mesh uniform with hj _ h _ __N_ j _ __ _ _ _ _ _N_ and summing

ASj

_ AMj

_ and _V_ f_j separately_ Thus_ summing _______

N X

j__

ASj

_

N X

j__

_dj___ dj        

p

h _ _ __

__ _ __ cj__

cj __

The _rst and last contributions have to be modi_ed because of the boundary conditions

which_ as noted_ prescribe c_ _ cN _ d_ _ dN _ __ Thus_

N X

j__

ASj

_ _d_  

p

h

__        _c_      _ _d__ d_      

p

h _ _ __

__ _ __ c_

c_ __ _ _ _

__dN___ dN__         

p

h _ _ __

__ _ __ cN__

cN__ __ _dN__         

p

h

__        _cN__ _

Although this form of the summation can be readily evaluated_ it obscures the need for the

matrices and complicates implementation issues_ Thus_ at the risk of further complexity_

we_ll expand each matrix and vector to dimension N _ _ and write the summation as

N X

k__

ASj

_ _d__ d__ _ _ _ _ dN__      

p

h

__

_ _

_____

__

c_

c_

___

cN__

_____

__ Introduction

__d__ d__ _ _ _ _ dN__       

p

h

__

_ __

__ _ _

_____

__

c_

c_

___

cN__

_____

__ _ _ _ _d__ d__ _ _ _ _ dN__       

p

h

__

_ __

__ _

______

__

c_

c_

___

cN__

_____

__d__ d__ _ _ _ dN__          

p

h

__

_

______

__

c_

c_

___

cN__

_____

Zero elements of the matrices have not been shown for clarity_ With all matrices and

vectors having the same dimension_ the summation is

N X

j__

ASj

_ dTKc_ _______a_

where

K _

p

h

__

 __

__  __

__  __

_ _ _

_ _ _

_ _ _

__  __

__

_________

_ _______b_

c _ _c__ c__ _ _ _ _ cN__     T _ _______c_

d _ _d__ d__ _ _ _ _ dN__    T _ _______d_

The matrix K is called the global sti_ness matrix_ It is symmetric_ positive de_nite_ and

tridiagonal_ In the form that we have developed the results_ the summation over elements

is regarded as an assembly process where the element sti_ness matrices are added into

their proper places in the global sti_ness matrix_ It is not necessary to actually extend the

dimensions of the element matrices to those of the global sti_ness matrix_ As indicated

in Figure ______ the elemental indices determine the proper location to add a local matrix

into the global matrix_ Thus_ the    element sti_ness matrix Kj is added to rows

____ A Simple Finite Element Problem __

AS

_ _ d_

p

h

__       

__z_

c_ AS

_ _ _d__ d_    

p

h _ _ __

__ _ _

_ _z _

_ c_

c_ _

AS

_ _ _d__ d_    

p

h _ _ __

__ _ _

_ _z _

_ c_

c_ _

K _

p

h

__

 __

__  __

__ _ _

___________ Figure ______ Assembly of the _rst three element sti_ness matrices into the global sti_ness

matrix_

j _ _ and j and columns j _ _ and j_ Some modi_cations are needed for the _rst and

last elements to account for the boundary conditions_

The summations of AMj

and _V_ f_j proceed in the same manner and_ using ________

and _________ we obtain

N X

j__

AMj

_ dTMc_ _______a_

N X

j__

_V_ f_j _ dT l _______b_

where

M _

qh

_

__

_ _

_ _ _

_ _ _

_ _ _

_ _ _

_ _ _

_ _

_______

_ _______c_

l _

h

_

__

f_ _ _f_ _ f_

f_ _ _f_ _ f_

___

fN__ _ _fN__ _ fN

_____

_ _______d_

__ Introduction

The matrix M and the vector l are called the global mass matrix and global load vector_

respectively_

Substituting _______a_ and _______a_b_ into ______a_b_ gives

dT __K _M_c _ l        _ __ _______

As noted in Section ___ the requirement that ______a_ hold for all V _ SN

_ is equivalent

to satisfying _______ for all choices of d_ This is only possible when

_K _M_c _ l_ ________

Thus_ the nodal values ck_ k _ __ _ _ _ _ _N _ __ of the _nite element solution are deter_

mined by solving a linear algebraic system_ With c known_ the piecewise linear _nite

element U can be evaluated for any x using _____a__ The matrix K _M is symmetric_

positive de_nite_ and tridiagonal_ Such systems may be solved by the tridiagonal algo_

rithm _cf_ Problem  at the end of this section_ in O_N_ operations_ where an operation

is a scalar multiply followed by an addition_

The discrete system ________ is similar to the one that would be obtained from a

centered _nite di_erence approximation of ________ which is __      

_K _ D__c __

l_ ______a_

where

D _ qh_

_

_

_

_ _ _

_

_____

_ _l

_ h_

_

f_

f_

___

fN__

_____

_ _c __

_

_c_

_c_

___

_cN__

_____

_ ______b_

Thus_ the qu and f terms in _______ are approximated by diagonal matrices with the

_nite di_erence method_ In the _nite element method_ they are _smoothed_ by coupling

diagonal terms with their nearest neighbors using Simpson_s rule weights_ The diagonal

matrix D is sometimes called a _lumped_ approximation of the consistent mass matrix

M_ Both _nite di_erence and _nite element solutions behave similarly for the present

problem and have the same order of accuracy at the nodes of a uniform mesh_

Example __      ___ Consider the _nite element solution of

_u__ _ u _ x_ _ _ x _ __ u___ _ u___ _ __

which has the exact solution

u_x_ _ x _

sinh x

sinh _

_

____ A Simple Finite Element Problem __

Relative to the more general problem ________ this example has p _ q _ _ and f_x_ _ x_

We solve it using the piecewise_linear _nite element method developed in this section on

uniform meshes with spacing h _ __N for N _ __ _ _ _ _ _ __ Before presenting results_

it is worthwhile mentioning that the load vector ________ is exact for this example_ Even

though we replaced f_x_ by its piecewise linear interpolant according to _________ this

introduced no error since f_x_ is a linear function of x_

Letting

e_x_ _ u_x_ _ U_x_ _______

denote the discretization error in Table _____ we display the maximum error of the _nite

element solution and of its _rst derivative at the nodes of a mesh_ i_e__

jej_ __ max

__j_N

je_xj_j_ je_j_ __ max

__j_N

je__x_

j _j_ ______

We have seen that U__x_ is a piecewise constant function with jumps at nodes_ Data in

Table _____ were obtained by using derivatives from the left_ i_e__ x_

j _ lim_

__ xj__ With

this interpretation_ the results of second and fourth columns of Table _____ indicate that

jej__h_ and je_j__h are _essentially_ constants_ hence_ we may conclude that jej_ _ O_h__

and je_j_ _ O_h__

N jej_ jej__h_ je_j_ je_j__h

_ ________ ________ ______ __ _____

 _______ ________ ________ _____

__ ________ ________ _________ ____

_ ________ _______ _________ ____

__ ________ _______ ________ _____

_ ________ _______ _______ ____

Table ______ Maximum nodal errors of the piecewise_linear _nite element solution and its

derivative for Example ______ _Numbers in parenthesis indicate a power of ____

The _nite element and exact solutions of this problem are displayed in Figure _____ for

a uniform mesh with eight elements_ It appears that the pointwise discretization errors

are much smaller at nodes than they are globally_ We_ll see that this phenomena_ called

superconvergence_ applies more generally than this single example would imply_

Since _nite element solutions are de_ned as continuous functions _of x__ we can also

appraise their behavior in some global norms in addition to the discrete error norms used

in Table ______ Many norms could provide useful information_ One that we will use quite

_ Introduction

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0

0.01

0.02

0.03

0.04

0.05

0.06

Figure ______ Exact and piecewise_linear _nite element solutions of Example _____ on an

_element mesh_

often is the square root of the strain energy of the error_ thus_ using _____c_

kekA __ pA_e_ e_ _ _Z _

_

_p_e___ _ qe_            dx____

_ ______a_

This expression may easily be evaluated as a summation over the elements in the spirit

of ______a__ With p _ q _ _ for this example_

kek_

A _ Z _

_

__e___ _ e_    dx_

The integral is the square of the norm used on the Sobolev space H__ thus_

kek_ __ _Z _

_

__e___ _ e_    dx____

_ ______b_

Other global error measures will be important to our analyses_ however_ the only one

____ A Simple Finite Element Problem __

that we will introduce at the moment is the L_ norm

kek_ __ _Z _

_

e__x_dx____

_ ______c_

Results for the L_ and strain energy errors_ presented in Table ____ for this example_

indicate that kek_ _ O_h__ and kekA _ O_h__ The error in the H_ norm would be

identical to that in strain energy_ Later_ we will prove that these a priori error estimates

are correct for this and similar problems_ Errors in strain energy converge slower than

those in L_ because solution derivatives are involved and their nodal convergence is O_h_

_Table _______

N kek_ kek__h_ kekA kekA_h

_ _______ ________ _________ _____

 _________ ________ _________ _____

__ _________ ________ ________ _____

_ _________ ________ ________ _____

__ _________ ________ _______ _____

_ ________ ________ ______ _____

Table _____ Errors in L_ and strain energy for the piecewise_linear _nite element solution

of Example ______ _Numbers in parenthesis indicate a power of ____

Problems

__ The integral involved in obtaining the mass matrix according to ________ may_ of

course_ be done symbolically_ It may also be evaluated numerically by Simpson_s

rule which is exact in this case since the integrand is a quadratic polynomial_ Recall_

that Simpson_s rule is

Z h

_

F_x_dx _

h

_

_F___ _ _F_h__ _ F_h_         _

The mass matrix is

Mj _ Z xj

xj__ _ _j__

_j ___j___ _j   dx_

Using ________ determine Mj by Simpson_s rule to verify the result _________ The

use of Simpson_s rule may be simpler than symbolic integration for this example

since the trial functions are zero or unity at the ends of an element and one half at

its center_

_ Consider the solution of the linear system

AX _ F_ ______a_

_ Introduction

where F and X are N_dimensional vectors and A is an N  N tridiagonal matrix

having the form

A _

__

a_ c_

b_ a_ c_

_ _ _

_ _ _

_ _ _

bN__ aN__ cN__

bN aN

_______

_ ______b_

Assume that pivoting is not necessary and factor A as

A _ LU_ ______a_

where L and U are lower and upper bidiagonal matrices having the form

L _

__

_

l_ _

l_ _

_ _ _

_ _ _

lN _

_______

_ ______b_

U _

__

u_ v_

u_ v_

_ _ _

_ _ _

uN__ vN__

uN

_______

_ ______c_

Once the coe_cients lj _ j _ _ __ _ _ _ _N_ uj_ j _ __ _ _ _ _ _N_ and vj_ j _ __ _ _ _ _ _N_

__ have been determined_ the system ______a_ may easily be solved by forward and

backward substitution_ Thus_ using ______a_ in ______a_ gives

LUX _ F_ ______a_

Let

UX _ Y_ ______b_

then_

LY _ F_ ______c_

___ Using _______ and ________ show

u_ _ a__

lj _ bj_uj___ uj _ aj _ ljcj___ j _ _ __ _ _ _ _N_

vj _ cj _ j _ _ __ _ _ _ _N_

____ A Simple Finite Element Problem _

__ Show that Y and X are computed as

Y_ _ F__

Yj _ Fj _ ljYj___ j _ _ __ _ _ _ _N_

XN _ yN_uN_

Xj _ _Yj _ vjXj____uj_ j _ N _ __N _ _ _ _ _ _ __

___ Develop a procedure to implement this scheme for solving tridiagonal systems_

The input to the procedure should be N and vectors containing the coe_cients

aj _ bj _ cj _ fj _ j _ __ _ _ _ _ _N_ The procedure should output the solution X_

The coe_cients aj _ bj _ etc__ j _ __ _ _ _ _ _N_ should be replaced by uj_ vj _ etc__

j _ __ _ _ _ _ _N_ in order to save storage_ If you want_ the solution X can be

returned in F_

___ Estimate the number of arithmetic operations necessary to factor A and for

the forward and backward substitution process_

__ Consider the linear boundary value problem

_pu__ _ qu _ f_x__ _ _ x _ __ u___ _ u____ _ __

where p and q are positive constants and f_x_ is a smooth function_

____ Show that the Galerkin form of this boundary_value problem consists of _nding

u _ H_

_ satisfying

A_v_ u_ _ _v_ f_ _ Z _

_

_v_pu_ _ vqu_dx _ Z _

_

vfdx _ __ _v _ H_

_ _

For this problem_ functions u_x_ _ H_

_ are required to be elements of H_ and

satisfy the Dirichlet boundary condition u___ _ __ The Neumann boundary

condition at x _ _ need not be satis_ed by either u or v_

___ Introduce N equally spaced elements on _             x          _ with nodes xj _ jh_

j _ __ __ _ _ _ _N _h _ __N__ Approximate u by U having the form

U_x_ _

N X

j__

ck_k_x__

where _j_x__ j _ __ _ _ _ _ _N_ is the piecewise linear basis ________ and use

Galerkin_s method to obtain the global sti_ness and mass matrices and the

load vector for this problem_ _Again_ the approximation U_x_ does not satisfy

the natural boundary condition u____ _ _ nor does it have to_ We will discuss

this issue in Chapter __

 Introduction

____ Write a program to solve this problem using the _nite element method devel_

oped in Part __b and the tridiagonal algorithm of Problem _ Execute your

program with p _ __ q _ __ and f_x_ _ x and f_x_ _ x__ In each case_ use

N _ __ _ ___ and __ Let e_x_ _ u_x_ _ U_x_ and_ for each value of N_ com_

pute jej__ je__xN_j_ and kekA according to ______ and ______a__ You may

_optionally_ also compute kek_ as de_ned by ______c__ In each case_ estimate

the rate of convergence of the _nite element solution to the exact solution_

__ The Galerkin form of _______ consists of determining u _ H_

_ such that ______ is

satis_ed_ Similarly_ the _nite element solution U _ SN

_ _ H_

_ satis_es _______

Letting e_x_ _ u_x_ _ U_x__ show

A_e_ e_ _ A_u_ u_ _ A_U_ U_

where the strain energy A_v_ u_ is given by _____c__ We have_ thus_ shown that the

strain energy of the error is the error of the strain energy_

Bibliography

__        I_ Babu_ska_ J_ Chandra_ and J_E_ Flaherty_ editors_ Adaptive Computational Methods

for Partial Di_erential Equations_ Philadelphia_ ____ SIAM_

_          I_ Babu_ska_ O_C_ Zienkiewicz_ J_ Gago_ and E_R_ de A_ Oliveira_ editors_ Accuracy

Estimates and Adaptive Re_nements in Finite Element Computations_ John Wiley

and Sons_ Chichester_ ____

__        M_W_ Bern_ J_E_ Flaherty_ and M_ Luskin_ editors_ Grid Generation and Adaptive

Algorithms_ volume ___ of The IMA Volumes in Mathematics and its Applications_

New York_ _____ Springer_

__        G_F_ Carey_ Computational Grids_ Generation Adaptation and Solution Strategies_

Series in Computational and Physical Processes in Mechanics and Thermal science_

Taylor and Francis_ New York_ _____

__        K_ Clark_ J_E_ Flaherty_ and M_S_ Shephard_ editors_ Applied Numerical Mathemat_

ics_ volume ___ _____ Special Issue on Adaptive Methods for Partial Di_erential

Equations_

__        R_ Courant_ Variational methods for the solution of problems of equilibrium and

vibrations_ Bulletin of the American Mathematics Society_ ____ __ _____

__        J_E_ Flaherty_ P_J_ Paslow_ M_S_ Shephard_ and J_D_ Vasilakis_ editors_ Adaptive

methods for Partial Di_erential Equations_ Philadelphia_ ____ SIAM_

_          A_ Hrenniko__ Solutions of problems in elasticity by the frame work method_ Journal

of Applied Mechanics_ ____ ____ _____

__        C_ Johnson_ Numerical Solution of Partial Di_erential Equations by the Finite Ele_

ment method_ Cambridge_ Cambridge_ ____

___      D_L_ Logan_ A First Course in the Finite Element Method using ALGOR_ PWS_

Boston_ _____

_

_ Introduction

___      D_ McHenry_ A lattice analogy for the solution of plane stress problems_ Journal of

the Institute of Civil Engineers_ ____ _ _____

__        J_C_ Strikwerda_ Finite Di_erence Schemes and Partial Di_erential Equations_

Chapman and Hall_ Paci_c Grove_ ____

___      M_J_ Turner_ R_W_ Clough_ H_C_ Martin_ and L_J_ Topp_ Sti_ness and de_ection

analysis of complex structures_ Journal of the Aeronautical Sciences_ ____ __

_____

___      R_ Verf!urth_ A Review of Posteriori Error Estimation and Adaptive Mesh_

Re_nement Techniques_ Teubner_Wiley_ Stuttgart_ _____

Chapter _