• Nie Znaleziono Wyników

Numerical solution of the shallow water equations by a finite element method

N/A
N/A
Protected

Academic year: 2021

Share "Numerical solution of the shallow water equations by a finite element method"

Copied!
159
0
0

Pełen tekst

(1)

Lab.

v. Scheepsbauwkuntic

Technische

Hogeschool

Delft

NUMERICAL SOLUTION

OF THE

SHALLOW WATER EQUATIONS

BY A

FINITE ELEMENT METHOD

NIEK PRAAGMAN

2 .1081

(2)

NUMERICAL SOLUTION

OF THE

SHALLOW WATER EQUATIONS

BY A

FINITE ELEMENT METHOD

(3)

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

(4)

Dit proefschrift is goedgekeurd

door de promotor

prof.dr. EL van Spiegel.

(5)

met dank aan:

Phil en Sue Anderson,

Bert Broere,

Gonny Scheurkogel - v. Werkhoven

en voor al

Guus Segal.

(6)

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 14

2.2

the linear heat equation 19

2.3

a linear transport equation 28

2.4

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 40

3.3

some examples 41

3.4

systems of ordinary differential equations 48

3.5

extended stability regions 51

3.6 observation 59

4 Stability and numerical results 60

4.0 introduction 6o

4.1

methods to investigate stability 60

4.2

lumping 68

4.3

results for system

2.(26)

69

4.4

results for system

2.(53)

72

4.5

results for system

2.(68)

82

2

Onedimensional incompressible flow 84

5.0

introduction 84

5.1 boundary conditions 84

1,2

the Galerkin formulation 87

(7)

VIII

the system of ODEs 88

5.4 numerical stability 89

LI

computations and results 90

la

remarks (other methods) 96

11

conclusions 96

6 A numerical treatment of two dimensional shallow water flow 97

6.0 introduction 97

6.1

boundary conditions 97

6.2 the Galerkin formulation 101

6.3

the system of ODEs 102

6.4 stability 103

computations and results 105

7' 'Conclusions and remarks' 131

Appendix 1 133

Appendix 2 138

References 143'

(8)

dit prodschrift werd gedrukt bij. drukkerij & uitgeverij

EGNER by Den Helder

(9)

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 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. (BODE.). Chapter 3; Survey and discussion of numerical methods for solving SODEs.

(10)

-1-

-2-Chapter 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

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 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]).

(11)

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 Eulerian

coordinate 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 the

a

equation for a fluid reads:

xaxis

fig. 1.1

+ 1- + 2- (pv), +

.(pw) =

0

at 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) = pF

at

ax,

ay az a , , tpv).

+' kpuv), +

ax

L (

ay Pv + az PVW )1 = pF

Lt

a, a. a.

ow/ + (puw) + (avw) + ,(pw2) pF .

atax

ay

as

undisturbed water surface. The z-axis is chosen perpendicular to

yaxis

and y-axis. (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-,

(12)

-4-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:

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 22-x ;:c F : =

IR

Y P aY Fz :

p az

with p : pressure

(6)

z

(13)

S97

PO,

412),

Friction (molecular viscosity)

Fx

vu

(8)

F vAv

= yaw

with ,2 ,2 A =

2-- + a

+ a2

ax

ay

ay2 322

v d, = 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)1, (3) and (4) leads to

a,

a , 2, 3 , , a , ,

43.7 IbuA + i; kpu )

+ 7 OUVJ 4- 7

kpuvi = p(201317-20w) - + pvnu

ay

az

32 and

a,

, a I , a , 2, a .

ov) +

tpuvm + ) + ,pvw) = pc2/1w-2111 3t ax ay az t+T1 1 r, 71,-

I

uo dt1=101

Here 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.02u-201 - pg - +

at ax

ay

az

aZ

Since we are only interested in the tidal motion, we apply, to equations, (1), (9), (10) and (11), a moving time-average 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 (1-i,7,-,;) and a fluctuating part (u' ,v' ) in such a way that

t+T 1 r

T

u (It = U. Fz : [18], -5-=

(14)

(17)

-6-Before we apply this time-averaging 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 well-nixedw

Results of the time averaging procedure are:

Equation '(1; aa aT;

ax

ay Equation (9) _

3uI

4. -2 a, -- 3

-

1 _a-P.

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 + 2-

f

(u,w, )1 dt/

T ax 3y 3z t t t

Equation 06)

_

-

1 li 1

(16')

av a

017) ,

(,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 3

f

,lt

r,

_ - {

tU V h uT +

7/. j Orkll dT + 2-an

f

(V'Wt ) diE 1 T ax ay t t 1 t Equation (11)1

Ai

a

(--)

a

(--,

a

-2

at

ax

''w, ji orw, (w ) = 2a2 _ 20 1 16 - P az Vt14

t+T

a

t+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

non-linearity 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

-

(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 - -2-p 33. ) (v+vT, Au ( 19 ) 77 + kuv) + -2

()

.57

(v.)

-- = 2Q1W 22 u -3

p ay

+ (v+vT) + +

00)

(22) K217

(16)

-8

423)1 aw af (--\i 4. 9

(--)

A_ (72, LzU2 ? -2 at

az

ay

'1

az

.

A in s 2Q2U - 2Q1 MO-4

m

s-2 = g - VO m 5-2 1

Ai=

-I LzU 1 -2 vTA1-4 m 10

ms

Lz

vA; m 10-6 LzU

1

Ins

2

L 2

Lz

Theorder of the values that should be substituted for a1' Q2 and v is well-known. 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 so-called hydro-static 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 x-y 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

(17)

.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. " w

The integration formula (n stands for x, y or t)

&

,

21ah

f

'

. dz = 1-

f

f dz - f(C) - f(-h) --an an an an -h -h is also used. (314) (35) I 323' ay

(18)

-9-A,

-

-10-Integration, of

ilk)

leads to

t

-

E J au r ay

dz +dz +

;(x,,y,E

,t ) - iTg(x,,y ,-h ,t = Dx ay -hi -h

Combination 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 Dx

If we neglect the slow variation in time of' the bottom, i.e. if we take

at

we find

ii

3 -, a

Kliu, + Kay)

at ax

Integration of the non-linear 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 ah

--Kliu)

u(x,y,C,t) at at r

-The 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)

(19)

(43)

- a -2 Tc. (u ) dz =

2- (

f

712 dz) 7.12(x,y,E,t) ax ax h -h -2 ah a

f

2 u (x,y,-h,t) + S ) dz = ax ax -h ar 2 -2

al

-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) 5

T- (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 right-hand side terms results in:

5

f.i.dz=fHv

h E (47) 5 g a dz = H g

ax

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 +

(20)

(50)

(51)

-12-Combination of (46), (47), (48),

(37)

and

(38)

gives:

(49)

(H71) + (HZ2) + -a--- (H7;) = at ax ay n2-

n2-=PH;

(v+vT)

f

(=-1 +

=a) dz

ax

-h ax ay3y2 a ( f& a2 dz) - a

f

dz)

(v+vT)11)

ax ay -h 3u - (v+v T az z=E z=-h

-The 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 two-dimensional wind velocity vector near to the air-water 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

+

(21)

-13-C 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 diffusion-like 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 so-called "shallow-water 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

(22)

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

(23)

-15-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:

(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 well-known 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

(24)

f(l6q

Xu,v)1 =

fi

(uv + vu.Vv)

dx = (u,v)0, + (01.1,Vv)0

Two norms, related with

(5)

and

(6)

are respectively

Nul I g =

(u,110'0

and

= (u,U)Iir

One can prove that boundary

sense for the spaces H1 ,

H1o

unstable ([24], (45]). Such

-16-conditions 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 of

0(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 these

spaces 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 ;12-norm:

(25)

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 orthogonal

eigen-vectors.

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 = (C-CT)

When the eigenvalues of Si and T are respectively

and

then

(i)1

min X! 5

fie(JS)

max A

icI 1

1E'

. 1

(ii) -

max A. 5 Im(),) 5 max lAil

iEI

iEI

= {1,2,...,N}

1.

V 5; c I

(26)

-17-(10)

-T x x

/ lail2

i=1

A real symmetrix matrix has real-valued 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

-18-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.

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 -TSx

i1

A7. lail2 x =

Remark: When A is an eigenvalue of a skew-symmetric matrix, then X = - A

is an eigenvalue too. Hence (11) estimates all the eigenvalues of T.

Consider Cx = Ax

x x

Tx x x

(27)

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 non-homogeneous 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] t

E [0,T]

then: -19-T x Cx xT(i(C+CT) + (C-CT))x

So non-homogeneous 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

(28)

-20-a = [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 Poincare-inequality (u,u)0 5 i(ux,ux)0

and the Schwarz-inequality

(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)20

-75 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 e-2tliu(0)110 ( e2(s-t) s)1

lie,11

k 10 ds V t

E [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 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).

2

(29)

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 approximate

u by the solution uh of the problem: (standard Galerkin approach,. see [L5],)n

. h

Find u such that

(uh vh) (uh vh)

in

vhk

x' 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=1

N, do. do.

V un r 1 --11 e, ) L i 'dx '0"ei'wj 0 1=1

Initial 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 Galerkin-weighted 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)

(30)

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 + 1

uh

(--1'' --'L) )0 i dx = (64j )0 ' V j E. I. 1 dt ((Pi4j -i=1 i=1

In 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

(31)

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,

the

0 (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 by

2-1

-1

2--1 0

\\\

-1 2 -1 0

1 2

=23-For the triangulation of( = {0,1]we use an equidistant stepsijize Ax and define jar x. = 1

e. = Ex.x.].

1 1 H1o(0).

(32)

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 = j-1 i=1 j=1 2 = (

x.(1).(x),

x.0.(x))

= II 1

x.0.(x)110

J J 0 j=1 1 1 i=1 J-1

Since the functions

4i(x)

form a basis of V1(0) they are linearly0

independent, hence

x.c1).(x) 0 0 when x

0

0

and we conclude that

xTMx 0

V x

o

The matrix S yields a similar property

-24-0

01

0 11N+1 11N+1

(36)

F =

M gi

g2

Ax 6 go o o Ax - 6 . gN 0 gN+1 . --i'.1o 1 0 I 0 . 1 Ax o o i=1 0

(33)

T, x ox

=I

x.( 2 s. .x.) = i=l J-1 N dO.

= Vx-)

i=1 j=1 dx dx 0 j dØ. N de. (

/

Xi

at,

2

x. .;.1)0 i=, dep. II Jc.

1112

dx 0

i=r

By Poincare: -

07)

TSx = ( ximi),a7(

xickino

i=1 i=l a. 2( x.01,, 1 '=1 when x Oi

Hence! xT's > 0

Vx

0 and the theorem is proved.

:Assume that the (real)eigenvalues of M satisfy::

M I

)5 NM 5

5 AN-11 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 -N-11 5 -11

In section 2.2.1 a stable property of the pde was found.

h 1 Since u c

Vo0),

(22) gives -25-i=1 =

-x.4).)

1=1 =

(34)

d1(uh uh) 4 + (uh uh) = uh) 10,0) dt 2 ' 0 x' x 0 ' 0, or in matrix-vector 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 M-1S

are important. These eigenvalues can be estimated if we use:

Theorem 2.5 The eigenvalues

xS, i

I of

-1 Q.= MI S

are real, positive. and

:(142 f( 44,)! or

xs

xs N

x9

1 AN AN Sx = AQMx ,Sx =

XYGx

ITST ) -1 1 SG Gx

This expression shows that the matrix

T -T

1(G ) SG

-26-V i c

I.

Proof: M is symmetric and positive definite hence there'exists a

Choleski decomposition (

43)

M =TG.

Consider -1 Q

M

Sx = A x (41) A

(35)

(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

M-1Sx =Qx

TSx =

AQT

x x

_

Mx 4-TSx xT x x Tx xT x MX

AQ

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.

(36)

-

-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)0

Equation (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 well-posed ([35]).

2.3.2 The Galerkin formulation of (2)

For problem (2) we assume that the space part of the solution u is

(37)

-29-A 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 the

unknowns

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=1

Initial 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:

(38)

and m.. = (o.4.) ji j 0 ,do. d.. =

dT4j)0

1

For 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

0

12

-30-Remark: 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 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 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

(39)

-31-Since the {cki(x)} , i c {102,...,N+1} form a basis of the space

V1 ([0 1])

"

not all Oi(x ) satisfy E

the 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:ET-L1

_

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)) 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

(57):

(55

+

(40)

+.32= d

hTh1h2

h T 1(58) ((1.1 ) 1N+1I = i(u 1 Y 1, .2 1

The 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 = M-1D

are important.

These eigenvalues satisfy:

, (60) =

1-i

or 1 Dx. = X.Ma. 1 T Define DS = -(D+D2 N and 1 T. DT = (D-D them T -T x.D x. x.D x.

-1

s-a -a T-a + -T -T x.Mx. x.Mx. 1 -a

-1 -1

The term T x.D x. S-3. T 1 -a

that 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 M-1AT, are estimated in the following way:

1(6N) a

1-i

2

-h,T M-1Dx. I

(41)

-.33--1

Theorem

2.6

: The eigenvalues Au of the matrix U = M

DT are imaginary and satisfy DT max IA I 1

(62)

max IX2,1 I

min IXM.I9.

Proof: Consider M-1D

x =

Aux

M is symmetric and positive definite hence there exists a Choleski decom-position M = GTG.

Then:

D x = AGTGx

Or

GT)-1DTG-1Gx = AGx

This proves that the matrices M-1DT and

T-1

, -1 (G

DTG have the same eigenvalues. The matrix

(GT)-1D G-1

is skew-symmetric and hence has imaginary eigenvalues. Besides that -T -T =

xDxxx

U

--A -T -T x x x Mx

Using theorem

2.3

we prove that the estimation

(62)

is valid.

2.4

The linear convection-diffusion 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

(42)

'466)

-34-K 64 4'

Kuu.)

+ n(u + b(u ,u ) =-,Kg,u)r

t

0

x

0

x

x 0

In 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, e-2btru.a. I I

)110

1 e2b(5-t)

II)11

g(s1.0 ds, ?ke t E

2.4.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 c [0,T] is an element of the space

Hence, if we use the finite dimensional subspace

v1(o)

of 11011(0,, then

the standard Galerkin formulation for (3) is

Find uh such that

hh

,

' Of

h h h h h

(u,v)

+ aku v ) + bku (s,v)0

vhr 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 (

(43)

(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 2

xx0

' 0 Since uh(x)]10 = 0 we find a(uh,uh) = -a(uh,uh) 0 x 0 x 0

Application of the Poincare and Schwarz inequalities shows that

Iluh(t)110 e-2btiluh(0)/1 a T x Mx + b =

-35-Now we can use theorem 2.3 to prove:

0 + f e-2b(

0

t-s)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 = M-1(aD + bS) Consider = x (aD + bS)x = AQMx T -T x Dx x Sx

_

T x Mx S d

(44)

and

with

D

max

A max = max AD.I

id

I 1

A max = max tA.)

id

1

M

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 well-known that the eigenvalues of the matrix

-36--1

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.5

(45)

M-1A

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 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.

(46)

-37-

-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 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 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.

Suppose that the solution y(t) of the system

{i

= f(y,t)

(47)

(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 p

if

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

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 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 "test-equation"' is utilized. K[12])' The test-equation is the differential equation

(48)

'0.6) = ASp.

AEC., Re00

S

This leads to:

Definition 3.3: The region of stability of the numerical method (3Y

is that set of values At (real, non-negative) 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 )(Z-Xn)

+ 1(in,tnY

+ ; n n(t-tri)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 Adams-methods,

(HY

The Runge-Kutta 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.

(49)

(9)

With respect to these properties the three classes mentioned above show the following characteristics:

Adams-methods (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 Adams-formulas is the one described in [44].

ii Runge-Kutta methods (S-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], [41] and 1:53] methods are derived that are based on the "R-K-idea" but which have error-estimators 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 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 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 method

is given. (See also [26])

3.3.1 Adams-Bashforth 2 (A-B-2) n+1 n At = 4. (3±'(t°)tn) f(un-1,tn-1)) ' s u u

_

_

p=

2

(50)

-.41,-

-42--3

-2

-1

REAL AXIS

fig. 3.1

Stability region of method

3.3.5.

3

-2

0

-3

-2

REAL AXIS

fig. 3.2

Stability region of methods

(51)

3.3.2 Adams-Bashforth 3 (A-B-3)

n At n n-2 n-2

un+1 = u + (23f(u,tn) -

16f(u,tn-1)

n-1 + 5f(11 ,t ))

(io)

P= 3

3.3.3 Adams-Bashforth-Moulton 2 (P-C) un+1,* = un + Atf(un tn), un+1 un At

,J1,

f(un+1,*,tn+1)) 2 (1.-"1 '6 )

p= 2

-43-n+1,** = un + Atf(un-2un+'*+2u114-1'*,t11+1)

3.34

Adams-3ashforth-Moulton 3 (P-C) f(un-1,tn-1)) (12) n+1,* n At n n

= u

+ 7

(3f(u,t) - f(un-1,tn-1 U )) n+1 un Alt2 (5f(un+1,*,tn+1 11 ) 8f(un,tn) = 3 3.3.5 R-K 1 (Euler explicit) n+1 n n n

u=i + Atf(ut)

(13) 1 3.3.6 R-K 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

R-K 3 (Heun method) un+,* un At f(un tn) 2 (15) un+1,* un iltf(un+i,*,tn44)

p

At

(52)

-3

-2

-1

REAL 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.3-3, 3.3.7, 3.3.11.

-3

0

of methods

(53)

(16)

(17)

n+1 1

T2un+,*+u

-u.)

P= 3

3.3.8 R-K 4 (Classical Runge-Kutta method)

un+3,* n At n n

_

= -1. + 1(11- ,t ) unig ,** = un + AL f(un+i'*,tn-ig)

-1 un+1,* un + At f(un+i'**,tn+)

_

_

un+1 '** = un + Atf(un+1'*,tn+1) -** n+1

un+1 -

-(2u

1.

' +4u

* n+'

2'

+2u

'

*+un+l

'**

-3un)

-

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

fie-1g

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 - u

p= 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

-

-45-n+i -u

(54)

-46--2

-1

REAL AXIS

1

0

fig. 3.5

Stability region

of method

3.3.1

3.3.4

2

1

0

fig. 3.6 Stability region

of method

3.3.2.

-2

-1

(55)

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 ftu

-n+1n

At , n n = u + 411)-U

P+1

n At 2

1

n+-4 -n+4 At ,.,n+;1

;I

3 ri. = u + -- ftu - tr1+41 3 i n+11

-nq

At

.4

n 4 ti = il + --i f(11 ,t )I n+1 1 (4-n+1 -n+n, u 3

p=

In figures

3.4

to

3.6

the stability regions of these !methods are plAtted4

RemArk: The stability region of method 3.3.9 is the. negative half-plane, 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,t

3)

1 3 + =

(56)

-48-Table 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 Gauss-Seidel method is defined by the formula

Name Formula Re-axis Im-axis 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 -., .

(57)

(un+1)T(

-49-n+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

(E-At(L+D flun+1 = (E+AtLT)un

with E : the identity matrix.

We multiply (25) by

(u1)T,

tranpose and add

r T , n+1

(un+1)2E-At(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 ) kz-At(D+L flu

From (26) and (27) we deduce:

n+1 n,T

(26) ku

-u )

,k2E-AtD)(un+1

-u)

n (u ,Tn) (AAt)un = (un+l)T(AtA)(un+1) In the case that A is negative definite, since then

2E-AtD 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

+

(58)

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 =

(59)

-51-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:

un+1 = un + At(A-D)um + AtDun+1 + Atfn

For stability we investigate the homogeneous part:

(E-AtD)un+1 = (E+AtA-AtD)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 ) (E-AtD+ A)u s (u ) (E-AtD+ 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 non-linear system is considered then the implicitness of the methods (23), (23) and (24), 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.

(60)

(34)

{

unfl'*

= u+

n Atf(un tn)

-52-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 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.

(61)

-5i-I

-8

-

-St

-4

-;

-/

-t

REAL AXIS

Stability

region method 3.(34)

I

1

for a = 8 , 5 _

--3

-8

-7

-6

-5

-4

-3

-2

REAL AXIS

fig. 3.7.1

Stability

region method 3.(34) 3

for a

= T ,

5 0 T.

-3

-2

-8

-;

-6

-5

-4

-3 - 2

-

; REAL AXIS

fig. 3..7.2

Stability

region method 3.(34)

(62)

-54-un+1'*

=

un

+

At(A-D)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 R-K4 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 R-K 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 R-K4 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

(63)

-It

Ateff = 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 R-K 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-

2-111-12.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

(64)

-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

criteria-But

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) shows

that the effective time step of

(16)

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)

2,42

4 4

At3 to, At

C30

PULAt)= 1 + AAt + =-S=- +

14-24

-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 of

this 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 + 24

(65)

lc I

-57-2) u4,st4b24 4 b A) +

1.12St2(-g

ie

576

. 1 * y OM4

1

f(r,

71

- 3

Consider the polynomial:

Q(x), b2

4(a2

-th 2 a b 1

576 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 (3

and when equality holds we also require

a2

0.n7-Tb-co

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)2

9 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 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 b b

(i)

(b-+ -=

(66)

(p)t

'LK

2 0

-58-I I

li

2i

3i

. 3.8 Dissipation of the R-104 method.

1

2i 3i

Cytaty

Powiązane dokumenty

Pierwsza charakteryzuje się utratą masy ciała oraz rozmieszczeniem mas tłuszczakowatych w obrębie górnej części ciała, a także wewnątrz jam ciała (najczęściej

Tego rodzaju transakcje mogły n a przykład zaistnieć wówczas, gdy strony były związane ze sobą bliskim i więzami rodzinnym i lub sąsiedzkim i, a po drugie m ogła

(im więcej dusz potępionych będzie m iał na sw ym koncie, tym większą zasługę zdobędzie u zw ierzch­ nictwa) raczej z obowiązku niż z entuzjazm u. Lecz

Verbeteringen worden voorgesteld, zoals het gebruik van gordels om de rug beter te ondersteunen en een grotere bewustwording van de kraanmeester ten aanzien van zijn eigen

Zostaây one ukazane na tle sĊdów innych specjalistów na temat specyfiki nauczania języka polskiego jako obcego (dalej: JPJO). Przyjęte przez Miodunkę zaâoůenia moůna

Segmental challenge with the individual top-dose of 100 μg carbon nanoparticles showed a significant relative increase of neutrophils (p = 0.05) in peripheral blood as compared

W celu zabezpieczenia się przed kon- sekwencjami uroków rzucanych przez czarownice zbierające się o pół- nocy na sabat, podczas wielu świąt trzaskano z bicza.. Miało to znacze-

The verification of the performance of the DARS algorithm for the RC structures re- vealed that the number of exact function evaluations (NLFEA based) depends on the target