w, the integral expression yields

Numerical solution of one-dimensional diffusion equation

236

b

g2 ( x) dx = 0 ’ g( x) = 0 on a ¤ x ¤ b. (14.7)

a

This follows immediately from the observation that the square of a function g( x)

is always greater than or equal to zero for any value of x, i.e. g2 ( x) ≥ 0, such that

the integral of g2 ( x) over the domain a ¤ x ¤ b must be greater than or equal to

zero, i.e.

b

g2 ( x) dx ≥ 0, (14.8)

a

and can only be equal to zero if g( x) = 0 for all a ¤ x ¤ b.

Effectively, the method of weighted residuals transforms the requirement that a

function, say g( x), must be equal to zero on a given domain at an in¬nite number

of points into a single evaluation of the integral, that must be equal to zero.

Using the method of weighted residuals, the differential equation, Eq. (14.1), is

transformed into an integral equation:

b d du

+f dx = 0,

w c (14.9)

dx dx

a

which should hold for all weighting functions w( x). The term between the square

brackets contains second order derivatives d2 u/dx2 of the function u. As has been

outlined in the introduction, approximate solutions of u are sought by de¬ning an

interpolation of u on the domain of interest and transforming the integral equation

into a discrete set of linear equations. De¬ning interpolation functions that are

both second-order differentiable and still integrable is far from straightforward,

in particular in the multi-dimensional case on arbitrarily shaped domains. Fortu-

nately, the second-order derivatives can be removed by means of an integration by

parts:

b b b

dw du

du

’ dx + wf dx = 0.

c (14.10)

wc

dx dx dx

a a

a

This introduces the boundary terms:

b

du du du

= ’w( a) c + w( b) c

wc . (14.11)

dx dx dx

a x=a x=b

At the boundary either u is prescribed (i.e. the essential boundary condition at

x = a) or the derivative c du/dx is prescribed (i.e. the natural boundary condition

at x = b). Along the boundary where u is prescribed the corresponding ¬‚ux, say

pu , with

du

pu = c , (14.12)

dx x=a

237 14.4 Polynomial interpolation

dx x=b = P is known. For

is unknown, while along the other boundary the ¬‚ux c du

the time being the combination of boundary terms in Eq. (14.11) is written as

b

du

B= w c = ’w( a) pu + w( b) P, (14.13)

dx a

realizing that pu is considered as yet unknown. The term B is introduced as an

abbreviation for the boundary contribution.

Consequently, the following integral equation, known as the weak form,

results:

b b

dw du

dx = wf dx + B.

c (14.14)

dx dx

a a

Eq. (14.14) is the point of departure for the Finite Element Method, and is called

the weak form, because the differentiability requirements imposed on u have been

reduced: in the original differential equation, Eq. (14.1), the second-order deriva-

tive d2 u/dx2 appears, while in the weak form, Eq. (14.14), only the ¬rst-order

derivative du/dx has to be dealt with.

14.4 Polynomial interpolation

Suppose that at a ¬nite number of points xi in the domain , the function values

ui = u( xi ) are known, then a polynomial approximation, denoted by uh , of u( x)

can be constructed. The polynomial approximation uh of degree n ’ 1 can be

constructed by

uh ( x) = a0 + a1 x + a2 x2 + · · · + an’1 xn’1 , (14.15)

if u is known at n points. The coef¬cients ai can be identi¬ed uniquely and

expressed in terms of ui , by solving the set of equations:

⎡ ¤⎡ ¤⎡ ¤

2 · · · xn’1

1 x1 x1 u1

a0

1

⎢ 1 x x2 · · · xn’1 ⎥ ⎢ a ⎥ ⎢ u ⎥

⎢ ⎥⎢ 1⎥ ⎢ 2⎥

2

⎢ 2 · · · xn’1 ⎥ ⎢ a ⎥⎢ ⎥

2 2

⎢ 1 x 3 x3 ⎥⎢ 2 ⎥ = ⎢ u3 ⎥ . (14.16)

⎢ ⎥⎢ ⎥⎢ ⎥

3

⎢. . . ⎥⎢ . ⎥ ⎢ . ⎥

.. . . . ¦⎣ . ¦ ⎣ . ¦

. .

⎣. . . .

. . .

···

2 n’1 an’1 un

1 xn xn xn

An example is given in Fig. 14.2, where, in the domain xi’1 ¤ x ¤ xi+2 , a third-

order polynomial (dashed curve) is used to approximate a given function (solid

curve) based on the function values ui’1 , ui , ui+1 and ui+2 .

Clearly, the coef¬cients ai are linearly dependent on the values ui , therefore the

polynomial may be rewritten in terms of ui by

Numerical solution of one-dimensional diffusion equation

238

u

x

xi “1 xi xi +1 xi +2

Figure 14.2

Solid line: u( x), dashed line: polynomial approximation of u( x).

n

uh ( x) = Ni ( x) ui , (14.17)

i=1

where the functions Ni ( x) are polynomial expressions of order n ’ 1 in terms of

the coordinate x. These functions Ni ( x) are called shape functions because they

de¬ne the shape of the interpolation of uh , for instance linear, quadratic etc. To

illustrate this, consider a ¬rst-order polynomial on the domain [ x1 , x2 ], with u( x)

known at x1 and x2 . In that case

x ’ x1 x ’ x1

uh ( x) = 1 ’ u1 + u2 , (14.18)

x2 ’ x1 x2 ’ x1

implying that

x ’ x1 x ’ x1

N1 = 1 ’ N2 =

, . (14.19)

x2 ’ x1 x2 ’ x1

Rather than approximating u( x) with a single polynomial of a certain degree over

the entire domain , the domain may also be divided into a number of non-

overlapping subdomains, say e . Within each subdomain e a local polynomial

approximation of u may be constructed. A typical example is a piecewise lin-

ear approximation within each subdomain e , see Fig. 14.3. Consider one of the

subdomains e = [ xi , xi+1 ]. Then within e the function u( x) is approximated by

uh ( x) | = N1 ( x) ui + N2 ( x) ui+1 , (14.20)

e

with, in conformity with Eq. (14.19):

x ’ xi x ’ xi

N1 = 1 ’ N2 = . (14.21)

xi+1 ’ xi xi+1 ’ xi

More generally, if within a subdomain e an n-th order polynomial approximation

of u is applied, the subdomain should have n + 1 points at which the function u is

239 14.5 Galerkin approximation

u ui

ui +1

x

xi “1 xi xi +1 xi +2

©e

Figure 14.3

Solid line: u( x), dashed line: piecewise linear approximation of u( x).

known. For instance, to use a second-order (quadratic) polynomial, the subdomain

should cover at least three consecutive points, for example e =[ xi , xi+2 ], such

that

uh ( x) | = N1 ( x) ui + N2 ( x) ui+1 + N3 ( x) ui+2 . (14.22)

e

The subdomain e within which a certain polynomial approximation is used is

referred to as an element. The points at which the values of u are de¬ned are

called nodes.

The shape functions Ni may not be chosen arbitrarily. The most stringent

requirement is that uh must be interpolated continuously over the total domain

(the ¬rst-order derivative of uh ( x) should exist). Suppose that the nodes xj within

an element are numbered j = 1, . . . , n and that the associated shape functions Ni

are numbered i = 1, . . . , n. In that case, for consistency, the shape functions must

be chosen such that

Ni ( xj ) = δij , (14.23)

with

δij = 0 if i = j

= 1 if i = j. (14.24)

14.5 Galerkin approximation

To transform the weak form into a linear set of equations in order to derive an

approximate solution the following steps are taken.

Step 1. Element division As shown in the previous section the domain may

be split into a number of subdomains e , elements, and within each element a

Numerical solution of one-dimensional diffusion equation

240

polynomial interpolation can be made of the function u. An example of such an

element division is given in Fig. 14.4. This distribution of elements is called a

mesh. Then, the integration over the domain can be performed by summing up

the integrals over each element. Consequently, Eq. (14.14) yields:

Nel Nel

dw du

dx = wf dx + B,

c (14.25)

dx dx

e e

e=1 e=1

where Nel denotes the number of elements.

Step 2. Interpolation Suppose that the domain has been divided into three

linear elements, as depicted in Fig. 14.4. Then the nodal values ui may be collected

in an array u:

∼

⎡ ¤

u1

⎢ ⎥

⎢ u2 ⎥

u=⎢ ⎥. (14.26)

⎣ u3 ¦

∼

u4

The unknowns associated with each of the elements are collected in the arrays

e

ue , such that

∼

u1 u2 u3

u1 = u2 = u3 =

, , . (14.27)

∼ ∼ ∼

u2 u3 u4

So, it is important to realize that each particular element array ue contains a subset

∼

of the total, or global, array u. Within each element array ue a local numbering

∼ ∼

may be used, such that for this particular example with linear elements: