Lab.
v. Scheepsbauwkuntic
Technische
Hogeschool
Delft
NUMERICAL SOLUTION
OF THE
SHALLOW WATER EQUATIONS
BY A
FINITE ELEMENT METHOD
NIEK PRAAGMAN
2 .1081NUMERICAL SOLUTION
OF THE
SHALLOW WATER EQUATIONS
BY A
FINITE ELEMENT METHOD
Numerical solution of the shallow water equations
by a finite element method.
Proefschrift ter verkrijging van
de graad van doctor in de
technische wetenschappen
aan de Technische Hogeschool Delft!,
op gezag van de rector magnificus
prof.dr.ir. F.J. Kievits,
voor een commissie aangewezen
door het college van dekanen
te verdedigen op
woensdag 19 september
1979'
te 14.00 uur door
Nikolaas Praagman
wiskundig doctorandus
geboren to Drachten
Dit proefschrift is goedgekeurd
door de promotor
prof.dr. EL van Spiegel.
met dank aan:
Phil en Sue Anderson,
Bert Broere,
Gonny Scheurkogel  v. Werkhoven
en voor al
Guus Segal.
Contents
0 Introduction
1 The shallow water equations 3
1.0 introduction 3
1,1 the system of partial differential equations 3
2 The finite element method to solve partial differential equations 14
2.0 introduction 14
2.1
mathematical tools 142.2
the linear heat equation 192.3
a linear transport equation 282.4
the linear convection diffusion equation 332.5
observations 363 Methods for solving ordinary differential equations 38
3.0
introduction 383.1 truncation error and stability 38
3.2
three classes of methods for ODEs 403.3
some examples 413.4
systems of ordinary differential equations 483.5
extended stability regions 513.6 observation 59
4 Stability and numerical results 60
4.0 introduction 6o
4.1
methods to investigate stability 604.2
lumping 684.3
results for system2.(26)
694.4
results for system2.(53)
724.5
results for system2.(68)
822
Onedimensional incompressible flow 845.0
introduction 845.1 boundary conditions 84
1,2
the Galerkin formulation 87VIII
the system of ODEs 88
5.4 numerical stability 89
LI
computations and results 90la
remarks (other methods) 9611
conclusions 966 A numerical treatment of two dimensional shallow water flow 97
6.0 introduction 97
6.1
boundary conditions 976.2 the Galerkin formulation 101
6.3
the system of ODEs 1026.4 stability 103
computations and results 105
7' 'Conclusions and remarks' 131
Appendix 1 133
Appendix 2 138
References 143'
dit prodschrift werd gedrukt bij. drukkerij & uitgeverij
EGNER by Den Helder
Q
Introduction.Mathematical models for two dimensional flow have been developed at various locations.
For instance, in each country that borders the North Sea, an operational
numerical model exists. (see [1], [14], [16], [17], [27], [37], [47]).
One of the reasons for the current interest in these models is the discovery of oil and gas in the North Sea.
For the exploitation of these minerals local currents and water levels are important. A mathematical model can be used to compute these quantities. Another reason for the interest is caused by the pollution of the North Sea.
The pollution of coastal seas is being studied extensively at the moment.
Numerical models have been developed for pollution problems too. However models that treat dispersion and advection of pollutants need data that
are generated by numerical models for two dimensional shallow water flow.
In this thesis we investigate a numerical method for two dimensional shallow water flow that is based on a combination of the finite element method
(f.e.m.) and the Runge Kutta method (RK4). The basis of this solution method is the combination of the f.e.m. with an 'ordinary differential equations solver" which can be utilized for all evolution problems in partial
differential equations.
Using a number of model problems special features of the method are studied and in the final chapters results are given for practical problems. In addition the influence of the various terms of the numerical model are
studied in chapter 6.
Below is a summary of the items treated in each chapter:
Chapter 1: Short derivation of the partial differential equations for the nonsteady motion of longwaves.
Chapter 2: Introduction of the Galerkin and finite element method on the basis of three model problems. Application of these methods results in Systems of Ordinary Differential Equations. (BODE.). Chapter 3; Survey and discussion of numerical methods for solving SODEs.
1
2Chapter 4: Solution of the SODEs of chapter 2 by use of the methods discussed in chapter 3.
Chapter 5: Computation of water levels and velocities in the case of
1dimensional flow.
Chapter
6:
Computation of water levels and velocities for the Southern Part of the North Sea.Investigation of the influence of the various terms.
Chapter 7: Remarks and conclusions.
A.F.E.P.
In chapters 4, 5, 6 we discuss the results of computations. All these
computations have been carried out by programs that are based on the finite element package A.F.E.P. ([42], [43]).
The shallow water equations.
1.0 Introduction
In this chapter we give a short derivation of the shallow water equations. These partial differential equations are derived by a vertical integration
of the equations of motion and the continuity equation
in
an Euleriancoordinate system.
In the final system of equations coefficients occur that are not completely
determined by theoretical considerations. Hence in the numerical ',computations
empirical values have to be used.
I1.t
The system of partial differential equations.Consider a coastal sea and define cartesian coordinates x and y in the
zaxis
"horizontal"plane of thea
equation for a fluid reads:
xaxis
fig. 1.1
+ 1 + 2 (pv), +
.(pw) =
0at ax ay
In
(11 11,17
and w are velocity components: in the x, y and z direction), and p is the density of the, fluid.The general form of the equatIons of motion is (in a fixed coordinate
system)
a2 a.
a(POI + 1(q1.1
) (puv) + r(puw) = pFat
ax,
ay az a , , tpv).+' kpuv), +
_{ax}
L (
ay Pv + az PVW )1 = pFLt
a, _{a.} a.ow/ + (puw) + (avw) + ,(pw2) pF .
atax
ay
asundisturbed water surface. The zaxis is chosen perpendicular to
yaxis
and yaxis. (See fig.1.1)The general form of the continuity
a,
F F and Fz represent the, components of the external forces per unit of mask
(1)
3(0u)
+ 2 ) x,
4In this thesis especially the southern part of the North Sea is considered.
In that case the following forces turn out to be important: (DO])
The Coriolisforce Gravity
The pressure force Friction forces Wind forces
The contribution of these forces to F , F and F is:
x y z Coriolis: Fx : = 2Q3v  2n2w F : = 2Q w  2Q3u 1 Fz : = 2Q 2u  2Q1 v with
= (0 ,2 n
).
the earth rotation vector..Gravity: Fx : = 0 F : = 0 Fz : = g with
g : the acceleration of gravity.
Pressure F 22x ;:c F : =
IR
Y P aY Fz :p az
with p : pressure(6)
zS97
PO,
412),
Friction (molecular viscosity)
Fx
vu
(8)
F vAv= yaw
with ,2 ,2 A =2 + a
+ a2ax
ay
ay2 322v d, = kinematic viscosity coefficient
Remark: Windforces and bottom friction will be treated after the
vertical integration of the horizontal components of, the equation of motion.
Insertion of the various terms in equations (2)1, (3) and (4) leads to
a,
a , 2, 3 , , a , ,43.7 IbuA + i; kpu )
+ 7 OUVJ 4 7
kpuvi = p(20131720w)  + pvnuay
az
32 anda,
, a I , a , 2, a .ov) +
tpuvm + ) + ,pvw) = pc2/1w2111 3t ax ay az t+T1 1 r, 71,I
uo dt1=101Here we assume that It is possible to make the time averaging interval I;
so large that the turbulent fluctuations vanish but also to make T so small that the important tidal variations are not significantly influenced by this procedure.

12
_{+ OvAli}ay
a
, 1a I
a a , 2, , + kpuw, + urn+, + Ow 4 = bi2.02u201  pg  +
at ax
ay
az
aZSince we are only interested in the tidal motion, we apply, to equations, (1), (9), (10) and (11), a moving timeaverage to separate the short term
turbulent fluctuations from the longterm tidal oscillations. (DO:, [49)IN.
So let the velocity field (u,v,w) be separated into a mean part (1i,7,,;) and a fluctuating part (u' ,v' ) in such a way that
t+T 1 r
T
u (It = U. Fz : [18], 5=(17)
6Before we apply this timeaveraging technique to equations (1), (9), (10) and (11) we make the assumption that the density o is a constant. Although
the density is a function of, for instance, temperature, salinity and microscopic movements and so consequently a function of time and space,
the assumption is reasonable if the sea under consideration is wellnixedw
Results of the time averaging procedure are:
Equation '(1; aa aT;
ax
ay Equation (9) _3uI
_{4.} 2 a,  3
1 _aP.Cu( ) + _{iT. (uv)} +
i; (a) .
2a3.7,  2a2w  +VU+
at ax
t+T a t+T t+T
I

( 
f(le2 di + 
f
(u'v' )
di + 2f
(u,w, )1 dt/T ax 3y 3z t t t
Equation 06)
_
1 li 1(16')
av a017) ,
(,7)2 +h
1(17;)= 20137'  2.03u 
7, ay + vA,7 +, at ax ay t+T t+T t+T I a' r ,,2
1 3f
,ltr,
_  {
tU V h uT +
7/. j Orkll dT + 2anf
(V'Wt ) diE 1 T ax ay _{t} _{t} 1 t Equation (11)1Ai
a()
a(,
a
2
at
ax
''w, ji orw, (w ) = 2a2 _ 20 1 16  P az Vt14t+T
at+T
t+T
3 t
=11(ew')
dr + (viw) at+f
dx,
(vT)2 dt),
ay t az t
The last term in each of the three averaged equations results from the
nonlinearity of equations (2),
(3)
and (4).In view of their similar effect onthe model, it is generally accepted to
express these terms as molecular diffusion.
Equations (15),
(16)
and (17) then become( respectively.:1(11()r,
A(15

7+
a()
+ a (W1.) + 2 (W2) =_{2Q2a _}
2Q1  g  + (v+vT) AW
at ax uw ay az
p az
The problems of the turbulence phenomenon have not been eliminated since the magnitude of VT is not known.
In some models VT is used as a control parameter. In other models VT is given a value in accordance with more or less theoretical considerations.
In section
6.6.3
we return to this problem.Equation (20) can be simplified. If we estimate and compare the orders the different terms it turns out that as a first approximation the equation
1
2i
_ _ 0
p az
holds.
The estimates are based on the fact that we study the tidal motion in
a shallow sea.
In such a case the characteristic horizontal length scales Lx and L are
much larger than the vertical length scale L.
If the horizontal velocities are characterized by U then it follows from
equation (14) that the characteristic vertical velocity W will be of the
order (if L = L = L ). x y
UL
Now Lz « L and hence W « U.The characteristic time scale T is not independent of U and L; the relation
between T, U and L
is
L r U T.
The terms of equation (20) can be estimated if we use characteristic values
for the various variables and appropriate magnitudes for the coefficients:
_{(CM}
(18) ail 3 2
T. 4' Tx (u )
a (717') a(u)
= 223v 22QW  2p 33. ) (v+vT, Au ( 19 ) 77 + kuv) + 2()
.57(v.)
 = 2Q1W 22 u 3p ay
+ (v+vT) + +00)
(22) K2178
423)1 aw af (\i 4. 9()
A_ (72, LzU2 ? 2 ataz
ay
'1az
.
A in s 2Q2U  2Q1 MO4m
s2 = g  VO m 52 1Ai=
I LzU 1 2 vTA14 m 10ms
LzvA; m 106 LzU
1
Ins2
L 2
Lz
Theorder of the values that should be substituted for a1' Q2 and v is wellknown. So these coefficients cause no problems. However the value of VT
is unknown and even the order is not quite clear. In [303 it is proposed
to take
3 2 1
vT = 10 m s, .
This rough estimate suffices at this stage. For if we compare the orders
of the different terms in equation (20) for the case that
1
(130): U in
we see that the pressure term is the only one that can balance the gravity term g. The socalled hydrostatic distribution of the pressure, equation
(211)1, is the approximation that results from this analysis.
This hydrostatic approximation provides a possibility to eliminate the
pressure from equations (18) and (19). First we define the waterdepth h
and the waterlevel !&e:
The depth t is the distance between the x9r plane and the bottom.
= The waterlevel is the distance between the xy plane and the
mater surface at the time considered.
Let pa be. the atmospheric pressure then integration of ,(21) gives:
P = Pg"°
Pa ,(210 'h25N '("261 (31) + L2 p 1.37)
The variations of the atmospheric pressure are small compared to the variations of pgC, so we take:
and
Equations (14), (18) and (19) will also be integrated over the vertical. Before carrying this out we define:
u =
u dz
H 11 v = dz h with H = h + C.It is also assumed that the surface of the fluid is a streamline. This means
that any particle which is on the surface remains on it. ([46]).
The interface sea air is given by the equation
36) C(x,y,t) = z
A particle will remain on this surface if at z = &:
iL
+ :
+at ax ay
The equation of the bottom
h(x,y,t) = z
provides, in the same way, in the relation at z =
ah  ah  ah
at
u T)c. " wThe integration formula (n stands for x, y or t)
&
,
21ah
f
'
. dz = 1f
f dz  f(C) _{ f(h) } an an an an h h is also used. (314) (35) I 323' ay9A,

10Integration, of
ilk)
leads tot

E J au r aydz +dz +
;(x,,y,E
,t )  iTg(x,,y ,h ,t = Dx ay hi hCombination of this relation with (37),
(38)
and_{(39)}
gives:7,

ah ah+ + ykx,yA,t) + u(x,y,,t) + +I
at ax Hy !at ax ay
ah , ah , an
it
+ ukx,y,h,t) 5.27 + vkx,y,h,t) i,. +C  iC
f
11=az +
f
ay  ax ay Az x :41 h J am , = + a (Hu) ay ley)  0 at DxIf we neglect the slow variation in time of' the bottom, i.e. if we take
at
we find
ii
3 , aKliu, + Kay)
at ax
Integration of the nonlinear equations /415)[ and (19) poses .more problems.
The various terms of (18) give:
r 2C ah .02) Du dz = 2
f
dz)  a(x,Y,t,t)]  a(x,y,h,t) =at at at h h a ahKliu)
u(x,y,C,t) at at rThe term ( j u2 dz is evaluated
3x taking
a
a
then, since dz h K4o)On')
0 + ah = ay = h = 0 we find: Tr(x,y,h,t)(43)
 a 2 Tc. (u ) dz =2 (
f
712 dz) 7.12(x,y,E,t) ax ax h h 2 ah af
2 u (x,y,h,t) + S ) dz = ax ax h ar 2 2al
2 ah 3 2= (Hu )  u (x,y,E,t) ax  u (x,y,h,t) + ( f S dz.)...
ax ax
h
In the same way
a  a ,
(44)
f
(uv) dz = (Hi)  uv(x,y,E,t) +ay ay ay h E 3h 3 U(x,y,h,t) + (
f
ST, dz) ay ay_{h}
and 3 (45) 5T (uw) dz = 71;(x,y,E,t) ZW(x,y,h,t)
h
Evaluation of the right hand side terms requires a few more asumptions. The first result of the order analysis was
W « Ti , so
«
In the Coriolis force we neglect therefore the term 222W. Besides that in the following we write f instead of
2523
in accordance with common practice. We study the tidal motion in a relatively small sea and
as a consequence f and g are assumed to be constants.
Integration of the righthand side terms results in:
5
f.i.dz=fHv
h E (47) 5 g a dz = H gax
ax h and: 2 2 u a u,f
(v+vT) AZ
dz =f
(v+v 2) dz h h T)k +ax2
ay(46)
au(v+vT)()
z=E au,
(v+v ), ) T az z=h u +(50)
(51)
12Combination of (46), (47), (48),
(37)
and(38)
gives:(49)
(H71) + (HZ2) + a (H7;) = at ax ay n2n2=PH;
_{(v+vT)}f
(=1 +=a) dz
ax
h ax ay3y2 a ( f& a2 dz)  af
dz)(v+vT)11)
ax ay h 3u  (v+v T az z=E z=hThe u and v stem from the variation of the velocity components over the
depth. The effect of this variations is, as was the case for the turbulent
fluctuations, similar to the effect of molecular diffusion. Hence the terms are represented in the model as diffusion. However as a consequence the actual value of the diffusion coefficient is undetermined.
The last two terms of equation (4 9) are used for wind and bottom influences.
Suppose that
is a twodimensional wind velocity vector near to the airwater interface.
Following [30] we put aa(v+vT) az = =2 =2 = K u(u +v ) z=E
The quantity K, which will be one of the empirical values in the final
system of equations, is a function of various variables, for instance 71,
V and H.
In the same way we put, following [9], [27]
(v+v ) 2L1
T az
 2 2
g u(u +v )
z=h C2
In this formula the quantity C, which is called the De Chezy value, still has to be specified.
If we use
(50)
and (51) and put+
13C 2 ,2 C li a Au = (+vT)
f
+ °
u) Az  3 (f
a
Az)  3 (f
dz), ax ay 2 2 h ax a1 h h(so the coefficient a represents all the diffusionlike processes:
molecular diffusion, turbulence and vertical fluctuations,) then we obtain
I(1.17.0 +
1 (H2) + 3 (Huv) = fH;  gH il + HaATA + K:(7.1242) +at 3x 3y ax
gu(u +v ) 32 32
 ,with A = +
C2 ax2 3y2
Using equation (40) and dividing by H then:
, ) 2 2, au

au  au
2i
1101 +v ic:(:242)i (52) + u  + v fv + g + g  aAu o
at ax ay ax C2(h+C) (h+C)The same procedure applied to equation (19) gives:
3;  6(11242)i ( 5 3
)u v
f17/ g aAv _{=0}
at ay g C2(h+C) (h+C)Equations (41), (52) and (53), the socalled "shallowwater equations",
together with appropriate initial and boundary values constitute a mathematica: model. In the next chapters a numerical technique for the solution of such a model is analysed and in the last chapters results of the application
of this technique to the model are discussed. I
The finite element method to solve partial differential equations.
2,S1 Introduction
In this chapter three examples will be given of the discretization of a partial differential equation of evolution type, by a combination of the standard Galerkin method and the finite element method. ([45]).
The partial differential equations, with initial and boundary values,
that will be treated are respectively; (1) The linear heat equation
I'L
Ut = ux.
+g
u(x,0) = u0(x) u(0,t) = f0(t) u(1,t) = fl(t)
(2) A linear transport equation..
ut + ux = g u(x,0) = u0(x) u(0,t) u f0(t) U= u(x,t) g = g(x,t) a Notation: ut = u X E COO] t c [0,T]. x c 0 = [0,1] t c [0,T]
(3) A linear convection diffusion equation
1
a,b 0 ut + aux  buxx = g u(x,0) = u0(x) u(0,t) = r0(t) u(1,t) = f1(t) x c 0 = [0,1] t a [O,T]Existence and uniqueness of the solution are assumed for each of the
problems (1), (2) and (3). Proofs can be found in the textbooks ([6],
[73
[28]).
2.1 Mathematical tools
In sections 2.2  2.4 the problems (1), (2) and (3) will be discretized with respect to the (space)variable x. Before we start this discretization
some mathematical tools are introduced.
xx
15The space part of each of the problems (1), (2) and (3) has to be found in some special spaces that represent the smoothness properties of solution and boundary conditions. These spaces are called spaces of admissible
functions.
Without loss of generality the boundary conditions may be assumed to be homogeneous (See section 2.2.1)
The special spaces that will be used in this chapter are:
(C2 c MiN, with boundary an)
The space L201): L2(0) =
{0(x)I
(04)0 < 0,1 with (u,v),v= f
uv dx the space HI k2): 1 "H. (0)= f4)(x)I
7(0 ' L2(n)} with vg,3z1 ".''axN
the space H1(2) ={0(x)I
0 c H1(1),0(x)I
= 01 0 IxEart the space 1 1 HE =(0(x)1
c H (0),0(x)I
= 01 with E c an.It is wellknown that
L2(a)
is a Hilbert space when the inner product(5)
is used.
1
The spaces Ho, H1 and HE are called Sobolev spaces. They are also Hilbert spaces with respect to the inner product
(4J
f(l6q
_{Xu,v)1 =}
fi(uv + vu.Vv)
dx = (u,v)0, + (01.1,Vv)0Two norms, related with
(5)
and(6)
are respectivelyNul I g =
(u,110'0and
= (u,U)Iir
One can prove that boundary
sense for the spaces H1 ,
H1o
unstable ([24], (45]). Such
16conditions containing derivatives do not make,
1
and HE, since these boundary conditions are boundary conditions are called "natural".
The spaces of admissible functions will be discretized using the f.e.m.
1
The discretization of H1, HI and HE will be called respectively V1,
Vo
o 1
and V.
The space V is defined in the following way
FirsttheregionSisdividedintofiniteelementse.,which is called
a' triangulation 'h of S (C 4
]), .
V1 01) is the space of functions satisfying:.
1,0
The restriction of0(x)
to an element ej is a polynomial of degree k,kIN
;(ii), The function
0(x) is
continuous over S.,The spaces V1(Q) and
V1(a)
are subspaces of V1(2). The elements of thesespaces satisfy respectively
0(x) = 0, X E DO if 0(x) a V10(S)
and
0(x) = 0,,
x E
E if 0(x)1 c 4(12)Frequently vectors x in the complex space are considered. For these
vectors the ;12norm:
111.112 =
will be used.
Eigenvalues of various matrices play a role in the analysis of the problems that will be treated. We shall use theorems that estimate the largest and smallest eigenvalues of matrices. The matrices are assumed to contain real valued elements.
Theorem 2.1 (Cerschgcrin): Every eigenvalue of a matrix A lies in the complex plane in at least one of the circular discs with centre a.. and radii
11
/
ji
Theorem
2.2
A Hermitian matrix A has a complete set of orthogonaleigenvectors.
Proofs of these two theorems can be found in [52].
C C
Theorem 2.3 Let C be a real matrix with eigenvalues Ai, X2,..., AN.
Let
s = ;(c+cT)
and
T = _{(CCT)}
When the eigenvalues of Si and T are respectively
and
then
(i)1
min X! 5fie(JS)
max AicI 1
1E'
. 1(ii) 
max A. 5 Im(),) 5 max lAiliEI
iEI
= {1,2,...,N}
1.
V 5; c I
17(10)
T x x
/ lail2
i=1
A real symmetrix matrix has realvalued eigenvalues, hence:
T
x Sx
T
x Sx
(ii) min A.
1 T1
s
5 max A.i
x x
i
In the same way we find for the skew symmetric matrix T
is real
is imaginary
5 max 'ATI
icI 1
18Proof: The matrix S is symmetric and the matrix T is skewsymmetric. S and iT are Hermitian, hence both have a complete set of orthogonal
eigenvectors.
Ss
Let x1)_2t...,_vx x_ be the normalized eigenvectors of S. Then every vector x
can be written as
X= La. x.
_{1 1}
i=1
Using the orthogonal properties of the set
we find
1 1E.I.
N TSxi1
A7. lail2 x =Remark: When A is an eigenvalue of a skewsymmetric matrix, then X =  A
is an eigenvalue too. Hence (11) estimates all the eigenvalues of T.
Consider Cx = Ax
x x
Tx x x
T T x x x x and T T x Sx x Tx A = T T
xx
xx
The theorem now follows from (10) and (11)
Remark: For the estimation of the AS and AT' i c I, theorem 2.1 can
utilized.
2.2 The linear heat eauation.
2.2.1 A property of the p.d.e. (1)
In this section we derive a qualitative property of the solution of problem (1). The boundary conditions of (1) are supposed to be homogeneous.
(This is not a restriction: consider a function v satisfying the nonhomogeneous boundary conditions.
Define w= u  v then 'satisfies' { wt = w.. 4 g 
vt 4
v.. (12) w(x,0) = u0(x)  v(x0) w(0,t) = 0 w(1,t) = 0 x c = [0,1] tE [0,T]
then: 19T x Cx xT(i(C+CT) + (CCT))xSo nonhomogeneous boundary conditions result in a different forcing
function g*:
ODE g = g  vt +V
In order to derive the property of the solution u(x,t) of problem (1) we multiply both sides of equation (1) with the function u(x,t), and integrate the result over the domain
20a = [0,1]:
(ut,u)0 = (uxx,u)0 (g'u)0
Partial integration gives:
(u u) + (u ,u ) = fu ,u]1 + (g,u)
t, 0
Since the boundary conditions are homogeneous we find:
1 d
((u,u)o) + (ux,ux)0 = (g,u)0
Using the Poincareinequality (u,u)0 5 i(ux,ux)0
and the Schwarzinequality
(u,g)2 s (u,u) (g,g)
0 0 0
we find
1 d
(19)(u,u)0 + 2(u,u)0
5 (u,u)
(g,g)2075 dt
Hence
(2o)
mull20)
_{211u1120 s Hullo}ilglio
Divide (20) by Hullo and multiply with e2t then:
d 2t
(e Hullo) 5 e2t11g110
Integration of this inequality leads to
(21) 1/11(t)110 e2tliu(0)110 ( e2(st) s)1
lie,11
k 10 ds V tE [0,T].
0)Equation (21)
expresses
that, in some sense, problem (1) is stable.In section 2.2.4 we shall
see
that this property is also valid for thesolution of the discretized form of (1) that we derive in section 2.2.2.
2.2.2 The Galerkin formulation of (1)
In this section we apply the standard Galerkin method in order to discretize problem (1).
2
f2a)
In [7] it is shown that the solution u of problem (1) is an element of
the infinite dimensional space H:(0) for all t
E [0,T].
Since V1(0) is an internal approximation of H1(1), ([
4 ]),
we approximateu by the solution uh of the problem: (standard Galerkin approach,. see [L5],)n
. h
Find u such that
(uh vh) (uh vh)
in
vhkx' 10' 70,
cv:(0I
, V t c [O,T]and
uh,(0)
Define a basis fOilNi..) in VIO(Q)
Since uh c V1(i1), we have: N 1(23Y
21(x,t)
=uco .(x)
Take with uh(t) V vh V1(0), d h(24).1
ui(t))(0 1 4) 0 1=1N, do. do.
V un r 1 11 e, ) L i 'dx '0"ei'wj 0 1=1Initial conditions for system (24] are found from the equation
N,
= u (x.) = ui(0) i(xj)
0 j, i=1
Remark: System (24) is the basis for all numerical computations. Another way to derive this system is to use the classical Galerkinweighted residual approach:
Let
uh(x,t) =y u.(r) 0.c.
in
V j, e I.
1
v =
i
h, , in 1(22). This. leads to the SODE with the unknown functions uitt1, i e
0 + given. 1=1 21: ( dt dx uh(x.,0)
Define: N du
Nd20.
_{h} R(x,t,ui) = 4).(x) dt u. 2  g. i=1 i=1 1 dx From(RAI) =0
_{j 0}VjEI
we find, after partial integration and usage of the boundary conditions,
h \ N dui(t) N
0.
dc1). dx 0 + 1uh
(1'' 'L) )0 i dx = (64j )0 ' V j E. I. 1 dt ((Pi4j i=1 i=1In the technical literature one often uses formulation (22) with vh
replaced by 5u.
2.2.3 The system of ODEs.
Equations (25) can be written:
.hd
(26)
Mu
+ S uh = F ( = )dt
M and S are N x N matrices with elements:
in.. (4).,0 )
J j 0
and _{dOi dO.}
s.
_{ji}= (
_{dx 'dx}1)
_{0}while
u(t) = Eu.1.11(t),11121(t),...,4(t)7T
and
= (g4j)0
TheevaluationofF.can be difficult for arbitrary g(x,t), therefore we
use the approximation
N+1
(31) g = gk(t) 0k(x)
k=0
Example_
Consider, for simplicity, the space Vo(Q) of linear elements, i.e. the
space of functions
0(x)
H1(S1) that are linear over each element e.(see section 2.2.1)
=
722.
1 1
034)I
A basis KO IN of _{V1(t0,11) is found if we define the functions}
Oi(x) by
K101 0W
is 'linear over each element e(ii)
01.(xj./dis
(Kronecker delta)Remark: Since 01([0,1]) c
HI({0,100,
the0 (x) c
Therefore the functions
00(x): and Own(x) do not belong to the space
1 ,
V't
(C
For this special case the matricesTI and S become (see also E241, [32], _{[1451}
and appendix 2). and
(,35/
S =
Ax M = 6 while F is given by21
1 21 0\\\
1 2 1 01 2
=23For the triangulation of( = {0,1]we use an equidistant stepsijize Ax and define jar x. = 1
e. = Ex.x.].
1 1 H1o(0).Remark: The last two vectors of (36) result from the boundary conditions. Both vectors contain only boundary values.
2.2.4 A property of the BODE
(26)
The elements of the matrices M and S are respectively
i,j I and AU.
=ijE
ij dx 'dx 0 ' T 'From these expressions itis immediately clear that both M and S are symmetric. Besides that we have:
Theorem 2.4 The matrices M and S are positive definite.
Proof: Consider for arbitrary x:
X
Mx = lam x. ( m..x.) =x.(
(0i(x),0.(x)) x.) =
j . J 0 J = j1 i=1 j=1 2 = (x.(1).(x),
x.0.(x))
= II 1x.0.(x)110
J J 0 j=1 1 1 i=1 J1Since the functions
4i(x)
form a basis of V1(0) they are linearly0independent, hence
x.c1).(x) 0 0 when x
0
0and we conclude that
xTMx 0
V x
oThe matrix S yields a similar property
240
01
0 11N+1 11N+1(36)
F =
M gig2
Ax 6 go o o Ax  6 . gN 0 gN+1 . i'.1o 1 0 I 0 . 1 Ax o o i=1 0T, x ox
=I
x.( 2 s. .x.) = i=l _{J1} N dO.= Vx)
i=1 j=1 dx dx 0 j dØ. N de. (/
Xiat,
2
x. .;.1)0 i=, dep. II Jc.1112
_{dx} 0i=r
By Poincare: 07)
TSx = ( ximi),a7(xickino
i=1 i=l a. 2( x.01,, 1 '=1 when x OiHence! xT's > 0
Vx
0 and the theorem is proved.:Assume that the (real)eigenvalues of M satisfy::
M I
)5 NM 5
5 AN11 5A/
2
.
N
M has a complete 'set of orthogonal eigenvectors and we find (theorem 2.31
(384) A xTx 5 xTMx S AMxTx
In the same way the expression
ST'
(39.)
1SxTxsxTSx
5Xxx
!holds.,when the eigenvalues Xs of S satisfy'
IS 5 N11 5 11
In section 2.2.1 a stable property of the pde was found.
h 1 Since u c
Vo0),
(22) gives 25i=1 =x.4).)
1=1 =d1(uh uh) 4 + (uh uh) = uh) 10,0) dt 2 ' 0 x' x 0 ' 0, or in matrixvector notation: d Ti( uh)TMuh), + Tuh ) = Cuh )TF
It is easy to show, making use of the inequalities of Poincar6
(OM
and Schwarx ((18)) that the stable property is still valid for the discrete formulation.
For the numerical solution of system (26) the eigenvalues of the, matrix M1S
are important. These eigenvalues can be estimated if we use:
Theorem 2.5 The eigenvalues
xS, i
I of1 Q.= MI S
are real, positive. and
:(142 f( 44,)! or
xs
xs Nx9
1 AN _{AN} Sx = AQMx ,Sx =XYGx
ITST ) 1 1 SG GxThis expression shows that the matrix
T T
1(G ) SG
26V i c
I.Proof: M is symmetric and positive definite hence there'exists a
Choleski decomposition (
43)
M =TG.
Consider 1 QM
Sx = A x (41) A(145)
and the matrix
1
M S
have the same eigenvalues. The matrix (GT) 1SG 1 Q is symmetric, so A is real. Moreover if
M1Sx =Qx
TSx =AQT
x x_
Mx 4TSx xT x x Tx xT x MXAQ
From (38) and (39) we find that (42) is valid and the theorem is proved.
2.2.5 The discretization error.
Using functional analytic tools it is possible to prove ([45]), assuming that u and g are smooth enough, that for t fixed,
(46) liu(x,t) uh(x,t)ilo = O(x)
where:
u(x,t) is the solution of
CO
uh(x,t) is the solution of (22)
and
p is a coefficient dependent on the degree it of the polynomials
that are used in the space
In the example treated in 2.2.3 p = 2 when linear polynomials are used
and p = 3 when quadratic polynomials are used to build Vo(2).
We mention the functional analytical approach since it is applicable here. However for more general problems it is difficult or impossible to use functional analysis in order to get error estimates.

28.We shall analyse the behaviour of the error using Taylor expansions and Richardson extrapolations.
In chapter 4 some examples will be given of this type of error
estimation.
2.3 A linear transport equation.
2.3.1 A property of the pde (2).
Multiply both sides of equation (2) by u(x,t) and integrate over = [0,1] to find:
(47) (ut,u)0 + (ux,u)0 = (g,u)0
Evaluation of the second term leads to
Id
(48)7
_{c dt}
(11u1120) u2]10 = (g,u)0Equation (48) shows that the quantity
huh 1
remains constant in time except for the effect of sources, i.e. g, and
the transport through the boundaries i.e. [u2 J.
0'
In section 2.3.4 we shall see that the discretized form of (2) that we
use still possesses this property.
Remark: It can be easily seen that the boundary condition
u(0,t) E 0
leads to a stable problem. If instead of u(0,t) the value u(1,t) is chosen one cannot prove stability.
Using the theory of characteristics one can prove that in that case the problem is not wellposed ([35]).
2.3.2 The Galerkin formulation of (2)
For problem (2) we assume that the space part of the solution u is
29A discrete formulation is found when the standard Galerkin method is
applied using the finite dimensional subspace liE1(c2) of 4(2).
Hence:
Find uh such that
t' 0 x 0 0 V vh E Vi(C2)E
(uh
vh)
+ (uh,vh)
= (g,vh)
1
with uh(t) c
VE(0,
V t c [0,T]and uh(0) given.
N+1 . 1
Define a basis
(0ii=1
in VE(Q).1
Since
ut.)
E V1(Q) we write:uh(x,t) = y u
4.(x)
i=,
Take j = 1,2,...,N+1 in
(51)
and the result is a SODE for theunknowns
u(t)
1 N+1 h d h./
tli(t))(,i4j)0
./ u.(,.)0 = (g,4)0
,j =
1,2,...,N+1. 1=1 1=1Initial conditions for system (52) follow from:
N
h, h
u kx.,0) = u (x.) = u.(0) 4.(x) ,
j
= 1,2,...,N+1.0 j j=1 1 j
2.3.3 The system of ODEs
Equation (52) can be written
Mu.h + Duh = F
with M and D (N+1) x (14+1) matrices with elements:
and m.. = (o.4.) ji j 0 ,do. d.. =
dT4j)0
1For the case that
VE(SO is the space of linear polynomials and when
equidistant nodal points are used then we find:
The vector F is:
1 82
F= M
012
30Remark: Notice the last row of both matrices.
2.3.4
Properties
of system (53)The matrix M is symmetric and positive definite. The matrix D is skewsymmetric, except for the element on the principal diagonal of the last row, due to
the lack of boundary conditions in x = 1. Consider do. do. 1 d.. = (1 O.) _{= (0.dx} 0 ,._..l)[O.(x)0.(x)] ij dc ' 0 j j 0 0 Ax 6 2 0 0 0 0
31Since the {cki(x)} , i c {102,...,N+1} form a basis of the space
V1 ([0 1])
"
not all Oi(x ) satisfy Ethe conditions:
03.(0) = 0i(1) = 0.
For instance, in the example treated in section 2.3.3 the basisfunction
N+1(x) yields 0N+1(1) = 1,
and ass. consequence
dN
= _{+ 1.}
+1,N+1 dN+1,N+1
In section 2.3.1 a conservative property of the pde was derived.
Since
uh c V (2)' (51) gives:
E
h h h h
(ut,u )0 + (ux,u )0 = (g,.. )0
or in matrix vector notation
d11,
h)Tmuh) + (uh)T011h = (uh)TF
ci:ETL1
_
13: partial integration
(uh,uh)0 ..(uh,uhx)0 Nu h)2310
or in matrix vector notation
hT
h= (uh)TDuh + [(uh)2110
It is easy to show, using
(56)
(or (57)) that(54)
(or (55)) expressthat the conservative property of the pde (see section 2.3.2) is still
valid.
For instance in the special case treated in 2.3.3 we find from
(55)
and(57):
(55
+
+.32= d
hTh1h2
h T 1(58) ((1.1 ) 1N+1I = i(u 1 Y 1, .2 1The term kuh ) has a stabilizing effect. Again we see that a boundary
2  N+
condition at x = 1 instead of a condition at x = 0 will lead to the
unstable problem:.
4
) 2 I,thTh
h (59) .(,(a Mu )  ,(u )0 = )For the numerical solution of equation (53) the eigenvalues, XR
_{i *}
i = 1,2,.%.,N+li of the matrix
R = M1D
are important.
These eigenvalues satisfy:
, (60) =
1i
or 1 Dx. = X.Ma. 1 T Define DS = (D+D2 N and 1 T. DT = (DD them T T x.D x. x.D x.1
sa a Ta + T T x.Mx. x.Mx. 1 a1 1
The term T x.D x. S3. T 1 athat represents the real part of the XR is, as showa,, the consequence of
the boundary conditions.
This term can be estimated using theorem 2.5.
The eigenvaluea of M1AT, are estimated in the following way:
1(6N) a
1i
2
h,T M1Dx. I
.331
Theorem
2.6
: The eigenvalues Au of the matrix U = MDT are imaginary and satisfy DT max IA I 1
(62)
max IX2,1 Imin IXM.I9.
Proof: Consider M1Dx =
AuxM is symmetric and positive definite hence there exists a Choleski decomposition M = GTG.
Then:
D x = AGTGx
Or
GT)1DTG1Gx = AGx
This proves that the matrices M1DT and
T1
, 1 (GDTG have the same eigenvalues. The matrix
(GT)1D G1
is skewsymmetric and hence has imaginary eigenvalues. Besides that T T =
xDxxx
U A T T x x x MxUsing theorem
2.3
we prove that the estimation(62)
is valid.2.4
The linear convectiondiffusion equation.2.4.1
A proPerty of the ode (3)We multiply equation
(3)
by u(x,t) and integrate over 12 = CO31] to find,after partial integration and usage of the boundary conditions, 1
'466)
34K 64 4'
Kuu.)
+ n(u + b(u ,u ) =,Kg,u)rt
0x
0x
x 0In the same way as it was proved in section 2%2.1: we prove that the
solution u of
(3)
satisfies: (63), Ild(t)110, e2btru.a. I I)110
1 _{e2b(5t)}II)11
g(s1.0 ds, ?ke t E2.4.2 The Galerkin formulation of (3)
For problem (3) it is proved
in [7]
that the space part of the solutionu for all t c [0,T] is an element of the space
Hence, if we use the finite dimensional subspace
v1(o)
of 11011(0,, thenthe standard Galerkin formulation for (3) is
Find uh such that
hh
,' Of
h h h h h
(u,v)
+ aku v ) + bku (s,v)0vhr v1,2)1.
t 0 x x'
with
1
uh _{c Vo(0,,,} V t c
10,T3
hand
uh40)
ziven
Initial conditions follow from
(67)
uh(x.,OD
u
=I u.(0)
0(x)
V
i=1 1
2.4.3 The system of ODEs.
Problem i,(66) can be written;
.
Mhu + aDuh + bSuh = F
The elements of the IN x N matrices M /1, ID and, S are
m = (0.,..) dhi d  = ji dx JI 0, d0i dC. s,. =
a)
dx 'dx 0, ,u)[0,T],
0 = 0 0 ((70)
The functions {0 ). are a basis of V1(St) and we see immediately that
i lEI
M is symmetric,
is skew symmetric and
is symmetric.
2.4.4 Properties of system
(68)
From (66):
1 h
h)0} + a(nhc,uh)0 + b(uh,uh) = (g uh)
(69)
,u
dt 2xx0
' 0 Since uh(x)]10 = 0 we find a(uh,uh) = a(uh,uh) 0 x 0 x 0Application of the Poincare and Schwarz inequalities shows that
Iluh(t)110 e2btiluh(0)/1 a T x Mx + b =
35Now we can use theorem 2.3 to prove:
0 + f e2b(
0
ts)lig,ssit as%
II
Is(
s)'And once more a property of the solution of the pde is also valid for the solution of the discretized problem.
For the numerical solution of system (68) we investigate the eigenvalues
of Q = M1(aD + bS) Consider = x (aD + bS)x = AQMx T T x Dx x Sx
_{_}
T x Mx S dand
with
D
max
A max = max _{AD.}I
id
I 1A max = max tA.)
id
1M
A min = min (A.)
icI 1
Observations
In sections 2.2  2.4 we have seen that the application of the Galerkin method and the f.e.m. in order to solve partial differential equations of evolution type leads to a system of ordinary differential equations.
For this SODE properties of the original pde concerning stability and/or conservation are still valid.
In this thesis we study partial differential equations with a convective and/or a diffusive part.
Suppose that the SODE, obtained after application of the Galerkin method to such a pde, can be written
(72) MU = Au + F.
It is wellknown that the eigenvalues of the matrix
361
Theorem 2.7 The eigenvalues of the matrix Q = M (aD + bS) belong to
the rectangle given by
0 5 Re(A)
5 a
and B 5 Im(Ai) 5 B,i
E I when a  b M As max A min  a A min . 2.5M1A
are almost imaginary if the convective terms dominate and alLiost real negative if the diffusive terms dominate ([32]).
Hence for the numerical integration of a SODE that results from a discretization of a stable transport type equation a method is needed that preserves stability properties near the imaginary axis.
When a heatdiffusion type equation is discretized the stability properties along the real axis are important.
A remark should be made concerning the accuracy of the numerical solution. The numerical solution of the pde is found after application of two
discretizations: the space part and the time part are discretized
seperately.
Both discretizations introduce an error. It is clear that for the
numerical integration of the SODE a method should be used that introduces an error of approximately the same size as the error introduced by
the discretization of the space part.
37
383 Methods for solving ordinary differential equations.
3.0 Introduction
In chapter 2 we have seen that the application of the f.e.m. to partial
differential equations of evolution type leads to a SODS.
In this chapter we give the properties of a few classes of methods for
solving numerically ordinary differential equations.
Furthermore some methods will be analyzed for stability and accuracy.
During the last few years several books on the numerical solution of
ordinary differential equations have been published, for example
[12], [15], [21], [25], [26], [44]. In all these books the solution
methods for the system
{Y
= f(y,t)Z.(0) =
yo
are mainly generalizations of the solution methods for the equation
= f(y,t)
y(0) =
yo
In the first part of this chapter we follow this approach too. However
in the second part some special methods to solve systems of ODEs will
be considered.
Although many useful implicit methods exist we shall almost completely
restrict ourselves to explicit methods. In chapters 5 and 6 large systems
of ODEs occur that are nonlinear. It is wellknown that nonlinearity causes difficulties if implicit methods are used, i.e. for every time
step a nonlinear system of equations has to be solved.
3.1 Truncation error and stability.
Before we start the analysis of the methods the definitions of order
(truncation error) and stability, that are used in this thesis, are given.
Suppose that the solution y(t) of the system
{i
= f(y,t)
(1>
16(5)
is approximated by the solution. un of the numerical method:;
i'. un+1 nun,tnoltn) with 0 ' = and At  stepsize
Remark: The notational convention that un denotes the numerical approximation of the exact solution z(t) at time tn is adopted.
Assume that the solution x(t) of the SODE is. "smooth enough", then we find using the differential equations, that
zetn+1),
F(z(e),tn,Atn) + dn
n
d is called the local truncation error.
Definition
3.:
The numerical method(3)
is said to be of order pif
the local truncation error satisfies dn = 0(AtP+1)
The definition of stability that is used is influenced by the type of problems that have to be solved. Our ultimate goal is to solve the
shallow water equations. For these equations a strong stability definition is required since the integration has to be carried out over a long
timeinterval and the solution is almost periodic. Therefore we use:
Definition 3.2: The numerical method (3) is said to be stable for a given
stepsize t and a given differential equation (2) if a positive constant K
exists such that the change due to a perturbation of size 6 in one of
the values un is no larger than 1(6 in all subsequent values um, m > 12,
m,n e IN..
This definition is dependent on the differential equation treated, so the idea of a "testequation"' is utilized. K[12])' The testequation is the differential equation
'0.6) = ASp.
AEC., Re00
SThis leads to:
Definition 3.3: The region of stability of the numerical method (3Y
is that set of values At (real, nonnegative) and A for which a perturbation in a single value un will produce a change in subsequent values which
does not increase from step to step.
Suppose that the SODE
= ftz,tj
has to be solved.
Such a system locally behaves like the linearized system
an
X a
_{az}
n )(ZXn)+ 1(in,tnY
+ ; n n(ttri)J.,t=t ,t=t
In such a case we choose the stepsixe At so that for all the eigenvalues
of the Jacobian matrix
the requirement that A. At is in the stability region, is satisfied.
n
3%2 'Three classes of methods for ODEs
In this section we consider briefly The Adamsmethods,
(HY
The RungeKutta methods,The Extrapolation methods.
Often it is desirable that a method offers some flexibility. In particular,
during the integration, one would like to be able to
I Estimate the error.
1:(8 2 Change the stepsize accordingly.
3. Change the order accordingly.
n z=z ,t=t.
_4o, 0.
(9)
With respect to these properties the three classes mentioned above show the following characteristics:
Adamsmethods (P(redictor)C(orrector)type.)
Since it is easy to change the order, it is also easy to estimate the error. However the stepsize can only be changed at the cost of a lot of
organization.
A package that is based on the Adamsformulas is the one described in [44].
ii RungeKutta methods (SK)
During the integration it is difficult to estimate the error. If the classical RK methods are used one can change the order only at the cost of additional function evaluations. As the RK methods are onestep methods it is easy to adjust the stepsize.
In [10], [11], [41] and 1:53] methods are derived that are based on the "RKidea" but which have errorestimators at the cost of little extra
work.
ii Extrapolation methods (Ex)
Compared to (i) and (ii), these methods were developed quite recently ([2 ], [13]). They combine the advantages of the Adamsmethods and the RK methods with respect to the points of (8).
Not only the properties mentioned in
_{(8)}
are important; stability andaccuracy play a role and also the number of times that the function f
(in the equation Sr = f(y,t)) has to be evaluated, especially if this
evaluation is difficult.
3.3 Some examples.
In this section we give the formulae and a picture of the stability regions
of some examples
of
the three classes. Also the order p of each methodis given. (See also [26])
3.3.1 AdamsBashforth 2 (AB2) n+1 n At = 4. _{(3±'(t°)tn)} f(un1,tn1)) ' s u u
_
_
p=
2.41,
423
2
1REAL AXIS
fig. 3.1
Stability region of method
3.3.5.
32
03
2
REAL AXIS
fig. 3.2
Stability region of methods
3.3.2 AdamsBashforth 3 (AB3)
n At n n2 n2
un+1 = u + (23f(u,tn) 
16f(u,tn1)
n1 + 5f(11 ,t ))(io)
P= 3
3.3.3 AdamsBashforthMoulton 2 (PC) un+1,* = un + Atf(un tn), un+1 un At,J1,
f(un+1,*,tn+1)) 2 (1."1 '6 )p= 2
43n+1,** = un + Atf(un2un+'*+2u1141'*,t11+1)3.34
Adams3ashforthMoulton 3 (PC) f(un1,tn1)) (12) n+1,* n At n n= u
+ 7
(3f(u,t)  f(un1,tn1 U )) n+1 un Alt2 (5f(un+1,*,tn+1 11 ) 8f(un,tn) = 3 3.3.5 RK 1 (Euler explicit) n+1 n n nu=i + Atf(ut)
(13) 1 3.3.6 RK 2 (Modified Euler) n+1'* = un + Atf(un,tn)(14) un+1 = un + {f(un,tn) + f(un+1'*,tn+1))
p = 2
3.3.7
RK 3 (Heun method) un+,* un At f(un tn) 2 (15) un+1,* un iltf(un+i,*,tn44)p
At3
2
1REAL AXIS
1
3
2
1
REAL AXIS
fig. 3.4 Stability region
3.3.8, 3.3.12.
^_{3}
0
fig. 3.3 Stability region :f methods
3.33, 3.3.7, 3.3.11.
3
0
of methods
(16)
(17)
n+1 1
T2un+,*+u
_{u.)}
P= 3
3.3.8 _{RK 4 (Classical RungeKutta method)}
un+3,* n At n n
_
_{=} 1. + 1(11 ,t _{)} unig ,** _{= un +} AL f(un+i'*,tnig) 1 un+1,* un + At f(un+i'**,tn+)_
_
un+1 '** = un + Atf(un+1'*,tn+1)_{ } ** n+1un+1 
(2u
1.' +4u
* n+'2'
_{+2u}_{'}
*_{+u}_{n+l}_{'**}
_{3u}_{n)}
6
p=
3.3.9 Trapezoid (Implicit) I n+1 n At u = u + (f(un,tn) + f(un+1,tn+1))2
p = 2 3.3.10 Ex 2 (Extrapolation)(18) n+1 an+4 pt
fie1g
tn+i)= un + AL f(un,tn) = un + Atf(u ,t ) 2
2 ' n+1 n n, n+1 n+1 n+1 U = 2u  up= 2
3.3.11 Ex 3 (Extrapolation) n+1 = un + Atf(11.n,t.n) a un + AL f(un,tn) 2 1 n+I n+; At ,n+ nigu + ,t ) 2
45n+i u462
1REAL AXIS
1
0
fig. 3.5
Stability region
of method3.3.1
3.3.4
2
1
0
fig. 3.6 Stability region
of method3.3.2.
2
1
C20) un+T 1 n+T n+1 n+1
=L
8u +9u ) , pr = 3 3.3.12, Ex 4 (Extrapolation) n At I n nn u=u +
,t )_n+1A ,n+2
= un + t ftun+1n
At , n n = u + 411)UP+1
n At 21
n+4 n+4 At ,.,n+;1;I
3 ri. = u +  ftu  tr1+41 3 i n+11nq
At.4
n 4 ti = il _{+ i f(11} ,t )I n+1 1 (4n+1 n+n, u 3p=
In figures
_{3.4}
to3.6
the stability regions of these !methods are plAtted4RemArk: The stability region of method 3.3.9 is the. negative halfplane, i.e. the region with Re(A) 5 04
Table
3.1
contains the boundaries, of the various stability regions_{on}
the real and imaginary axes together with the number of function evaluations required in one full integration step At.419)
5
1 n 3 U n+g  3 u zin+11 n At , n n= u +
3 f ,t, ) 1 n  At ,n n = u +5
fku ,t 2 2 2 n+ n+At f(u 3,t3)
1 3 + =
48Table 3.1 Stability boundaries on real and imaginary axes and number of function evaluations per step At.
3.4 Systems of ordinary differential equations.
Up to now all the methods that have been considered were in fact generalizations of methods based on the scalar equation
= f(x,t)
Here we shall consider briefly some methods constructed particularly for solving SODEs.
Consider the system of differential equations
Sr = Az + f(t)
z(0)
= yo
t , [to,T]The matrix A can be written as
A =L+D+ U
with L a lower diagonal, U the upper diagonal and D the principal diagonalmatrix. The GaussSeidel method is defined by the formula
Name Formula Reaxis Imaxis Evaluations
AB2 (9) 1.0 .35i 1 AB3 (10)  .55 .701 1 ABM2 (11) 2.0 2 ABM3 (12) 2.4 1.21 2 RK1 (13) 2.0  1 RK2 (14) 2.0  2 RK3 (15) 2.52 RK4 (16) 2.8 1.71 2.81 3 4 Trap ' (17) .., Ex2 (18) 2.0
wi
1 2 Ex3 (19) 2.5 1.7i 4 Ex4 (20) 2.8 i 2.8i 5 ., .(un+1)T(
49n+I n+1
= un + At(L+D)u + AtUun + Atfn.
Remark: (23) may be followed by the formula
un+2 = un+1 + At(U+D)un+2 + AtLun+1 + Atfn+1
We investigate the stability of method (23) for two special cases:
i) Suppose A = AT, hence L =T
Then since f(t) is not of interest in the stability analysis, we consider
(EAt(L+D flun+1 = (E+AtLT)un
with E : the identity matrix.
We multiply (25) by
(u1)T,
tranpose and addr T , n+1
(un+1)2EAt(D+L)QtkD+L ))u
E+AtLT)un+(un)T(E+AtL)un+1
Multiply (25) by (un)T, transpose and add:
(27) _{(un)T(2E+AtLT+AtL)un} =
(un)T(E....6t(134.1,))un+1 (un+1,T. T, n ) kzAt(D+L flu
From (26) and (27) we deduce:
n+1 n,T
(26) ku
u )
,k2EAtD)(un+1u)
n (u ,Tn) (AAt)un = (un+l)T(AtA)(un+1) In the case that A is negative definite, since then2EAtD is positive definite, it follows that
(29) (un+1)T (tA)un+1 (un)T(...6tA)un
So if the norm 11.11A is defined by
(23)
(2/0
426)1
+
fir
00)
xT(A)x
follows from (29) that
11/11/+1
I'
A 1111iln1111A,and unconditional stability has been proven for both methods (231 and (23) followed by (24), provided that the matrix A is symmetric and negative definite.
ii)
Suppose A = AT,, hence If = LT and D = 0.Proceeding along the some lines as in case, CO results
iD:
(un+1)T(2E+At(L+LT))ufl"
nT
) (2E+At(L+LT))2
If the matrix
2E + At(L+LT) =
is positive definite we find from (3011
Ilten1B,
=I ID
Hence the method is stable and conservative. The property of conservation is also present for the analytical solution, for if,
with A =' AT then: dt (xT x)' = 2x Ax
= d
d hence xTx = 11x112 is a constant in time.The order of method i(23) is easily found
by
means' of u Taylorrserles: 5Q,(30) pi= V = it = 3 =
51In the same way as the GaussSeidel method is related to the GaussSeidel iterative procedure, it is possible to relate a Jacobi method to the Jacobiiterative method:
un+1 = un + At(AD)um + AtDun+1 + Atfn
For stability we investigate the homogeneous part:
(EAtD)un+1 = (E+AtAAtD)un
Suppose that A = AT and that A is negative definite: Transpose, multiply and add to obtain
, n+1 T At _{n+1} , n,T At n
(u ) (EAtD+ A)u s (u ) (EAtD+ A)u .
Hence method (32) is stable if the matrix
At
Q = E  AtD + A
2
is positive definite.
The order of method (32) is one: p = 1.
When instead of the linear system (21) a nonlinear system is considered then the implicitness of the methods (23), (23) and (24), and (32)
causes some problems. In that case a nonlinear solver, like the NewtonRaphson method ([22]) has to be used.
Remark: It is possible to construct methods of order p > 1. Compared
to the methods of the preceding sections these methods have the advantage of good stability properties but the disadvantage of being implicit. It is completely dependent on the problem under consideration which of the methods is to be preferred.
3.5 Extended stability regions.
When a SODE resulting of the discretization of a partial differential equation is solved numerically by an explicit method a stability
requirement with respect to At has to be satisfied. Accuracy also imposes restrictions on the stepsize At. In a lot of ordinary differential
equation problems stability and accuracy permit stepsizes At that are of the same order of magnitude.
(34)
{
unfl'*
= u+
n Atf(un tn)52However in some systems (for instance: stiff systems) the stability
condition imposes a far heavier restriction on At than is necessary for the required accuracy.
In such cases implicit methods (for instance method (17)) may be used since they have in general mild stability conditions.
It is also possible to construct methods that have a relatively large
stability region but that are still explicit. Additional function
evaluations, which may be expensive, is the price that has to be paid
for this desirable behaviour.
In this section formulae of three methods, with good stability properties,
are given. (For more methods of this type see [21]).
The first formula is a modified Euler formula. The extra function evaluation per step is used to extend the stability region. The following problem is solved:
Find a and S such that the method
un+1 = un + at(af(un,tn) + Bf(un+10*,tn+1)
has maximal stability along the real axis and is of order p = 1.
Remark: In this example we take the real axis. Of course any other line
in the complex plane could have been taken.
In chapter 2 we have seen that the solution of the problem stated above
will lead to a method that is appropriate for the solution of a BODE resulting from a partial differential equation of diffusion type.
The solution of (34) is (see [33])
7 1
a
In figure 3.7 stability regions of method (34) for some values of a and 5 are plotted.
The second method of this section is designed for the solution of SODEs.
5iI
8

St4
;
/
t
REAL AXIS
Stability
region method 3.(34)I
1for a = 8 , 5 _
3
8
7
6
5
43
2
REAL AXIS
fig. 3.7.1
Stability
region method 3.(34) 3for a
= T ,
5 0 T.3
2
8
;
6
5
4
3  2
; REAL AXISfig. 3..7.2
Stability
region method 3.(34)
54un+1'*
=un
+At(AD)un
+AtDun+1'*
+Atfn
while the formula
n+1 n At n n+1,* n+1,
= u + LAu + Au +
fn
+ f )2
is used as a corrector.
For (35) Pi = 2.
In chapter 4 it will be shown that for a diffusion problem method (35) has a stability boundary which is more than 3 times the stability boundary of method (14).
Methods (34) and (35) are stable if the eigenvalues A of the matrix A lie near the real axis in the negative half plane.
If the eigenvalues are near to the imaginary axis other methods are
preferable. The last method that we consider in this section is a
method that is appropriate for convection type problems.
It is well known that the RK4 method, (16), has good stability properties
along the imaginary axis (see also table 3.1), however it has order
p = 4. Hence we investigate if it is possible to construct a RK type method of order p < 4 that has an extended stability region along
the imaginary axis and that is more economic in use than the RK4 method.
Before we start this investigation we give two definitions:
Definition 3.4: The amplification factor of a numerical method
un+1 = F(untn,At)
is the quantity P(XAt) in the expression
un+1 P(XAt)un=
resulting from application of this numerical method to the test equation
Y= xy
Definition 3.5: The effective time step Ateff of a method, when X E C
is given, is (35)
=
(36
It
At_{eff} = sup
tiR
N '
AT E&
$ is the stability region of this method and NI the number of function evaluations per integration step.
Remark:: The effective time step obviously depends on A..
The investigation is done in the following way, Firstly we try to establish an amplification factor that has maximal stability along
the imaginary axis, given the order of the method and the number of
function evaluations. After that we try to construct a RK method that possesses this amplicication factor.
First we try three functions evaluations with second order accuracy, hence an amplification factor
X37) P(At)1 = I + At + 2At2 +
Aat3
2 a
'Remark: _{If partial differential equations of convective type}
are discretized
in space by the f.e.m. the discretization _{error is at least of order}
p = 2. Hence we investigate second order methods for the time discretization.
Theorem IT. The amplification !actor
_{(37)}
_{is maximally stable. along}the imaginary axis for a = 4.
Proof: Suppose X .=.41
.2 2 _{3} _{3}
IniutIt3 2
= X1
211112.g_41.)2
2 a '
14 1 2
=1i + (PAt) _{}
y2At21
_{.}
2
a
4
Xi)
When _{u} _{p, stability requires}1 2
4  a
hence: 0 < a
< 8
A
55.. 1 2 p2At2 (11) Ta a+ a2 < 2 a 42At2 < 2a
So the value a that maximizes the expression
2
f(a) = 2a
in the interval (01,8) satisfies the
criteriaBut
fi(a) = 0 if a = 4,
which completes the proof.
If we can find a method with the amplification factor (37) with a
= 4
then the stability region along the imaginary axis will be
2 5 11XII At 5 2
Comparing this to the stability region of
(16)
(see table 3.1) showsthat the effective time step of
(16)
along the imaginary axis is longerthan that for such a method.
Hence such a method would not be satisfactory and we therefore investigate
also the amplification factor: (second order accuracy)
2,42
4 4At3 to, At
C30
PULAt)= 1 + AAt + =S= +1424
Note that method
(16)
has an amplification factor of this type (a=b=1)In fact we will show that
(16)
has the optimal amplification factor ofthis type.
Theorem 3.2: The amplification factor (38) is maximally stable along the imaginary axis if a=b=14
Proof: Suppose X = iu 1P(iuAt)12
=2At2
I. h 4At4 2 ) + (aAt /111PA,2 6 m 56 0 + 24lc I
572) u4,st4b24 4 b _{A) + }1.12St2(g
ie
576. 1 * y OM4
1f(r,
71 3
Consider the polynomial:
Q(x), b2
4(a2
th 2 a b 1576 x
4.
'ig  24'  4. 12 L'
A necessary condition for stability is Q(x) s Of, and in particular this
should hold in the neighbourhood of x = 0, hence:
4
5 (3and when equality holds we also require
a2
0.n7Tbco
From the results of table 3.1 we know that, if a=b=1
q(x) s Or x2 s 8
Q(x)1 > 0,,
xa> 8
We next intend, to show that, for all a and b satisfying (i), we have
416/F.1 2: 0 with equality only if a=b=1,.. This will prove the theorem.
On substitution we find, after some rearrangements
Q(4) =
1)2 7_{2(S}
 74
using ti)i leads to
GIVE) t
(b_
2)2,t_/2\2_
_{.}I (1,
1)29 8 '"1,21 iTE 8
and this shows that .41(07g) 0 if and only if a=b=ll.
Theorems 3.1 and 3.2 show that
(16)
'(the classical fourth orderRungeKutta method) is the best one , among the second order methods,
that we tried to extend for stability reasons. Besides that, the method is very accurate: p b b
(i)
(b+ =(p)t
'LK
2 0
58I Ili
2i
3i. 3.8 Dissipation of the R104 method.
1
2i 3i