NUMERICAL SOLUTION
OF THE
SHALLOW WATER EQUATIONS
BY A
FINITE ELEMENT METHOD
SS a V£BVAIL£N !
)Z S a — » — I... ^ . i
«o «o
00 -4NUMERICAL SOLUTION
OF THE
SHALLOW WATER EQUATIONS
BY A
FINITE ELEMENT METHOD
NIEK PRAAGMAN
1376
419
0
U DelftNumerical 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 li+.OO uur door
Nikolaas Praagman
wiskundig doctorandus
geboren te Drachten
Dit proefschrift is goedgekeurd
door de promotor
met dank aan:
Phil en Sue Anderson,
Bert Broere,
Gonny Scheurkogel - v. Werkhoven
en voor al
Contents
0^ Introduction • ' 1
^_ The shallow water equations 3
1 .0 introduction 3 1•1 the system of partial differential equations 3
2 T h e finite element method to solve partial differential equations 11*
2.0 introduction l i t
2.1 mathematical tools ]h
2.2 the linear heat equation I9
2.3 a linear transport equation 28
2.k the linear convection diffusion equation ' 33
2.5 observations 36
3^ Methods for solving ordinary differential equations 38
3.0 introduction 38 3.1 truncation error and stability 38
3.2 three classes of methods for ODEs hO
3.3 some examples 1*1 3.*t systems of ordinary differential equations 1(8
3-5 extended stability regions 51
3.6 observation 59
U_ Stability and numerical results 60
U.0 introduction 60 U.1 methods to investigate stability 60
U.2 lumping 68 U.3 results for system 2.(26) 69
h.k results for system 2.(53) 72
14.5 results for system 2.(68) 82
^ One-dimensional incompressible flow 8U
5.0 introduction 8U 5.1 boundary conditions 81* 5.2 the Galerkin formulation 87
5-3 the system of ODEs 88 5.^ numerical stability 89 5.5 computations and results gO 5.6 remarks (other methods) 95
5.7 conclusions 96
A numerical treatment of two dimensional shallow water flow 6.0 introduction
97 97
6.1 boundary conditions 9j 6.2 the Galerkin formulation 101 6.3 the system of ODEs 102
6.1* stability IO3 6.5 computations and results I05
J Conclusions and remarks 13I
Appendix 1 133
Appendix 2 I38
References 11^3
T U Delft
Technische Universiteit Delft
Memorandum
'K
•yu
'i
1_l
I^^V-vVA^
t
dit proefschrift werd gedrukt bij: drukkerij & uitgeverij
EGNER bv Den Helder
Introduction.
Mathematical models for two dimensional flow have been developed at various locations.
For instajice, in each country that borders the North Sea, an operational numerical model exists, (see [1], [H*], [l6], [17], [27], [37], [^7]).
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 ntimerical 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 (RK^). 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 non-steady motion of long-waves.
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. (SODEs).
Chapter h: 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 1-dimensional 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 U, 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. ([U2], [1*3]).
-3-The shallow water equations.
1.0 Introduction
1.1
(1)
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 Eulerian coordinate system.
In the final system of equations coefficients occur that are not completely determined by theoretical considerations. Hence in the ntimerical computations empirical values have to be used.
The system of partial differential equations.
Consider a coastal sea and define cartesian coordinates x and y in the "horizontal"-plane of the undisturbed water surface. The
z-axis is chosen perpendicular to X-, and y-axis. (See fig.1.1 ) The general form of the continuity equation for a fluid reads:
y-axis
x-axis
fig. 1.1
In (1) u, V 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)
(2)
h ^p"^ * fc ^p"^^ * fc ^"''''^ ^ h ^^"""^ = p^x
;3)
f^ (pv) + 1^ (puv) + 1^ (pv^) + I Y (pw) = pFy
(M
It (pw)
*
li (P^w) + I7 (PVw) +
J^
(pw^) K pF^.
In 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 Coriolis-force - Gravity
- The pressure force - Friction forces - Wind forces
The contribution of these forces to F^, F and F^ is:
Coriolis: F^ : = 2«3V - 2«2"
;5) Fy : = 2n^w -
2a^u
(7)
F : = 2n„u - 2n V z 2 1 withfl = (n ,np,n^): the earth rotation vector.
Gravity: F : = 0 X
(6) ^y = = °
F^ : = -g withg : the acceleration of gravity.
Pressure ^ x •• . 'y--\ •• with P :
1 iE
p 3x1 iE
^ " P 3y_ _ 1
l2
~ " p 3z pressureFriction (molecular viscosity) F : = vAu (8) F : = vAv F^ : = vAw z with
A : =
^
. 3 _ . l l
-2 -2 -2 3x 3y 3zV : = kinematic viscosity coefficient
Remark: Wind-forces 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), (3) and (1+) leads to:
(9) -^ (pu) + 1^ (pu^) + 1^ (puv) + 1^ (puw) = f,{2n^^-2a^%,) - 1^ + PvAu
(10) 1^ (pv) + 1^ (puv) + 1^ (pv^) + 1^ ,pvw) = P{29.^V-2Ü^M) - | ^ + PvAv
(11) 1^ (pw) + 1^ (puw) + 1^ (pw) + j ^ (pw^) = p(2S22U-2n^v) - pg - |£ + pvAw
Since we are only interested in the tidal motion, we apply, to equations (l), (9), (10) and (11), a moving time-average to separate the short term turbulent fluctuations from the longterm tidal oscillations. (130], [l8], 0*9]). So let the velocity field (u,v,w) be separated into a mean part (u,v,w) and a fluctuating part (u' ,v' ,w' ) in such a way that
t+T (12) - ƒ u dT = u t and t+T (13) - ƒ u' dT = 0 t '
Here we assume that it is possible to make the time averaging interval T 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.
-6-Before we apply this time-averaging technique to equations ( 1 ) , ( 9 ) , (10) and (11) we make the assumption that the density p 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 asstaaption is reasonable if the sea under consideration is well-mixed.
Results of the time averaging procedure are: - .
Equation (1) 3x dy 3z Equation (9) ^ '
'15)
I T
+
h (^^)
+ T7
(^^)
+
h (^^)
=
2«,^
-
2n.w
- ^ | £ +
VAG
+
dt dx oy dZ 5 i^ p oX "t+T "t+T t+T ' - T ^ f c ^ ("'^^'i^+l^ / ( u ' v ' ) d T + | ^ ƒ (u'w')dT} Xi X/ % Equation (lO) ' ' - • 16) I T - T- (üv) + I- (v^) + |- ( w ) = 2n,w - 2n,ü - 1 -^ + vAv + dt dx dy dz 1 3 P By t+T t+T t+T - ^ { | ^ ƒ ( u ' v ' ) d T + | ^ I ( v ' ) ^ d T + | ^ J (v'w')dT} t t t Equation (11) , ' • ' (17) 1^ + |- (üw) + |- (vw) + — (w^) = 2n„ü - 2n.v _ g - ^ |£ + vAw + 3t dx dy oz c: 1 p dz t+T t+T t+T - i { | ^ ƒ ( u ' w ' ) d T + | ^ ƒ (v'w')dT+|^ ƒ ( W ) 2 d T }The last term in each of the three averaged equations results from the non-linearity of equations (2), (3) and (U).
In view of their similar effect on the nodel, it is generally accepted to express these terms as molecular diffusion.
-7-(^ö) I - fc (ü') - 1 ^ (üv) + f^ (üw) = 2«3V - 2«2W - I f ^ + (v+v^) Aü
(^9) If - fc (-) - 1^ (v^) - fc (^) = 2«^w - 2^3^ - 1 0 + (v+v^) Av
The problems of the turbulence phenomenon have not been eliminated since the magnitude of v is not known.
In some models v is used as a control parameter. In other models \>„ is given a value in accordance \rith 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
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 L and L are x y much larger than the vertical length scale L .
If the horizontal velocities are characterized by U then it follows from ' equation (1^*) that the characteristic vertical velocity W will be of the order (22) i ^ ( i f L = L ^ = L^). Now L << L z 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 = 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: ([3CQ)
L U^
'23) ' I T * |- (^^) ^ |- (^) ^ |- (^^) ~ ^ m s-2
3t 3x 3y 3z .2 - -1+ -2 (2U) 2ü^u - 2n V 5s 10 m s (25) - g K - 10 m s"2 (26)1^^ ,
p 3z ~ • _. L U (27) v„Aw « 1 0 - 2 - — m s" T L ^2 , z fi L U , „ (28) vAw = lO"'' - ~ \ ^ ^'^
ir
zIte order of the values that should be substituted for ü , ü and v is well-known. So these coefficients cause no problems. However the value of v„ is unknown and even the order is not quite clear. In [30] it is proposed to take
(29) v^ = 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
(30) U =; 1 m s"^
we see that the pressure term is the only one that can balance the gravity term g. The so-called hydro-static distribution of the pressure, equation
(21), is the approximation that results from this analysis.
This hydrostatic approximation provides a possibility to eliminate the pressure from equations,(l8) and (19). First we define the waterdepth h and the waterlevel E,:
h - The depth h is the distance between the x-y plane and the bottom. 5 - The waterlevel C is the distance between the x-y plane and the
water surface at the time considered.
Let p be the atmospheric pressure then integration of (21) gives:
-9-The variations of the atmospheric pressure are small compared to the variations of pgS, so we take:
and
Equations ( 11*) , (18) and (19) will also be integrated over the vertical. Before carrying this out we define:
~ 1 ^ -(3it) u = - ƒ u dz -h ~ 1 ^ -(35) V = - / V dz - h _ , • with H = h + 5.
It is also assumed that the surface of the fluid is a streeiiline. This means that any particle which is on the surface remains on it. ([^6]).
The interface sea air is given by the equation -• •'
(36) 5(x,y,t) = z
A particle will remain on this surface if at z = ^:
(37) |f--F*^|^=w
3t 3x 3y , - , , The equation of the bottom
-h(x,y,t) = z
provides, in the same way, in the relation at z = -h: ''
,,a\ 3h ^ - 3h ^ - 3h 3t 3x 3y
The integration formula (n stands for x, y or t )
(39) J^ If dz=f- j S d z - f ( £ ) | J - f(-h)|J
_h -h
i s also used.
-10-Integration of
{^h)
leads to:
ƒ 1^ dz + ƒ 1^ dz + w(x,y,£,t) - w(x,y.-h,t) = 0
-h -h
Combination of this relation with (37), (38) and (39) gives:
ff + ü(x,y,5,t) If + v(x,y,5,t) g ^ H + ü(x,y,-h,t) |^ + v(x,y,-h,t) |^ +
| f + ü ( x , y , - h , t ) g + v ( x , y , - h , t ) g +
-h '^ -h '^
(^0) =lf^l7(«^)^l7^"^' = °
If we neglect the slow variation in time of the bottom, i.e. if we take
we find
c*^) |{^fc(»-)^i(«^) =
°-Integration of the non-linear equations (18) and (19) poses more problems.
The various terms of (l8) give:
(U2) ƒ
^^^=~(!
iï dz) - ü(x,y,5,t) If - ü(x,y,-h,t) | J =
-h -h
= 1 ^ (HÜ) - a(x,y,S,t) II - ü(x,y,-h,t) |^
3 ^ -2
The term
-r~ { j
ti dz) is evaluated
^"^ -h
taking
u = u + u
then, since
ƒ Ü dz = 0 we find:
-h
(1*3) J — (u ) dz = — ( J u d z ) - u (x,y,5,t) T ^ r^ 9 ^ - 2 ^ ^ 9 ^ f^ ~ 2 ^ , - 2 , . , N at, -h 3x 3x 3x 3h _!_ 3_ 3x 3x - u^(x,y,-h,t) i r * T T ( ƒ ö ) dz = -h
= 1^ (Hu2) -
Ü^U,j,^,t)
If - ü2(x,y,-h,t) g + I7 ( ƒ Ü^ dz).
-h
(1»U)
In the same way
5 .
3y u v ) dz = | - (Huv) - uv(x,y,C,t) f ^ + 8y 3y
- uv(x,y,-h,t) T - + T — ( J üv dzl 3y 3y _^
and
(1*5) ƒ r — (uw) dz = uw(x,y,£,t) -uw(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 << U , so w << Ü
In the Coriolis force w e neglect therefore the term 2n w.
Besides that in the following we write f instead of 2ü in accordance with common p r a c t i c e . W e study the tidal motion in a relatively small sea and as a consequence f and g are assumed to be constants.
Integration of the right-hand side terms results in:
(1*6) ƒ f v d z = f H v -h (1»7) - h
11
3x dz = H 3x a n d : e 5 . 2 . 2 -ƒ (v+v ) AG dz = -ƒ ( v + v ^ ) ( ^ + ^ ) dz -h -h ^ 3x 3y * ( v > ) ( f t ) z=£ (y+v ) ( i ü ) ^ T ' ^ 3 Z ^ z=-h
-12-Combination of (U6), (1*7), (1*8), -(37) and (38) gives:
'M)
| ^ ( H u ) + | ^ ( H u 2 ) + | ^ ( H u v ) = = f H V - g If H + (v+vj 1^ ( ^ + ^ ) dz 3x T i \ 2 ^2' -h 3x 3y ,3u, fc (_ƒ' a" dz) - 1 ^ (_/' ur dz) + (v+v,)(f^) z=5 (v+v ) ( — ) ^ T ' ^ 3 Z ^ 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 (1*9) are used for wind and bottom influences, Suppose that
(u,7)
is a two-dimensional wind velocity vector near to the air-water interface. Following [30] we put
:50) (v+v, 3u 'T' 3z
K Ü(ü2+v2)5 z=5
The quantity K, which w i l l be one of the empirical values in t h e final system of equations, i s a function of various v a r i a b l e s , for instance u, V and H.
(51)
In the same way we put, following [9]> [27]
3Ü (^^^T) 3Z 2. ~2 J _ g u(u +v ) 2 z=-h C
In this formula the quantity C, which is called the De Chlzy value, still has to be specified.
S 2- 2- 5 5 H a Au = (v+v^) ƒ ( ^ + ^ ) dz - |- ( ƒ ^2 dz) - |- ( ƒ Gv dz),
^ -h 3x2 g^2 8x _^ 3y _j^
(so the coëfficiënt a represents all the diffusion-like processes:
molecular diffusion, turbulence and vertical fluctuations,) then we obtain
|- (Hu) + |- (Hu2) + f- (Hw) = fHv - gH 1^ + HaAu + Ku(u2+v2)^ +
Ob dX dy dx
gu(u tv'^)^ -^v, A 3 . 9 - ^^ ^ , with A = — ^ + — ^
C 3x 3y
Using equation (l*0) and dividing by H then:
•co^ 3u ^ ~ 3u _^ ~ 3u ~ ^ 35 ^ ü(u2+7^)^ ~ K5(S2+v2)g ^* ^"^ ^y ^"^ c2(h+5) (h+5)
The same procedure applied to equation (19) gives:
•c:-,\ 3v ^ ~ 3v _^ ~ 3v ^ ~ _^ 3^ _^ v(ü2+v2)2 Kv(u2+v2)5
^^
^"^ ^^ ^^
c2(h+Ü (h+O
Equations (1*1), (52) and (53), the so-called "shallow-water equations", \ together with appropriate initial and boundary values constitute a mathematical 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.
2_ The finite element method to solve partial differential equations.
2.0 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. ([1*5]).
The partial differential equations, with initial and boundary values, that will be treated are respectively;
(1) The linear heat equation
(1) u = u(x,t) g = g(x,t) Notation: u = — u t d t X £ fi = [0,1] t e [0,T].
(2) A linear transport equation.
(2) ^t * ^x = 8 u(x,0) = U Q ( X ) u{0,t) = fQ(t) X £ n = [0,1] t £ [0,T]
(3)
(3) A linear convection diffusion equation u^ + au - bu = g t X XX ° u(x,0) = U Q ( X ) u(0,t) = fQ(t)
Lu(l,t) = f^(t)
a,b > 0 x t n = [0,1] t £ [0,T]Existence and uniqueness of the solution are assumed for each of the
problems (I), (2) and (3). Proofs can be found in the textbooks ([6], [7] ,
[28]).
2.1 Mathematical tools
In sections 2.2 - 2.1* 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.
The 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: N
(SI c IR , with boundary 3ü)
The space L2(f2): L2(n) = {<ti(x)| (((1,(1))^ < "} with (u,v)|^ = ƒ uv dx the space H ( n ) : H ^ n ) = {<t,(x)| If., V<l, e L2(n)} with
V. = (il_ li_)
* ^ax^ '•••• 3xj,^ the space H ( n ) : H ^ n ) = {())(x)| <t, £ H \ n ) , (t,(x)| .- = 0} O I X C Oit the space H (fj): H^(n) = {4,(x)| <!> £ H \ n ) , •(x)l ^g = 0} with E c an.It is well-known that L„(n) is a Hilbert space when the inner product (5)
is used.
The spaces H , H and H^ are called Sobolev spaces. They are also Hilbert
;6) (u,v) = ƒ (uv + Vu.Vv) dx
= ( U , V ) Q + ( V U , V V ) Q
Two norms, related with (5) and (6) are respectively:
l|u|lo = ( U , U ) Q and
|u|I^ = (u,u)^.
One can prove that boundary conditions containing derivatives do not make sense for the spaces H , H and H^, since these boundary conditions are unstable ([2l*], [1*5]). Such boundary conditions are called "natural".
The spaces of admissible functions will be discretized using the f.e.m. The discretization of H , H and H„ will be called respectively V , V
' o E -^ •' ' o and V„.
E
The space V is defined in the following way
First the region R is divided into finite elements e., which is called a triangulation T, of H ([!*]).
V (n) is the space of functions satisfying:
(i) The restriction of (j)(x) to an element e. is a polynomial of degree k J
k £ IN
(ii) The function ())(x) is continuous over fi.
The spaces V (fi) and V (n) are subspaces of V ( n ) . The elements of these o E
spaces satisfy respectively
and
i(x) = 0, X £ 3fi if (j>(x) £ V (n)
(()(x) = 0, X £ E if (t>(x) £ Vg(fi)
Frequently vectors x in the complex space C are considered. For these vectors the «.„-norm:
-17-(7) l|x||^=i^x
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 (Oerschgorin): 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
j^i '^
Theorem 2.2 A Hermitian matrix A has a complete set of orthogonal eigen-vectors.
Proofs of these two theorems can be found in [52].
C C C Theorem 2.3 Let C be a real matrix with eigenvalues X^, X_,..., \^.
Let
(8) S = J(C+C^) and
(9) T = ^(C-C^)
When the eigenvalues of S and T are respectively
X^ , i £ I = {1,2 N} and T x: , i £ I 1 then V j £ I : S C s (i) min X. £ Re(X.) £ max X.
iel ^ '^ iel ^
(ii) - max |xT| < Im(X^) < max |x'^| i£l ^ -^ i£l ^
- I t
Proof: The matrix S is symmetric and the matrix T is skew-symmetric. S and iT are Hermitian, hence both have a complete set of orthogonal eigenvectors.
S S S
Let 21i «io>• • • »2^ ^^ ^^^ normalized eigenvectors of S. Then every vector x can be written as
N
x = y a. x.
— .'", 1 — 1 1=1
Using the orthogonal properties of the set
{x.} . -r we find — 1 l£l
-T^
I
X^ la.|2
X Sx > 1 ' 1'-T
"
~1
~
i=l ' (10)A real symmetrix matrix has real-valued eigenvalues, hence:
(i) -T„ X Sx -T X x is real -T g X Sx g ii) min X. < —^— < max X.
, . 1 -T . 1 ^ 1 X x 1
In the same way we find for the skew symmetric matrix T
:ii)
i
( i )( i i )
-T„ x Tx X X -T X Tx IS imaginary -T X X s majc I X . I i£l 'Remark: When X is an eigenvalue of a skew-symmetric matrix, then X = - X
is an eigenvalue too. Hence (11) estimates all the eigenvalues of T.
t h e n :
--|9-x^Cx x'^{nc+C^) + ^(C-c'^)}x
and -T X X -T X Sx ^ = -T X X + -T X X - T „ X Tx -T X XThe theorem now follows from (10) and ( l l ) S T
Remark: For t h e e s t i m a t i o n of t h e X. and X. . i e l , theorem 2.1 can be
1 1 '
utilized.
2.2 The linear heat equation.
2.2.1 A property of the p.d.e. (l)
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
non-homogeneous boundary conditions. , Define w = u - V then w satisfies ' ' w ^ = w + g - v . + v t XX t XX (12) -S w(x,0) = UQ(X) - v(x,0) w(0,t) = 0 ^ w(l ,t) = 0 X £ fi = [0,1] t £ [0,T]
So non-homogeneous boundary conditions result in a different forcing *
function g :
(13) e = S - ^ + v^^
).
In order to derive the property of the solution u(x,t) of problem (l)
we multiply both sides of equation (l) with the function u(x,t), and integrate the result over the domain
-20-n = [0,1]:
(^'*) (^'^)o = (\x'")o " (6'")o
Partial integration gives:
(^5) (^.^)o^ ('^x'\)o = ^\'^]o* (S'^)o
since the boundary conditions are homogeneous we find:
(^6) ^|^((u,u),) + (u^,u^)^=(g,u)^
Using the Poincare-inequality
(17) (u,u)„. ^(^'\)o
and the Schwarz-inequality (18) (u,g)Q < (u,u)g (g,g)Q
we find
(19) f i r (''''')o •"
2(U,U)Q< (u,u)^ (g,g)^
Hence
(20) I|_(||,||2),2||u||2^. ||u||„ ||g||^
2t
Divide (20) by ||u| L and multiply with e then:
d ^ 2 t I I II V 2 t I I II
^ (e M U M Q ) < e llgll^
Integration of this inequality leads to
(21) ||u(t)||Q < e-2^||u(0)||Q + ] e^('-*)||g(s)||Q ds V t £ [ 0 , T ] .
Equation (21) expresses that, in some sense, problem (l) is stable. In section 2.2.1* we shall see that this property is also valid for the
solution 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).
In [7] it is shown that the solution u of problem (1) is an element of the infinite dimensional space H (n) for all t e [0,T].
o
Since V (fj) is an internal approximation of H (n), ([!*]), we approximate u by the solution u of the problem: (standard Galerkin approach, see [1*5]):
Find u such that
(<'^')0 * (-x'^')o = (^'^')0 ^ -' ^ ^ > ) [22) with u^(t) e v\u) , V t £ [0,T] and \ u (0) given. N 1 Define a basis {<|i.}-_, in V (0)
(23) Since u £ V (fl), we have: N
I
i=l u''(x,t) = I u j ( t ) <t,.(x) Take , j «^ Iin ( 2 2 ) . This leads to the SODE with the unknown functions u . ( t ) , i £ I:
(21*)
I (JTu;;(t))(*.,*.)^ + I ui ( ^ , ^ ) ^ = (g,*.)^ V j £ I.
1=1 i=1
Initial conditions for system (2l*) are found from the equation
u^(x.,0) = u j x ) = I u ^ O ) *.(x ) V j £ I
1=1
Remark: System (2l*) is the basis for all numerical computations. Another way to derive this system is to use the classical Galerkin-weighted residual approach:
" h
Let u (x,t) = I u. (t) (|)^(x) i=l
-22-Define: , P N duV N d •.
R(x,t,u^) = I i r * i ( ^ ) l ^ ^ s
-i= 1 -i= 1 dx From (R,4.J)Q = 0 V j £ Iwe find, after partial integration and usage of the boundary conditions, N du^(t) N d(ti. d*.
.1 -IT-
(*i'*j)o
^ .1 \
(ir'l/)o = («'*ó)o '
^ ^ ^
^-1=1 " 1=1
In the technical literature one often uses formulation (22) with v
replaced by 6u.
2.2.3 The system of ODEs.
Equations (25) can be written:
(26) M u^ + S u"" = F (• = T T )
— — Qt
f
M and S are N x N matrices with elements:
(27) -ji=(*i'*j)o
^^ d*. d*.
(28) ^ji = ( i r ' ^ ) o
while
(29) u^(t) = [u'^(t),u5(t) uj(t)]^ and
(30) F^ = (g,*j)o
The evaluation of F. can be difficult for arbitrary g(x,t), therefore we J
use the approximation N+1
(31) i =
I
g,(t) * (x)
k=0 ^ "Example
Consider, for simplicity, the space V (f!) of linear elements, i.e. the space of functions (t)(x) s H (n) that are linear over each element e. (see section 2.2.1)
For the triangulation of il = [0,1]we use an equidistant stepsize Ax and define X. = lAx 1 e. = C x . _ ^ , x . ] , N 1
A basis {*.}. . o f V ([0,1]) is found if w e define the functions 1 1=1 o
()>^(x) by
(i) ({i-(x) is linear over each element e.
(ii) (t).(x.) = 5.- (Kronecker delta) 1 J 1J
Remark: Since V ([0,1]) c H \ [ 0 , 1 ] ) , the è. (x) c uUn) .
o o 1 o Therefore t h e functions
iji_(x) and (j) (x) do not belong to t h e space
V ^ ( [ 0 , 1 ] ) .
For this special case the matrices M and S become (see also [2l*]» [32]> [1*5] and appendix 2 ) (31*) (35) M - M M - - g and S =
Ax
2 -1 -1 2 -1 -1 2 -1 -1 2 while F is given by-21*-(36) F = M Ax Ax 6 Ax ^•N =N+1
V i
0 0 0Vi
Remark: The last two vectors of (36) result from the boundary conditions. Both vectors contain only boundary values.
2.2.1* A property of the SODE (26)
The elements of the matrices M and S are respectively
m.. = (*j,*i)o , i,j ^ I
and d((i. difi.
"ij = (ïï^'dr)o ' i'J ^ I
From these expressions it is immediately clear that both M and S are symmetric. Besides that we have:
Theorem 2.1* The matrices M and S are positive definite.
Proof: Consider for arbitrary x:
N N N N
x V = i X. ( ^ ^ 1 ^ 1 ) = I ^,-( I (<t>i(x),*.(x)) X :
i=1 j=1 "^ ^ i=1 j = l J " J= ( I x.*.(x), I X *
( X ) ) Q =I I I x.4.^
1=1 j=l 1=1(x)l
'i' •''0Since the functions ((i-(x) form a basis of v'(n) they are linearly independent, hence
N
I x ((>.(x) ^ 0
i=l
and we conclude that
m
x M x > 0 V x # 0 .
when X ^ 0
2 5 -N -N s . - , J J x s x = y x . ( y S . . X . - - i = i ^ j = i ^-1 ' N N d()i. d((>.
= J, ^i(J(ir'ïï^)o ^j) =
1 = 1 J = 1 N d*. N d(|).^ J, ""i d ^ ' . S ''j dx )o
1=1 j = l N d(|.. -- I I Y J--l | 2 ' I I / , ' ' i dx I ' o 1 = 1 By P o i n c a r é ( 1 7 )-X Sx = ( ^ ( I x . * . ) , 5 - ( y x . ^ ) ) _
— — d x . ' - , i i d x > , 1 1 O 1=1 " 1=1 N N . 2( y X . * , , y X . * . ) 1=1 1=1 T = 2 x Mx when x ^ O THence x Sx > O ^ ]L ^ ^ s^^d the theorem is proved. M
Assume that the (real) eigenvalues X . of M satisfy:
x < x < < x < x
"l - ^2 - • K-1 ~ NM has a complete set of orthogonal eigenvectors and we find (theorem 2.3) T T M T
(38) X^x,x ^ xMx,< X'x X
In the same way the expression
S T T S T
(39) ^1X X ^ X Sx £ X^x X holds, 3
when the eigenvalues X. of S satisfy
In section 2.2.1 a stable property of the pde was found. Since u £ V (n), (22) gives
-26-C'^)
i-t
^i(-'.u')o> ^ (^'^)o = (^.-')o
or in matrix-vector notation:
n ,\ d ,1, h,T,^ h, , h.T^ h , h.T^ (1*1) -r- {-(u ) Mu } + (u ) Su = (u ) F
d t 2 — — — — — —
It is easy to show, making use of the inequalities of Poincaré ((17)) and Schwarx ((l8)) that the stable property is still valid for the discrete formulation.
For the numerical solution of system (26) the eigenvalues of the matrix M S are important. These eigenvalues can be estimated if we use:
Theorem 2.5 The eigenvalues X., i e I of Q = M " ^ S
are real, positive and
~x^ x ^
(U2) ^ , x Q , ^ V i £ l .
M 1 . M ^1 '^N
Proof: M is symmetric and positive definite hence there exists a Choleski decomposition (1*3) M = G'^G. Consider
(1*1*) M " ^ S X = xSc
or
S X = \%x -f Q T S X = X^G Gx •* T -1 -1 0 (G ) SG Gx = X X x .This expression shows that the matrix
and the matrix
M ^S
have the same eigenvalues. The matrix (GT)-^SG-^ is symmetrie, so X is real. Moreover if M ^Sx = x'^x T Q T X Sx = X X Mx ->• T T x Sx X X (1*5) X X X Mx
From (38) and (39) we find that (1*2) is valid and the theorem is proved.
2.2.5 The discretization error.
Using functional analytic tools it is possible to prove ([1*5]), assuming
that u and g are smooth enough, that for t fixed, /
(1*6) ||u(x,t) - u^(x,t)||g = O(AxP)
where:
u(x,t) is the solution of (1)
and
u (x,t) is the solution of (22)
p is a coefficient dependent on the degree k of the polynomials that are used in the space V [U].
o
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 V ( n ) .
We mention the functional analytical approach since it is applicable here. However for more general problems it is difficult or impossible to use
We shall analyse the behaviour of the error using Taylor expansions and Richardson extrapolations.
In chapter 1* some examples will be given of this type of error estimation.
A linear transport equation.
A property of the pde ( 2 ) .
Multiply both sides of equation (2) by u(x,t) and integrate over a = [0,1] to find:
(^t'^)o ^ (^x'^)o = (^'^)o
Evaluation of the second term leads to
1 d , I I I |2> r1 2-,l , .
2 i r (ll^llo) * ^ i " ^ = (S'^)o
Equation (1*8) shows t h a t t h e q u a n t i t y
ll-llo
remains constant in time except for the effect of sources, i.e. g, and 2 1
the transport through the boundaries i.e. [ju ] . '
In section 2.3.1* 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) i 0
leads to a stable problem. If instead of u(0,t) the value u(l,t) is chosen one cannot prove stability.
Using the theory of characteristics one can prove that in that case the
problem is not well-posed ([35]).
The Galerkin formulation of (2)
For problem (2) we assume that the space part of the solution u is
A discrete formulation is found when the standard Galerkin method is applied using the finite dimensional subspace V (fi) of H (n).
Hence:
: 5 i )
' Find u such that
(\'^ )o ^ (\'^ )o = (S'^ )o
^ ^ '
^E(")with u'^(t) £ vj,(n), V t £ [0,T]
and u (O) given.
r, ^ • , • r , 1 W + 1
Define a basis ii-i- ,
1 1 = 1 in M^iü). Since u (t) t V (fj) we write: N
I
i=1 u^(x,t) = I uj(t) *.(x)Take v =((>., j = 1 ,2,. . . ,N+1 in ( 51 ) and the result is a SODE for the unknowns u'?(t) : 1 ( 5 2 ) N+1 ,d h N ^ d(j>. h, ^1
y (fr u ^ ( t ) ) ( * . , * J o ^ y u ^ ( ^ , * . ) o = (g,*^)^ ,
J=
I , 2 , . . . , N + I , 1=1 i=1Initial conditions for system (52) follow from:
h " h
u"(x.,0) = u (x ) = y u (0) *.(x.) , j = 1,2,...,N+1.
J ^ J n = 1 - ' - J
2.3.3 The system of ODEs
Equation (52) can be written
: 5 3 ) Mu*^ + Du^ = F
m. . = (4 . ,4 . )„ ji '^i'^j'O
Ji
= (^i
dr'*j)o
For the case that V (n) is the space of linear polynomials and when E
equidistant nodal points are used then we find:
M = ^ 1* 1 1 1. 1 1 1* 1 1 2 and D = 2 0 1 -1 0 1 -1 0 1 -1 1
Remark: Notice the last row of both matrices.
The vector F is:
S' 'N+l Ax
6
0 Ax 6 J 0 0 2.3.1* Properties of system (53)The matrix M is symmetric and positive definite. The matrix D is skew-symmetric, except for the element on the principal diagonal of the last row, due to
the lack of boundary conditions in x = 1.
Consider
d(ti. d * . .
Since the {<(i.(x)} , i c { 1 ,2,. . . ,N+1} form a basis of the space
V„([0,1]), not all (|i.(x ) satisfy
E 1
the conditions:
*i(0) = 41.(1) = 0.
For instance, in the example treated in section 2.3.3 the basisfunction
V i ( ^ ) y^^^"^^ V i ( ^ ) = ^'
and as a consequence
d = - d + 1 . N+l,N+l N+1,N+1
In section 2.3.1 a conservative property of the pde was derived. Since
u*^ £ V^(n), (51) gives:
(51*) - (u^,u )^ + (u^,u ) ^ = (g,u )^
or in matrix vector notation .
-i^^\ d ,1, h>T, h, , h,T„ h , hvT^
(55) 1T(^('^ ) Mu ) + (U ) D U = (U ) F
d t 2 — — — — — —
By partial integration
(56) ( U ; ; , U ^ ) Q = - ( U \ U ; ; ) ^ + [ ( U ^ ) 2 ] ;
or in matrix vector notation
(57) (u ) DU = -(u ) Du + [(u ) ] Q
It is easy to show, using (56) (or (57)) that (5'*) (or (55)) express that 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
-32-/^a\ • <i ^z h^T„ h< ^ 1 , h^2 , h>T^ (58) •^ ((u ) % ) + 2 (- )N+1 " ( Ü ) ^
The term —(u ) has a stabilizing effect. Again we see that a boundary condition at x = 1 instead of a condition at x = 0 will lead to the unstable problem:
(59) -^ (-^(u ) Mu ) - -(u ) Q = (u ) F.
For the numerical solution of equation (53) the eigenvalues X. , i = 1,2,...,N+1 of the matrix
R = M"^D
are important.
These eigenvalues satisfy:
(60) M~^Dx. = X^x. i £ I — 1 1—1 Dx. = X^Mx. — 1 1 — 1 ; 6 i ) 1 T Define D = ^(D+D ) and
°T
then: X« = 1 The term = ^(D-D^) x.D X. x.D^. -1 S-L -1 T-i -T -T x.Mx. x.Mx. — 1 — 1 — 1 — 1 x.D„x. —1 S—1 -T x.Mx 9 L pthat represents the real part of the X. is, as shown, the consequence of the boundary conditions.
This term can be estimated using theorem 2.5.
Theorem 2.6 : The eigenvalues X. of the matrix U = M D are imaginary (62) and s a t i s max i min j fy
lA,
u^i
-
T
K
Proof: Consider M ^ D ^ = x"xM is symmetric and positive definite hence there exists a Choleski decom-. decom-. T position M = G G. T h e n : - • - ' D ^ = XG'^GX or • (G'^)"^D G"^GX = XGx
This proves that the matrices M D and T -1 -1
(G ) D O have the same eigenvalues. The matrix
(0^)-^D,G-^
is skew-symmetric and hence has imaginary eigenvalues. Besides that
1
_T^ _T .' U X D ^ X X^ = ^T
I T 7
X X X MxUsing theorem 2.3 we prove that the estimation (62) is valid.
2.1* The linear convection-diffusion equation.
2.1*. 1 A property of the pde (3)
We multiply equation (3) by u(x,t) and integrate over Ü = [0,1] to find, after partial integration and usage of the boundary conditions,
(61*)
(^'^)o * ^(^x'^)o* ^(^x'\)o = (S'^)o
In the same way as it was proved in section 2.2.1 we prove that the solution u of (3) satisfies:
(6=
|u(t)| < e-2^^|u(0)|L + ƒ e2^(^-^)||g(s)|L ds, V t e [ 0 , T ] .2.1*.2 The Galerkin formulation of (3)
For problem (3) it is proved in [7] that the space part of the solution u for all t £ [0,T] is an element of the space H iü).
° 1 1 Hence, if we use the finite dimensional subspace V (fi) of H iil), then the standard Galerkin formulation for (3) is
(66)
Find u such that
, h h , . , h h , , ^ , h h , , h>
(^'^ )o ^ ^ ( \ ' ^ )o * ^(\'^)o = (^'^ )o
with
u^(t) . v^(n), V t £ [o,T]
^ and u (O) given.
V v*" e v\n). o
(67)
Initial conditions follow from N
u^(x.,0) = u.(x.) = y u^(0) *.(x.), V j £ I
2.1*.3 The system of ODEs.
Problem (66) can be written;
(68)
...h „ h . ^ h „ Mu + aDu + bSu = FThe elements of the N x N matrices M, D and S are:
m.. = U-A.)^,
d(Ji.
•^ji = (dr'*j)o'
d((i. d * .
^ji " (dF'dx )o
The functions {*.}. - are a basis of V (fi) and we see immediately that
1 ' i£l o ''
M is symmetric,
D is skew symmetric and S is symmetric.
2.1*.1* Properties of system (68) From (66):
i^^\ d r 1/• h h> , , , h h> , ^, h hv , hv
(69) ^ {^(u ,u )Q} + a(u^,u ) Q + b(u^,\)Q = (g,u ) Q
Since u (x)] = 0 '' ; _ ' ' we find
, h hv , h hv
a(u^,u ) Q = -a(u^,u )„ = 0
Application of the Poincaré and Schwarz inequalities shows that
(TO) ||u^t)||Q < e-2^^||uNo)||Q + ƒ e-2^(^-")||g(s)||Q ds.
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 = M~\aD + bS) Consider Q (aD + bS)x = x'^lx ^--T ^--T X Dx X Sx „ (71) a i ï ~ "" ^ 3 f - = X X Mx X Mx
-36-Theorem 2.7 The eigenvalues of the matrix Q = M (aD + bS) belong to the rectangle given by
0 < Re(X.) < a 1 and when and with < Im(X.) < 6, i £ I ,S X max - D
x"
x^
x^
min max min X max = max X. l£l vS /,S> X max = max (X.) i£l ^ ,M . . ,,M, X mm = mm ( X . ) i £ l ^ ObservationsIn sections 2.2 - 2.1* 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
MÜ = Au + F.
M A
are almost imaginary if the convective terms dominate and aliiost 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 heat-diffusion 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.
-38-_3 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 SODE.
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], [1*1*]. In all these books the solution methods for the system
: i )
(2)
'± = Kz.t)
. ^ ( ° ) = ^
are mainly generalizations of the solution methods for the equation
y = f(y,t)
y(0) = yg
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 non-linear. It is well-known that non-linearity causes difficulties if implicit methods are used, i.e. for every time step a non-linear 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.
(3)
is approximated by the solution u of the numerical method: n+1 „, n ^n ,^ x u = F(u ,t ,At^) with > u = and
^
At - stepsize n -^Remark: The notational convention that u denotes the numerical approximation of the exact solution j£(t) at time t is adopted.
Assume that the solution ^(t) of the SODE is "smooth enough", then we find using the differential equations, that
Z(t""^ = F.;i(t'^),t^AtJ + d"
;
d is called the local truncation error.
:5)
Definition 3.1: The numerical method (3) is said to be of order p if the local truncation error satisfies
d"= 0(Atf')
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
time-interval 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 u is no larger than K6 in all subsequent values u , m > n, m,n £ IN .
This definition is dependent on the differential equation treated, so the idea of a "test-equation" is utilized. ([12]) The test-equation is the differential equation
-1*0-y = X-1*0-y , X £ C , Re(X) < 0.
This leads to:
Definition 3.3: The region of stability of the numerical method (3)
is that set of values At (real, non-negative) and X for which a perturbation in a single value u will produce a c
does not increase from step to step.
in a single value u will produce a change in subsequent values which
Suppose that the SODE
i = f ( x , t )
has to be solved.
Such a system locally behaves like the linearized system
^ - ' ^
3f I 3f
Z=Z , t = t
n n(^-^")'
In such a case we choose the stepsize At so that for all the eigenvalues X. of the Jacobian matrix
3f
the requirement that X. At is in the stability region, is satisfied.
Three classes of methods for ODEs In this section we consider briefly
(i) The Adams-methods, (ii) The Runge-Kutta methods, (iii) 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
1 Estimate the error.
2 Change the stepsize accordingly.
With respect to these properties the three classes mentioned above show the following characteristics:
_i Adams-methods (P(redictor )-C(orrector )-t^-pe.)
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 Adams-formulas is the one described in [1*1*].
ii Runge-Kutta methods (R-K)
During the integration it is difficult to estimate the error. If the classical R-K methods are used one can change the order only at the cost of additional function evaluations. As the R-K methods are one-step methods it is easy to adjust the stepsize.
In [10], [11 ] , [1*1 ] and [53] methods are derived that a.re based on the "R-K-idea" but which have error-estimators at the cost of little extra work.
iii Extrapolation methods (Ex)
Compared to (i) and (ii), these methods were developed quite recently ([ 2 ] , [13]). They combine the advantages of the Adams-methods and the R-K methods with respect to the points of (8).
Not only the properties mentioned in (8) are important: stability and accuracy play a role and also the number of times that the function f ' (in the equation y = f(y,t)) has to be evaluated, especially if this evaluation is difficult.
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 method is given, (See also [26])
Adams-Bashforth 2 (A-B-2)
/ n + l n At ,-_, n .n. „, n-1 .n-K, it Ü •*" ~^ (3f(u ,t ) - f(u ,t ))
LO
ü)
ON U ) LO O )I) U ) to5?
?L H H. <+«<«
•1 n> H' O 3 O til B (1) c + o Xm
>
r->
X U) 1 1 K J-IMAGINARY AXIS
1IMAGINARY AXIS
3.3_._2 Adams-Bashforth 3 (A-B-3) n+1 n At 10) n-1 ,n-l n-2 , n - 2 , u + T ^ (23f(u",t") - I6f(u"-',t"-') + 5 f ( u " - % t " - ^ ) ) 12 P = 3 3.3.3 Adams-Bashforth-Moulton 2 (P-C) n+1 ,* n ^ ,, „/ n n> u ' = u + Atf(u ,t ) 11 \ <* n+l n , ü-c ,_, n n, ^ „/ n + 1 , * , n+1 . . ,11) \ u = u + -- (f (u ,t ) + f (u ' ,t ) ) n+1 _ n . At 2 "-"-p = 2 3.3-1* Adams-Bashforth-Moulton 3 (P-C) (12) \ n+1,* n , At ,_„, n ,n, , n-1 .n-l,, u = u + — (3f(u ,t ) - f(u ,t )) n+1 n At , ^ „, n+1,* .n+1, „ , n n% „/ n-1 n-1,.
u = ü + "f^ (5f(u ' ,t ) + 8f(u ,t ) - f(u ,t ))
P = 3 3.3.5 R-K 1 (Euler explicit) (13) S n+1 n _^, n .n. u = u + Atf(u ,t ) 1 3.3.6 R-K 2 (Modified Euler) n+1,* n ^ , , _ , n n , u = u + Atf(u ,t ) (11*) u = ü + "^ t£(li >* ) "*•n+1 n ^ At ,„, n .n, ^ „, L^R '^ n+1,* .n+1,, )j p = 2 3.3-7 R-K 3 (Heun method) n+5 ,* n At _, n ^n, u = u + — f(u ,t )
(15) 1 u"*^'* = u"^ + Atf(u"*^'*,t'^-'5)
n+1,** n .^^, n ^ n+5,*,_ n+1,* .n+1, u = u + Atf(u -2u +2u ' ,t )
^—
\
1-3
1-3
-1*1*->
^
^
1 1 1-2 -1
3
2 X < o — IMAGINAR YREAL AXIS
fig- 3.3 Stability region rf methods 3.3-3, 3.3.7.
3.3.11-y^
.
1 13
2 X < > -• IMAGINAR Y-2 -1 "
REAL AXIS
fig. 3.1* Stability region of methods 3.3.8,
3.3.12-^n+l 1 (2u"*^'*+l*u"+^\u"^1.**_u")
P = 3
3.3.8 R-K 1* (Classical Runge-Kutta method)
(16)
u " * ^ ' * = u " + ^ f ( u " , t " )
u"^^'**=u"+^f(u"^^.*,t"-^)
u"^^'* = u'^ + At f(u"*^'**,t"^^)
u""^'** = u" + Atf(u""l'*,t""^
u"*1 = ^(2u""^ '*+l*u'^+^ '**+2u"^^ '*+u'^^1 '--3u")
( p = 1* 3-3-9 Trapezoid (implicit) (17) ^ n+1 P = 2 n At (f(u'\t") . f(u"^',,n-')) n+1 n+1, 3-3-10 Ex 2 (Extrapolation) -n+1 ( 1 8 ) u" + A t f ( u " , t " )
f*i = u " + ^ f ( u " , t " )
A ü"^1 = fl"+5+ A | f(2n+.^^^n+^ ,n+1 ^m+1 -n+1 u = 2u - u p = 2 3-3-11 Ex 3 ( E x t r a p o l a t i o n ) ü""^ = u " + A t f ( u ^ t " ) Ü " " ^ = u " + % f ( u " , t " ) gn+l ^ ^ n + ^ , A | f ( ö n + ^ , ^ n + i ^-l»6.
fig. 3.5 S t a b i l i t y region cf methcd
3.3.1
3.3.U
fig- 3.6 Stability region of method 3.3.2.
( 1 9 ) \ ~ 3 n At , n n , 2 1 1 1 n+-r n+-r , . n + r n+— u 3 = ^ ^ ^^LQ \t 3)
u"- = f ^ ^ ^ £(f V ^ )
n+1 1 , - n + 1 o->n+l ~n+1 , u = s ( u - 8 u +9u ) P = 3 3 . 3 . 12 Ex 1* ( E x t r a p o l a t i o n ) ( 2 0 ) ( -n+l n At „, n , n , u = u + At f ( u , t "') ^ n + 5 ^ ^ n + ^ f ( ^ n , t " ) 2"+^ ^ ^n ^ A t f ( a n + i ^ ^ n + i ^ -. 1* ..n+B , At „ , * n + 5 ^ n + j , u = ü + - 5 - f ( u , t ^) ..n+1 m + ^ _^ At „ , - ' ' ' T ^ / T " ^ u = u + — f (u . t ) n+1 1 ,, ^n+1 - n + 1 , u = T (i+u - u ) P = 1*.In figures 3-1 to 3.6 the stability regions of these methods are plctted.
Remark: The stability region of method 3.3.9 is the negative half-plane, i.e. the region with Re(X) < 0.
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.
Name AB2 AB3 ABM2 ABM3 RK1 RK2 Formula
(9)
(10) (11) (12) (13) (1U) Re-axis -1 .0 - -55 -2.0 -2.1* -2.0 -2.0 T Im-axis i .351 .TOi -1.2i -; Evaluations 1 1 2 2 1 2 RK3 RKl* Trap Ex2 Ex3 Exl*I (
(15) (16) (17) (18} 19) (20) -2.52 -2.8 — 00 -2.0 -2.5 -2.8 1.7i 2.8i ooi 1.7i 2.8iI 2.8i
3 1* 1 2 1* 5Table 3.1 Stability boundaries on real and imaginary axes and number of function evaluations per step At.
3.1* 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
y = f(x,t)
Here we shall consider briefly some methods constructed particularly for solving SODEs.
(21)
Consider the system of differential equations
y = Ai^ + f(t)
X(0) = Zo t £ [tQ,T]
The matrix A can be written as
(22) A = L + D + U
with L a lower diagonal-, U the upper diagonal- and D the principal diagonalmatrix. The Gauss-Seidel method is defined by the formula
-U9-(23) u"""^ = u" + At(L+D)u"''^ + AtUu" + Atf".
Remark: (23) may be followed by the formula
(2^) u'^*'' = u'^^^ + At(U+D)u"^2 ^ ^^^^n+1 ^ ^^,.n+1
We investigate the stability of method (23) for two special cases:
T T i) Suppose A = A , hence L = U
Then since J^(t) is not of interest in the stability analysis, we consider
(25) (E-At (L+D ))!!"•"^ = (E+AtL^)u"
with E : the identity matrix. n+ 1 T
We multiply (25) by (u ) , tranpose and add
(26) • (u"'^^^(2E-At(D+L)-At(D+L^))u"''^ = '
, n+l,T,„ . T, n , n,T,„ . ,, n+1 (u ) (E+AtL )u +(u ) (E+AtL)u
Multiply (25) by (u ) , transpose and add:
(27) (u")^(2E+AtL'^+AtL)u" =
• (u'')^(E-At(D+L))u"*^ + (u""^^^(E-At(D+L^))u"
From (26) and (27) we deduce:
(28) (u"*^-u")^2E-AtD)(u""^-u") + (/)^(AAt)/ = (u"^^^( AtA) (u'^*^ )
In the case that A is negative definite, since then .2E-AtD is positive definite, it follows that
(29) (u"*^)^(-AtA)u"^^ i (/')^(-AtA)u"
-50-/
i|x|lA = / ( - A ) x
it follows from (29) that I I n+1 I I I I nI I I i_ 11^ I i_ I 1^
and unconditional stability has been proven for both methods (23) and (23) followed by (2l*), provided that the matrix A is symmetric and negative definite.
T T
ii) Suppose A = -A , hence U = -L and D = 0.
Proceeding along the same lines as in case (i) results in;
(30) (u""'^^(2E+At(L+L^))u""'^ =
(u'')^(2E+At(L+L^))u" If the matrix
2E + At(L+L'^) = 3
is positive definite we find from (30) I I n+1 I I I I nI 1
l|u
| I B = | | U W^.Hence the method is stable and conservative. The property of conservation is also present for the analytical solution, for if,
T X = Ax with A = -A then: d T T X T (x x) = 2x Ax = 0 dt — — — — hence X X. ~ I \—\ I p is a constant in time.
The order of method (23) is easily found by means of a Taylor-series:
In the same way as the Gauss-Seidel method is related to the Gauss-Seidel iterative procedure, it is possible to relate a Jacobi method to the Jacobi-iterative method:
(32) u"""^ = u " + At(A-D)u" + AtDu"""^ + Atf"
For stability we investigate the homogeneous part:
(E-AtD)u"'^ = (E+AtA-AtD)u"
T . . . .
Suppose that A = A and that A is negative definite: Transpose, multiply and add to obtain
, n+1,T, At ., n+1 , n>T, „ At ., n (u ) (E-AtD+ — A)u < (u ) (E-AtD+ — A)u .
Hence method (32) is stable if the matrix
(33) Q = E - AtD + ^ A
is positive definite.
The order of method (32) is one: p = 1.
When instead of the linear system (21) a non-linear system is considered then the implicitness of the methods (23), (23) and (2l*), and (32) causes some problems. In that case a non-linear solver, like the Newton-Raphson 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.
However 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 oc and 6 such that the method
f n+1 ,* n ^ _ „, n ,n, u ' = u + Atf(u ,t )
(3M
h+1 n . / „/ n ^n, „„, n+1,* ^n+1 , u = u + At(af(u ,t ) + ef(u ' ,t ) Vhas 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 SODE resulting from a partial differential equation of diffusion type.
The solution of (3l*) is (see [33])
'>= 8 ' ^ =ÏÏ
In figure 3.7 stability regions of method (3l+) for some values of a
and B are plotted.
The second method of this section is designed for the solution of SODEs.
ODH *< 0 3 e C+- o-O p . U J
,
* r-—'
JO m>
r->
X(/)
1 4/1 r-IMAGINARY AXIS o -l^ji v; c r . |; 3 5 )
\
-51*-( u " ^ ^ ' * = u " + A t -51*-( A - D ) u " + AtDu"*^ ' * + At f" w h i l e t h e f o r m u l a n+1 n At / . n . n + 1 , * „n „n+1. u = u + - r ( A u + A u + f + f ) is used as a corrector, For (35) P = 2.In chapter l* 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 (ll*).
Methods (3I*) and (35) are stable if the eigenvalues X 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 R-Kl* method, (I6), has good stability properties along the imaginary axis (see also table 3.1), however it has order
p = 1*. Hence we investigate if it is possible to construct a R-K type method of order p < 1* that has an extended stability region along
the imaginary axis and that is more economic in use than the R-Kl* method. Before we start this investigation we give two definitions:
Definition 3.^: The amplification factor of a numerical method n+1 „, n ,n _ ,
u = F(u ,t ,At)
is the quantity P(XAt) in the expression
(36) u"""^ = P(XAt)u"
resulting from application of this niimerical method to the test equation
y = Xy
Definition 3.5= The effective time step At of a method, when X € C is given, is
-55-" e f f = f^P, N '
XteS' ' ^ ^
S is the stability region of this method and N the number of function evaluations per integration step.
Remark: The effective time step obviously depends on X.
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 R-K method that possesses this amplicication factor.
First we try three functions evaluations with second order accuracy, hence an amplification factor
2 2 3 3 :37) P(XAt) = 1 + XAt + ^ - | ^ + ^-^^.
id 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 3.1. The amplification factor (37) is maximally stable along the imaginary axis for a = 1*.
Proof: Suppose X = iy -»•
|P(iyAt)|2= (i.i4^)2, (^,,_ü!^)2
a
2.^2
1 + (uAt)'{i-f*^}
(i) When \i -*• 0, stability requires
r - - < °.
1* a ' hence: 0 < a < 8
-56-I* a 2 a 2»..2 ' ' a^
y At < 2a - r—.
So the value a that maximizes the expression
2 f(a) = 2a -
|-in the |-interval (0,8) satisfies the criteria. But
f'(a) = 0 if a = 1*, which completes the proof.
If we can find a method with the amplification factor (37) with a = 1* then the stability region along the imaginary axis will be
-2 < |X| At < 2
Comparing this to the stability region of (l6) (see table 3-1) shows that the effective time step of (l6) along the imaginary axis is longer than that for such a method.
Hence such a method would not be satisfactory and we therefore investigate also the amplification factor: (second order accuracy)
(38) P(XAt) = 1 + XAt + A ! | £ + 1 ^ + ^ A V .
Note that method (l6) has an amplification factor of this type (a=b=1). In fact we will show that (l6) has the optimal amplification factor of
this type.
Theorem 3.2: The amplification factor (38) is maximally stable along the imaginary axis if a=b=l.
Proof: Suppose X = iy ->•
2 2 1+1* 3 3
\-nt • .j.>|2 /- y At . bu At ,2 . , .. ay At ,2 |P(lyAt)| = (1 - ^^— + - ^ — ) + (yAt - - ^ — ) =
-57-1 + P At {(j^ + - ^ - -) + y At (- ^ + ^ ) + —yj^}
Consider the polynomial:
2 2
Q(x) = ^ X + (^ - ^ ) X - 3 + ^ + IT
A necessary condition for stability is Q(x) i 0, and in particular this should hold in the neighboiirhood of x = 0, hence:
(i) -f>l|4^°
and when equality holds we also require
From the results of table 3.1 we know that, if a=b=1
Q(x) 5 0, x2 S 8
Q(x) > 0, x2 > 8 •
We next intend to show that, for all a and b satisfying (i), we have Q{/8) 2 0 with equality only if a=b=l. This will prove the theorem.
On substitution we find, after some rearrangements
Q ( / ^ ) = I ( b - f ) 2 + 2 ( | . ^ ) 2 _ ^
using (i) leads to
Q ( / 8 ) . l ( b - f ) 2 + 2 ( ^ ) 2 - ^ = l ( b - l ) 2
and this shows that Q(/8) = 0 if and only if a=b=l.
Theorems 3.1 and 3.2 show that (l6) (the classical fourth order Runge-Kutta 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 • 1».
