# Numerical Methods for the Solution of Boundary Value Problems for Ordinary Differential Equations

by user

on
Category: Documents
3

views

Report

#### Transcript

Numerical Methods for the Solution of Boundary Value Problems for Ordinary Differential Equations
```Chapter 4
Numerical Methods for Solution of
Boundary Value Problems
for
Ordinary Dierential Equations
1 Introduction
Boundary Value Problems (BVP) for higher order ordinary dierential equations are frequently encountered in applications. These require the determination of a function of a
single independent variable satisfying a given dierential equation and subject to specied
values at the boundaries of the solution domain. Second and fourth order BVP are most
common in engineering applications.
Initial value problems require determination of the function subject to specied value(s)
at one end of the domain (typically t = 0). In contrast, in BVPs involving second order ODEs
and function values are specied at the two ends of the solution domain (typically x = 0 and
x = L). Many problems in engineering and science can be formulated as BVPs. Examples
include: steady state conduction heat transfer in a thin heated wire, electric potential inside
a thin conductor, deection of a thin elastic thread under load and many others.
In contrast with initial value problems, boundary value problems often have multiple
solutions or even fail to have a solution.
For instance, consider the problem of determining y(x) such that
d2 y + y = y00 + y = 0
dx2
By inspection (or by use of Frobenius's method), the general solution of this dierential
equation is of the form
y(x) = A sin(x) + B cos(x)
Now, if the boundary conditions are
y(0) = 0
y(=2) = 1
1
the unique solution to the BVP is y(x) = sin(x).
Alternatively, if the boundary conditions are
y(0) = 0
y() = 0
multiple solutions of the form y(x) = A sin(x) with A arbitrary exist.
Finally, if the boundary conditions are
y(0) = 0
y() = 1
no solution y(x) exists.
Three commonly used numerical methods for BVPs are: shooting, nite dierence and
nite element. All three can be used both for linear and nonlinear problems. The principles
of the three methods are described below.
2 Shooting Methods
Shooting methods are based on the transformation of the original boundary value problem
into an equivalent set of initial value problems. The solution of the original problem is
obtained directly from the solutions of the equivalent set of initial value problems.
2.1 Higher Order Equations and Systems of Equations
A higher order IVP can always be rewritten as a system of rst order IVPs.
dy = f (t; y)
dt
with
y(a) = (1 ; 2; :::; n)
where y(t) = (y1; y2; :::yn) and f (t; y) = (f1 (t; y); f2(t; y); ::fn(t; y))
For instance, the third order equation
d3y = f (x; y; dy ; d2y )
dx3
dx dx2
subject to the initial conditions
y(0) = 1
y0(0) = 2
y00(0) = 3
2
by introducing the new variables
y1 = y
dy
y2 = dx
d2y
y3 = dx
2
can be transformed into the equivalent system of three rst order equations
dy1 = y
dx 2
subject to
y1(0) = 1
subject to
dy2 = y
dx 3
y2(0) = 2
and
subject to
dy3 = f (x; y ; y ; y )
1 2 3
dx
y3(0) = 3
Uniqueness Theorem. If f is continuous and Lipschitz in all the variables y, the system
of IVP's has a unique solution y(t).
Therefore, any of the various Initial Value Problem algorithms used to solve single equations can be readily extended to deal with systems of equations. In Maple the solution of
systems is implemented with the dsolve command.
Recall that the fourth order Runge-Kutta method is a single step procedure which produces an approximation wi to the solution y(ti) of the initial value problem represented by
the equation
dy = f (t; y)
dt
on a x b subject to the condition
y(a) = 3
given by
wi+1 = wi + 61 (k1 + 2k2 + 2k3 + k4)
where
k1 = hf (ti; wi)
k2 = hf (ti + h2 ; wi + k21 )
k3 = hf (ti + h2 ; wi + k22 )
k4 = hf (ti+1 ; wi + k3 )
and h = t is the step size.
Runge-Kutta Method for Systems Algorithm.
Give number of equations m, the function f (t; y), initial condition i , endpoints a and
b, and the number of subintervals N .
For each variable yj , j = 1; 2; :::; m compute K1;j ; K2;j ; K3;j ; K4;j and the approximations wj
2.2 Linear Shooting
A boundary value problem (BVP) involving an ordinary dierential equation of order two
consists of determining the function y(x) for x 2 [a; b] such that
y00 = f (x; y; y0)
subject to the boundary conditions
y(a) = and
y(b) = Theorem of Uniqueness. If the function f (x; y; y0) as well as @[email protected] and @[email protected] are
continuous in the set D = f(x; y; y0) : a x b; ,1 y +1; ,1 y0 +1g and
if @[email protected] > 0 for all (x; y; y0) 2 D and [email protected][email protected] M for all (x; y; y0) 2 D where M is a
constant, then the BVP has a unique solution.
The shooting method is based on the idea of replacing the original BVP by an equivalent
set of initial value problems (IVPs) which are solved using standard IVP methods.
4
Corollary on Uniqueness of the Linear Problem. If the BVP is linear , i.e.
y00 = p(x)y0 + q(x)y + r(x)
in
axb
subject to
y(a) = and
y(b) = and the coecients p(x); q(x) and r(x) satisfy
1. p(x); q(x), and r(x) are continuous on [a; b],
2. q(x) > 0 on [a; b],
then the problem has a unique solution given by
y(x) = y1(x) + ,y (yb1)(b) y2(x)
2
where y1(x) and y2(x) are, respectively, the solutions to the following two IVPs
y100 = p(x)y10 + q(x)y1 + r(x)
in
axb
subject to
y1(a) = and
y1(a)0 = 0
and
y200 = p(x)y20 + q(x)
5
in
axb
subject to
y2(a) = 0
and
y2(a)0 = 1
Linear Shooting Algorithm
Give functions p(x); q(x); r(x), endpoints a; b, boundary conditions ; and the number
of subintervals N .
Set h = (b , a)=N ; u1;0 = , u2;0 = 0, v1;0 = 0, v2;0 = 1.
For i = 0; 2; :::; N , 1 set x = a + ih and use RK method for systems to determine
u1;i+1; u2;i+1; v1;i+1; v2;i+1.
Set w1;0 = and w2;0 = ( , u1;N )=v1;N .
For i = 1; 2; :::; N produce w1 = u1;i + w2;0v1;i y(x) and w2 = u2;i + w2;0v2;i y0(x)
2.3 Nonlinear Shooting
If the BVP is nonlinear, its solution cannot be obtained by linear combination of the solutions of the associated IVPs. Instead one must produce a sequence of solutions y(x; t)
corresponding to an IVP of the form
y00(x; t) = f (x; y(x; t); y0(x; t))
in
axb
subject to
y(a; t) = and
y0(a; t) = t
6
where t is a parameter of the method whose values must be selected so that
y(b; t) , = 0
If the preceeding dierential equation is dierentiated with respect to t and z(x; t) =
(@[email protected])(x; t) the second associated IVP is obtained
@f
z00 (x; t) = @f
z
(x; t) + 0 z0 (x; t)
@y
@y
in
axb
with
z(a; t) = 0
and
z0 (a; t) = 1
Therefore, nonlinear shooting proceeds by solving the above two IVP's and Newton's
method to determine subsequent iterate values for tk as follows
tk = tk,1 , y(b;z(tb;k,t 1) ,) k,1
Nonlinear Shooting Algorithm.
Give function f (x; y; y0), endpoints a; b, boundary conditions ; , the number of subintervals N , the tolerance TOL and the maximum number of iterations M .
Set h = (b , a)=N ; k = 1, TK = ( , )=(b , a),
For k M set w1;0 = , w2;0 = TK , u1 = 0, u2 = 1.
For i = 1; 2; :::; N set x = a + (i , 1)h and use RK method for systems to determine
w1;i; w2;i, u1; u2.
If jw1;N , j TOL stop else set TK = TK , (w1;N , )=u1 (Newton's method) and
continue until k = M .
7
3 Finite Dierence Methods
In the nite dierence method one starts with the dierential formulation and by a process
of discretization transforms the problem into a system of interlinked simultaneous algebraic
equations that then must be solved in order to determine an approximation to the desired
solution. Discretization consists of rst introducing a mesh of nodes by subdividing the
solution domain into a nite number of sub domains and then approximating the derivatives
in the boundary value problem by means of appropriate nite dierence ratios which can
be obtained from a truncated Taylor series expansion. As a result, system of interlinked
simultaneous algebraic equations is obtained that then must be solved in order to determine
an approximation to the desired solution.
3.1 Linear Problems
Consider the linear BVP
y00 = p(x)y0 + q(x)y + r(x)
in
axb
subject to
y(a) = and
y(b) = Introduce a mesh in [a; b] by dividing the interval into N + 1 equal subintervals of size
h. This produces two boundary mesh points x0 and xN +1 , and N interior mesh points xi ,
i = 1; 2; ::; N . The values y(x0) = y(a) and y(xN +1) = y(b) are known but the values y(xi),
i = 1; 2; ::; N must be determined.
From the Taylor expansions for y(xi+1) and y(xi,1) the following centered dierence
formula is obtained
h2 y(4)( )
y00(xi ) = h12 [y(xi+1) , 2y(xi) + y(xi,1)] + 12
i
for some i 2 (xi,1 ; xi+1).
Similarly, an approximation for y0(xi) is
2
y0(xi ) = 21h [y(xi+1) , y(xi,1)] , h6 y(3) (i)
8
for some i 2 (xi,1 ; xi+1 ).
If the higher order terms are discarded from the above formulae and the approximations
are used in the original dierential equation, the second order accurate nite dierence
approximation to the BVP becomes a system of simultaneous linear algebraic equations
,( wi+1 , 2hw2i + wi,1 ) + p(xi )( wi+1 2,h wi,1 ) + q(xi )wi = ,r(xi )
or, alternatively
,(1 + h2 p(xi))wi,1 + (2 + h2 q(xi))wi , (1 , h2 p(xi ))wi+1 = ,h2 r(xi)
for i = 1; 2; ::; N .
In matrix notation the system becomes
Aw = b
where A is a tridiagonal matrix.
Theorem of Uniqueness. If p(x); q(x) and r(x) are continuous and q(x) > 0 on [a; b]
the problem Aw = b has a unique solution provided h < 2=L where L = maxaxb jp(x)j.
Further, if y(4) is continuous on [a; b] the truncation error is O(h2).
Linear Finite Dierence Method Algorithm.
Give functions p(x); q(x); r(x), endpoints a; b, boundary conditions ; , and the number of subintervals N + 1.
Set h = (b , a)=(N + 1); x = a + h a1 = 2 + h2 q, b1 = ,1 + (h=2)p, d1 = ,h2 r + [1 +
(h=2)p].
For i = 2; ::; N , 1 set x = a + ih ai = 2 + h2 q, bi = ,1 + (h=2)p, ci = ,1 , (h=2)p,
di = ,h2 r.
Set x = b , h, aN = 2 + h2 q, cN = ,1 , (h=2)p, dN = ,h2 r + [1 , (h=2)p]
Solve tridiagonal system for zi .
w0 = , wN +1 = , wN = zN .
For i = N , 1; :::; 1; set wi = zi , uiwi+1
Stop.
9
3.2 Nonlinear Problems
Let the BVP
y00 = f (x; y; y0)
subject to
y(a) = and
y(b) = be such that
1) f , @[email protected] = fy and @[email protected] = fy0 are continuous on D;
2) @[email protected] > 0;
3) Constants k and L exist such that k = [email protected][email protected] and L = [email protected][email protected] j
Then, a unique solution exists. Introducing nite dierences and since f is nonlinear,
the following system of nonlinear algebraic equations is obtained
,( wi+1 , 2hw2i + wi,1 ) + f (xi; wi; wi+1 2,h wi,1 ) = 0
for each interior mesh point i = 1; 2; ::; N . The nonlinear system has a unique solution as
long as h < 2=L.
The Jacobian matrix associated with this system is tridiagonal and it is given by
8
>
< ,1 + h2 fy0 i = j , 1;
J (w1; :::; wN )ij = > 2 + h2fy i = j;
: ,1 + h fy0 i = j + 1;
2
j = 2; ::; N
j = 1; ::; N
j = 1; ::; N , 1
Newton's method applied to this problem requires rst the solution of the linear system
J v(k) = b
to produce v(k) which is then added to w(k,1) to obtain w(k), i.e.
w(k) = w(k,1) + v(k)
for each i = 1; 2; ::; N .
Nonlinear Finite Dierence Method Algorithm.
Give functions p(x); q(x); r(x), endpoints a; b, boundary conditions ; , the number
of subintervals N + 1, the tolerance TOL and the maximum number of iterations M .
10
Set h = (b , a)=(N + 1), w0 = , wN +1 = .
For i = 1; ::; N set wi = + i[( , )=(b , a)]h.
Set k = 1, and while k M compute ti; ai; bi; ci; di for i = 1; 2; :::; N .
Solve tridiagonal system for zi .
Set vN = zN , wn = wN + vN .
Compute vi = zi , uivi+1 and wi.
If kvk TOL stop, else continue until k = M .
3.3 Exercises
3.4
Solve numerically the following boundary value problem:
in x 2 [1; 2] and subject to
d2y + 2 dy , 2 y , sin(ln x) = 0
dx2 x dx x2
x2
y(1) = 1
and
y(2) = 2
3.5
Solve numerically the following boundary value problem:
in x 2 [1; 3] and subject to
d2y + 1 y dy , 1 x3 , 4 = 0
dx2 8 dx 4
y(1) = 17
and
y(3) = 43
3
11
4 Numerical Approximation Methods based on Variational Formulations
4.1 Introduction
Mathematical representation of physical laws is possible in two alternative but equivalent
formulations
Dierential formulation
Variational formulation
The dierential form is usually obtained by rst applying the conservation principle on a
small box of material followed by a limiting process in which the volume of the box is made
innitesimally small. The variational formulation is the integrated form of the conservation
principle over the entire volume of the material under consideration. Each formulation can
be obtained from the other one by a straightforward process.
For instance, mechanical equilibrium in an elastic body of volume , loaded by conservative force at its boundary is described by the equilibrium equations
@ij = 0
@xj
Where ij (x; y; z) is the stress tensor eld which is to be determined by solving the above
equations subject to the imposed force vector (per unit area) Fi at its boundary ,. In
practice, one always determines rst the displacement eld ui, then the strain tensor eld
eij (x; y; z) and last the associated stress tensor eld.
Alternatively, the variational representation of the same problem involves a functional given by
=
Z
Z
Wd
, Fi uid,
,
Where W = @ij [email protected] is the strain energy function.
The desired solution ui(x; y; z) is selected from amongst an appropriately chosen set of
functions and is such that the value of is minimal. This is the principle of minimum
potential energy.
As a second example, for a solid material of volume , inside which thermal energy is
generated at a rate g per unit volume and which is maintained at xed temperature B = 0
at its boundary ,, under conditions of thermal equilibrium, the dierential representation is
the steady state heat conduction equation
,r kr + g = 0
12
inside , where k and are the thermal conductivity and the temperature inside the solid.
To determine the function (x; y; z) one must solve the above dierential equation subject
to the condition
= B = 0
at the boundary ,.
The corresponding variational representation of the same system involves a functional given by
Z
Z
= (r) krd
+ gd,
,
The desired solution (x; y; z) is selected from amongst an appropriately chosen set of functions and is such that the value of is minimal. This is the principle of conservation of
thermal energy.
Except for the most elementary situations, the solution of most engineering problems
cannot be obtained using classical analytical methods of applied mathematics and numerical
solution procedures must be employed. Two commonly used methods of numerical solution
of engineering problems are the nite dierence method and the nite element method.
In the nite element method, one starts with the variational formulation and then by
a process of nite element function representation transforms the problem into a system of
interlinked simultaneous algebraic equations that then must be solved in order to determine
an approximation to the desired solution. The process of nite element function representation is based on an idea rst introduced by Rayleigh and further developed by Ritz and
Galerkin. It is based in rst representing the desired solution as a linear combination of
suitable basis functions and then proceeding to the minimization of the functional resulting
in a set of algebraic equations for the unknown coecients.
Although both methods transform the original innite dimensional problem into a nite
dimensional systems of algebraic equations, the specic form of the resulting systems may
or may not be the same when the two methods are applied to the same problem.
Therefore, the implementation of the nite element method for the solution of engineering
problems involves two key steps:
The variational statement of the problem.
The approximate solution of the variational problem by means of nite element functions.
T
4.2 Basic Theory of the Finite Element Method
Consider the following dierential formulation of a boundary value problem
2
, ddxu2 = ,u00 = f
13
for x 2 [0; 1] and subjected to the boundary conditions
u(0) = 0
and
u0(1) = 0
The task is to determine the function u(x).
The equivalent variational
formulation is obtained by introducing a function v(x) satisR1 2
fying v(0) = 0 and 0 v dx < 1, but otherwise rather arbitrary and then integrating the
dierential equation, i.e.
Z
1
0
,u00 vdx =
Z
1
0
fvdx
Integrating the left hand side by parts readily yields the variational or weak form
a(u; v) = (f; v)
where
a(u; v) =
and
(f; v) =
Z
1
u0v0dx
1
fvdx
0
Z
0
We say that the function v belongs to a space of functions V (v 2 V ) dened as follows
V = fv 2 L2 (0; 1) : a(v; v) < 1; v(0) = 0g
Then, it can be shown that, as long as f is a continuous function of x (i.e. f 2 C 0 ([0; 1]))
and u possesses continuous derivatives up to the second order (i.e. u 2 C 2 ([0; 1]), the desired
solution u(x) 2 V is the function satisfying
a(u; v) = (f; v)
for all v 2 V .
The Ritz-Galerkin approximation of the above problem is obtained by introducing a
subspace of functions S V and taking the function uS from there rather than from V , such
that
a(uS ; v) = (f; v)
14
for all v 2 S . It can be shown that this problem has a unique solution which is obtained by
solving the nite dimensional system of algebraic equations
KU = F
Where K is the stiness matrix which can be regarded as a nite dierence operator.
The estimated error involved in the Ritz-Galerkin approximation process is given by
k u , uS kE = mink u , v kE : v 2 S
where the energy norm k : kE is dened as
q
k v kE = a(v; v)
Moreover, one can also show that
k u , uS k 2 k f k
q
where k v k= (v; v) is the L2 norm, and is a small number.
In the nite element method one uses piecewise polynomial spaces of functions. For
instance, partition the interval [0; 1] into elements joined at the nodes xi such that 0 = x0 <
x1 < x2 ::: < xi::: < xn = 1. Now introduce the basis functions i with value of 1 at node i
zero at all other nodes and linear between nodes i , 1 and i and i and i + 1. The interpolant
vI is dened as
vi =
n
X
i=1
v(xi )i
Then, it can be shown that
k u , uI kE Ch k u00 k
where h is the size of the largest element and the constant C is independent of h and u.
The equivalence between the solution of a boundary value problem and that of the associated variational problem follows from the uniqueness theorem (Rayleigh-Ritz):
Let p 2 C 1[0; 1], q; f 2 C [0; 1] with p(x) > 0 and q(x) 0 for 0 x 1. Then, the
function y(x) 2 C 2[0; 1] with y(0) = y(1) = 0 is the solution to the dierential equation
d (p(x) dy ) + q(x)y = f (x)
, dx
dx
if and only if y(x) is the unique function which minimizes the integral
I [y] =
Z
0
1
fp(x)[y0(x)]2 + q(x)[y(x)]2 , 2f (x)y(x)gdx
15
In nite element practice, one searches for the solution among a set of functions produced
by linear combination of a set of linearly independent andPwell dened basis functions i(x),
i = 1; 2; :::; n. The solution y(xP) is rst approximated by ni cii(x) and then the coecients
ci are chosen so that I [y] = I [ ni cii(x)] is minimized, i.e.
@I = 0
@cj
for each j = 1; 2; ::; n. P
Substituting y(x) ni ci i(x) into I [y], dierentiating and introducing the minimization
Ac = b
where
aij =
Z
1
0
[p(x)0i (x)0j (x) + q(x)i(x)j (x)]dx
are the entries of the n n stiness matrix A,
bi =
Z
1
0
f (x)i(x)dx
is the load vector and c = (c1; c2 ; ::; cn)t is the vector of unknown coecients.
The simplest kind of basis functions is the set of piecewise linear polynomials dened by
rst introducing in [0; 1] a collection of nodes x0; x1 ; :::; xn+1 with hi = xi+1 , xi and then
making
8
>
0;
0 x xi,1
>
>
< x,h x ,1 ; xi,1 < x xi
i(x) = > x +1,,1 x ; x < x x
> h
i
i+1
>
: 0;
xi+1 < x 1
i
i
i
i
for each i = 1; 2; ::; n. This approximation is continuous but it is not dierentiable. Note
that i(x)j (x) = 0 and 0i(x)0j (x) = 0 except when j is equal to i , 1, i, or i + 1 and this
makes A tridiagonal. Evaluation of the matrix entries requires the calculation of 6 integrals
for each node. One can approximate the functions p, q, and f by rst order interpolating
polynomials and then integrate numerically.
Piecewise Linear Rayleigh-Ritz Algorithm.
Give functions p(x); q(x); f (x), endpoints 0; 1, boundary conditions = 0; = 0, and
the number of nodes n + 2.
For each node i = 1; :::; n dene the basis function i(x).
16
For i = 1; :::; n , 1 compute the integrals Qj;i for j = 1; :::; 6.
Compute entries in the tridiagonal matrix.
Solve tridiagonal system for the unknown coecients ci.
Compute (x) = Pni=1 cii(x).
Piecewise linear basis functions produce a continuous solution which is however nondierentiable. Dierentiability of the approximation can be obtained by using spline functions and a Cubic Spline Rayleigh-Ritz algorithm can be obtained.
4.3 Rayleigh-Ritz-Galerkin Finite Element Method
Recall that the boundary value problem consisting of nding the function u(x) such that
2
, ddxu2 = ,u00 = f
subject to the boundary conditions u(0) = u(1) = 0 is equivalent to the variational problem
consisting of determining the function u(x) such that
a(u; v) = (f; v)
where
a(u; v) =
and
Z
1
0
,u00dx =
(f; v) =
Z
1
0
Z
1
0
u0v0dx
fvdx
and also to the problem of nding the function u(x) that minimizes the functional
I (u) = 12 a(u; v) , (f; v)
where v is a reasonably well behaved but rather arbitrary function.
When one looks for approximate solutions to the integrated forms of the problem by
the Galerkin method or by the Ritz method, one seeks for solutions in a nite dimensional
subspace of the actual solution space.
In the Ritz method, one searches for an approximation to the function that minimizes
the functional. The method starts by introducing a set of basis functions fi(x)gNi=1 and
expresses the approximate solution as the trial function uN u as
uN =
N
X
i=1
ai i
17
The functional obtained using this approximation then becomes
I (uN ) = 12 a(uN ; v) , (f; v)
Now, the function v being arbitrary, is selected as the test function
N
X
v = vN = uN =
i=1
ai i
and is substituted into the functional expression to give
N
N
N
X
X
X
I (uN ) = 12 a( aii; ai i) , (f; aii)
i=1
i=1
where
N
X
a( ai i;
i=1
N
X
i=1
aii) =
and
(f;
N
X
i=1
ai i) =
Z
i=1
Z
1
0
1
0
(ai
N d
X
i
i=1
dx ) dx
2
N
X
f ( ai i)dx
i=1
Finally, the unknown coecients ai ; i = 1; 2; :::; N are determined by solving the system
of algebraic equations obtained from the extremum conditions. i.e.
@I = @I = :::: = @I = 0
@a1 @a2
@aN
which in matrix notation are simply written as
Ka = F
In the Galerkin method, one seeks directly for an approximate solution of a(u; v) =
(f; v). Like in the Ritz method, in the Galerkin approach one starts by introducing a set
of basis functions fj (x)gNj=1 satisfying the stated boundary conditions and expresses the
approximate solution as the trial function uN u as
uN =
N
X
j =1
aj j
However, the test functions are selected so as to be identical to the basis functions, i.e.
fvigNi=1 = fi(x)gNi=1 .
18
The resulting discrete form of the problem is then
N
X
i=1
aj a(j ; i) = (f; i)
for all i = 1; 2; :::; N . Here
a(j ; i) =
Z 1 d d
j
i
and
(f; i) =
dx dx dx
0
Z
1
0
fidx
for i = 1; 2; :::; N
Using matrix notation, the above is simply written as
Ka = F
Solving the algebraic problem yields the values of the coecients ai; i = 1; 2; :::; N .
The system of equations obtained with the Galerkin method is identical to the one obtained using the Ritz method. However, the Galerkin method is of more general applicability
since it works even when the problem of minimizing the functional has no solution. The nite
element method is obtained by implementing the Galerkin (or Ritz) procedure with a very
particular choice of basis functions.
4.4 Numerical Implementation of the Finite Element Method
The nite element method is a direct implementation of the Galerkin (or Rayleigh-Ritz)
procedure in which the chosen basis functions are nite element basis functions (also called
sometimes nite element shape functions). These functions are of very simple form (piecewise
polynomials of low order are most common) and their most distinctive feature is that they
are nonzero only in a small subregion of the computational domain, i.e. they have a local
character. Another important feature of the nite element basis functions is that the basis
function associated with a particular element partially overlap with those of immediately
element method by simply increasing the number of sub domains (i.e. by decreasing their
size) without increasing the complexity of the polynomials used to represent the solution
inside each individual subregion. Moreover, the matrix associated with the resulting system
has a sparse structure and ecient solution methods are applicable. This is the well known hformulation of the nite element method. Last, but not least, since the coordinate functions
of adjacent subregions overlap partially with each other, the resulting system of algebraic
equations can be constructed and assembled together by the computer just as it is created.
19
In implementing the nite element method one starts by using simple interpolation functions as basis functions at the element level and then proceeds to represent the solution in the
entire domain by collecting the contributions associated with each element. This is possible
thanks to the mathematical property that allows computation of integrals over domains as
sums of integrals over sub domains spanning the original domain.
The basis functions involved at the element level are called local basis functions while
the ones that produce the solution in the entire domain are called global basis functions.
Clearly, the two sets of basis functions are closely related.
The role of the local basis functions is to generate the value of the approximated quantity
inside the element from the values at the nodes by interpolation. Assuming again a one dimensional system and the simple subdivision introduced above, consider an arbitrary interior
element ei where the index i denotes the global node number. If the values of the required
solution at xi and xi+1 are u(xi) = ui and u(xi+1) = ui+1, respectively, the approximate
values of u(x), u(x)e inside the element are calculated by simple interpolation as follows
ue (x) = uii1 + ui+1i2
Here, the local nite element basis functions for element e , i1 and i2 are dened as
i1 = xi+1h , x
i
i
i
i
i2 = x ,h xi
i
where hi = xi+1 , xi is the element size and all the positions are measured in the global
coordinate system. Note that the indices 1 and 2 on the local basis functions refer to the
local node numbers for the element.
For simplicity, assume all nite elements are of equal size (uniformly spaced mesh) so
that h1 = h2 = ::: = hi = h. Now, the approximate solution on the entire domain uh(x) is
represented as a linear combination of the nite element basis functions, i.e.
uh(x) = u11 + u22 + ::: + uii + uN N =
N
X
i=1
uii
However, note that here, i; i = 1; 2; :::; N are global nite element basis functions.
There is a simple relationship between local and global nite element basis functions,
specically, in element 1,
1 = 11
in element 2
( 1
2 = 22 ifif xx 22 ee1
2
1
20
in element ei,
(
i,1
i = 2i ifif xx 22 eei,1
i
1
and in element eN
N = N2
Using the introduced notation the Galerkin nite element method equations can be expressed as
N
X
i=1
uj a(j ; i) = (f; i)
for all i = 1; 2; :::; N . Here
a(j ; i) =
Z 1 d d
j
i
dx dx dx
0
and
(f; i) =
Z
1
0
fidx
for i = 1; 2; :::; N .
Using matrix notation, the above is simply written as
Ku = F
Solving the algebraic problem yields the values of the nodal values ui; i = 1; 2; :::; N . The
matrix K is called the nite element stiness matrix, the colum vector F is the force vector
and the column vector u is the vector of unknown nodal values of u.
The same thing can be done to express the Ritz nite element method equations using
the new notation. Specically, the functional I (uh) is given by
N
N
N
X
X
1 X
I (uh) = 2 a( uii; uii) , (f; uii)
i=1
i=1
i=1
where
N
X
a( uii;
i=1
N
X
i=1
uii) =
and
(f;
N
X
i=1
ai i) =
Z
1
0
Z1X
N
di
(
0 i=1
N
X
ui dx )2dx
f ( uii)dx
21
i=1
Finally, introduction of the extremum conditions
@I = @I = ::: = @I
@u1 @u2
@uN
yields the system
Ku = F
which is identical to the one obtained using Galerkin's method.
4.5 Examples
4.5.1 One Dimensional Problem in Cartesian Coordinates
Consider the boundary value problem
, ddxu2 = 20x3
2
subject to the boundary conditions u(0) = u(1) = 0. Find an approximate solution using
the Galerkin nite element method with three equal size elements and the following global
nite element basis functions:
(
x 2 [0; 1=3]
1 (x) = 10 , 3x for
elsewhere
8
>
for x 2 [0; 1=3]
< 3x
2(x) = > 2 , 3x for x 2 [1=3; 2=3]
:0
elsewhere
8
>
< 3x , 1 for x 2 [1=3; 2=3]
3(x) = > 3 , 3x for x 2 [2=3; 1]
:0
elsewhere
(
x 2 [2=3; 1]
4 (x) = 30x , 2 for
elsewhere
The variational representation of the problem is readily obtained by rst multiplying the
given dierential equation by a function v and subsequently integrating from x = 0 to x = 1.
After integration by parts, the result is
a(u; v) = (f; v)
22
where
a(u; v) =
and
(f; v) =
Z 1 du dv
0
Z
1
0
dx dx dx
20x3vdx
Now, let us proceed to the implementation of the nite element method. Here there are
three elements e1 ; e2; e3 and four nodes i = 1; 2; 3; 4. The nodal values of u are ui; i = 1; 2; 3; 4
and since u is specied at x = 0 and x = 1, then u1 = u4 = 0 and the only unknowns are u2
and u3. The representation of the trial function u in terms of the nodal values then becomes
"
4
X
uh = uii = u22 + u33 = [2; 3] uu2
3
i=1
#
In the Galerkin method the test function v is given directly in terms of the basis functions.
In the present case, since u1 = u4 = 0, v is simply given by
"
vh = 2 + 3 = [1; 1] 2
3
#
Substituting the above into the variational formulation yields
"
a(uh; vh) = [1; 1] kk22 kk23
32
33
and
"
(f; v) = [1; 1] FF2
3
where
Z
#"
u2
u3
#
#
Z 2=3 d
d
2 2
k22 = ( dx ) dx =
( 2 )2 dx = 6
dx
0
0
k23 = k32 =
k33 =
1
Z 1 d d
2
3
0
Z
1
0
Z
dx dx dx =
( d3 )2 dx =
dx
1=3
Z
23
2=3
d2 d3 dx = ,3
dx dx
( d3 )2dx = 6
1=3 dx
1
F2 =
and
F3 =
Z
0
1
Z
1
0
20x32dx =
20x3 3dx =
Z
Z
1=3
0
2=3
1=3
20x3 3xdx +
Z
2=3
1=3
20x3(3x , 1)dx +
Z
20x3(2 , 3x)dx = 4 + 20 = 10
81 81 27
20x3(3 , 3x)dx = 49 + 131 = 60
81 81 27
2=3
1
Therefore, the nite element equation is
"
#" # "
6
,
3
u2 =
Ku = F = ,3 6
u3
Finally, solving for u2 and u3 yields
80
u2 = 243
u3 = 130
243
10
27
60
27
#
4.5.2 One-Dimensional Problem in Cylindrical Coordinates
Consider the boundary value problem consisting of nding the function u(r) satisfying
1 @ @u
r [ @r (r @r )] = 0
subject to
,( @u
@r )r1 = q1
at r = r1, and
u(r2) = uS
The problem is a mathematical model of steady state conduction through a hollow cylindrical
shell of a material with unit thermal conductivity, exposed to heat ux input at its inner
radius r1 and maintained at xed temperature u2 = u(r2) at its outer radius r2 . Here u(r)
is the temperature, r is the radial coordinate and k is the thermal conductivity.
The analytical solution to this problem is
u(r) = u2 , q1 r1 ln( rr )
2
Using the calculus of variations one can show that the solution to the above problem is
the function that minimizes the following functional
Z 2 Z r2 du
Z 2
1
2
I (u) = 2
k( dr ) rdrd , q1 u1r1d
0
r1
0
24
where u1 = u(r1) is the yet to be determined temperature on the inner surface of the
cylindrical shell and is the angular coordinate.
To implement the nite element Ritz-Galerkin method, nite element basis functions are
required. The key is to assumes a simple, easily integrable form for the solution in terms of
the spatial coordinate r and still to be determined parameters. This is then substituted into
the expression for the functional and the minimum of the functional is sought by equating
the derivatives of the functional with respect to the unknown parameters to zero. The nal
result is a system of algebraic equations that must then be solved for the values of the
parameters. The procedure is exactly as the classical Ritz method but the basis functions
used are nite element basis functions.
The simplest thing to do at rst is to consider a single element representation. In this
case, the nite element is a radial segment spanning the entire shell thickness of the hollow
cylinder. The ends of the nite element are the locations of the inner and outer surfaces.
These are also the two nodes associated with the nite element. Since there is only one
element, global and local quantities are the same and no distinguishing superscripts are
needed.
Assume now that the solution can be approximately expressed as linear combination of
the nodal values of u
uh(r) = 1 (r)u1 + 2(r)u2 = u
where 1; 2 are the local nite element basis functions,
= [1; 2]
is the global basis functions vector (a row matrix) and
"
u = uu12
#
is the (column) vector of nodal values of u.
The basis functions are selected to be simple linear functions of the radial distance r and
are given by,
1 = rr2 ,,rr
2
1
2 = rr ,,rr1
2
1
As mentioned before, since the element spans the entire domain, 1 and 2 are also global
basis functions associated with the (nodal) positions r1 and r2 .
25
du = Ba
dr
where the row matrix B is
1 d2
B = [ d
dr ; dr ] = [B1 ; B2]
Using the given vector notation and substituting into the expression for the functional
yields
Z 2 Z r2
Z 2
T BT Burdrd ,
u
u1q1 r1d
I (uh) = 12
0
r1
0
One can now proceed to carry out the integrals. The integral over is direct since none
of the quantities in the integrands depends on the angular coordinate. Further, the vector
operations can be performed to split the integrand in the remaining integral into four terms,
i.e.
Z r2
1 2 2
) u1 + ( d1 )( d2 )u1u2+
(1)
I = [( d
dr
dr dr
r1
d d
d
( 2 )( 1 )u2u1 + ( 2 )2 u22]rdr , 2u1q1 r1
(2)
dr dr
dr
The individual integrals are easily obtained and are
Z r2 d
r +r
( 1 )2u21 = 2 1 u21
2(r2 , r1)
r1 dr
Z
r2
r1
Z
r2
r1
( d1 )( d2 )u1u2 = , r2 + r1 u1u2
dr dr
2(r2 , r1)
( d2 )( d1 )u2u1 = , r2 + r1 u2u1
dr dr
2(r2 , r1)
Z
( d2 )2u22 = r2 + r1 u22
2(r2 , r1)
r1 dr
Therefore, the functional nally becomes
I = 24((rr2,+rr1)) (u21 , 2u1u2 + u22) , 2u1q1 r1
2
1
r2
26
or, in matrix notation
where K, the stiness matrix is
I = 12 uT Ku , uT F
(r + r ) 1 ,1 2
1
11 K12
K= K
K21 K22 = (r , r ) ,1 1
2
and the distributed force vector is
F=
2q1r1
0
1
!
Since the u2 = uS is given, the only independent parameter aecting the value of the
functional is u1. Therefore the condition for the functional to attain its minimum value is
simply
@I = 0
@u1
Carrying out the dierentiation and solving for u1 yields
r2 , r1 )
u1 = uS + 2q1(rr1(+
2 r1 )
In this case, the nite element equation
Ku = F
reduces to a single equation in a single unknown.
27
```
Fly UP