Statistical mechanics and numerical modelling of geophysical fluid dynamics - Thesis

(1)

UvA-DARE is a service provided by the library of the University of Amsterdam (http

s

://dare.uva.nl)

Statistical mechanics and numerical modelling of geophysical fluid dynamics

Dubinkina, S.B.

Publication date

2010

Document Version

Final published version

Link to publication

Citation for published version (APA):

Dubinkina, S. B. (2010). Statistical mechanics and numerical modelling of geophysical fluid

dynamics.

General rights

It is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s)

and/or copyright holder(s), other than for strictly personal, individual use, unless the work is under an open

content license (like Creative Commons).

Disclaimer/Complaints regulations

If you believe that digital publication of certain material infringes any of your rights or (privacy) interests, please

let the Library know, stating your reasons. In case of a legitimate complaint, the Library will make the material

inaccessible and/or remove it from the website. Please Ask the Library: https://uba.uva.nl/en/contact, or a letter

to: Library of the University of Amsterdam, Secretariat, Singel 425, 1012 WP Amsterdam, The Netherlands. You

will be contacted as soon as possible.

(2)

Statistical Mechanics and

Numerical Modelling of

Geophysical Fluid Dynamics

Svetlana Dubinkina

Sta

tistical M

echanics and Numer

ical M

odelling of G

eoph

ysical F

luid D

ynamics S

vetlana D

ubink

ina

Uitnodiging

tot het bijwonen van

de openbare verdediging

van mijn proefschrift,

getiteld

Statistical Mechanics and

Numerical Modelling of

Geophysical Fluid Dynamics

op vrijdag

28 mei 2010

om 10.00 uur

in de Agnietenkapel,

Oudezijds Voorburgwal 231,

te Amsterdam

U bent ook

van harte welkom

op de receptie na afloop

Svetlana Dubinkina

CWI, Postbus 94079

1090 GB Amsterdam

(3)

Numeri al Modelling of

Geophysi al Fluid

Dynami s

(4)

PrintedandboundbyIpskampDrukkers.

CIP-DATALIBRARYUNIVERSITEITVAN AMSTERDAM

SvetlanaBorisovnaDubinkina

Statisti alMe hani sandNumeri alModellingofGeophysi alFluidDynami s/

bySvetlanaBorisovnaDubinkina . Amsterdam: UniversiteitvanAmsterdam,

2010. Proefs hrift.

(5)

ACADEMISCH PROEFSCHRIFT

terverkrijgingvandegraadvando tor

aandeUniversiteitvanAmsterdam

opgezagvandeRe torMagni us

prof.dr. D.C.vandenBoom

tenoverstaanvaneendoor

het ollegevoorpromotiesingestelde ommissie,

inhetopenbaar teverdedigen indeAgnietenkapel

opvrijdag28mei2010,te 10:00uur

door

Svetlana Borisovna Dubinkina

(6)

Promotor: prof. dr.ir.J.E.Frank

Co-promotor: prof. dr.J.G.Verwer

Overigeleden: prof. dr.R.P.Stevenson

prof. dr.C.A.J.Klaassen

prof. dr.P.W.Hemker

prof. dr.B.Leimkuhler

prof. dr.S.Rei h

prof. dr.W.Hundsdorfer

Fa ulteitderNatuurwetens happen,WiskundeenInformati a

T

HOMAS

S

TIELTJES

I

NSTITUTE

FOR

M

ATHEMATICS

Hetonderzoekdattotditproefs hriftheeftgeleidwerdmedemogelijkgemaakt

doorhetprogrammaKlimaatvariabiliteitvandeNederlandseOrganisatievoor

Wetens happelijk Onderzoek, gebied Aard- en Levenswetens happen, onder

(7)

(8)

(9)

Preface

Statisti alme hani sisapowerfulapproa htounderstanding omplexphysi al

systems. The purpose of statisti al me hani sis to onstru t methods whi h

an handlein ompletelyknown systems;todes ribethemostlikelybehaviour

ofasystem;toworkwithtimeseriesdatafromexperimentsornumeri al

mod-els. Inthisthesisweusestatisti alme hani sasatooltoshowtheimportan e

of onservation laws of a numeri al method. Depending on what onserved

quantitiesthedis retesystemhas,thestatisti altheoryofthisdis retesystem

varies. Therefore,the hosenmethodhasaninuen e onthestatisti alresults

of thesimulatedmodel. Inthe thesisIshowthat statisti alme hani s anbe

alsoemployedasatoolbyanumeri alanalysttoverifythestatisti ala ura y

of a numeri al method. This is very important issue for appli ations su h as

limate variability, sin e in these appli ations long numeri al simulations are

run for dynami al systems that are known to be haoti , and for whi h it is

onsequentlyimpossible tosimulate aparti ularsolutionwithanya ura yin

the usual sense. Instead, the goalof su h simulations is to obtain a data set

suitable for omputing statisti al averages or otherwise to sample the

proba-bilitydistributionasso iatedwiththe ontinuousproblem. Dierentnumeri al

dis retizations havedierent dis rete dynami s. Therefore it is ru ial to

es-tablish the inuen e that a parti ular hoi e of method hason the statisti al

resultsobtainedfrom simulations.

Anotherattra tiveappli ation forstatisti al me hani s is themodeling of

sub-s ale motion. Invis iduidmodelsarenatural in a number of appli ation

areas, su h as atmosphere and o ean s ien e. These ows are hara terized

by onservation of energy, sensitive dependen e on initial onditions, and the

as adeofvorti itytoeverners ales. Forthesimulationofsu hows

numer-i ally,thevorti ity as adepresentsthe hallengethatanydire tdis retization

of the equation of motionmust eventually be omeunderresolved. To address

thisproblemee tivelyrequiresmodellingthesub-grids aledynami sandtheir

inuen eonthe oarses ale. Inthethesis,Ishowforapointvortexmodelthat

these ee ts an be parameterized using an adapted 'mathemati al

thermo-stat', ate hniqueusedin mole ularsimulationstomodela systemofparti les

intera ting ata onstanttemperature. I believethat this methodology anbe

extended for more omplex models with feasible appli ations in limate

vari-ability.

This proje t was funded in NWO Earth& Life S ien es Coun ilClimate

Variability program,anditaddress theproblemsarisingin numeri al

simula-tionsof geophysi al uiddynami sstatisti al a ura yof methods, inuen e

(12)

a newte hniqueformodellingthesub-s alemotion,whi h anbeemployedin

atmosphereando eans ien e.

Thisthesisre ordsthenumeri almathemati sresear hI ondu tedbetween

August2005andJanuary2010intheModeling,AnalysisandSimulation(MAS)

departmentoftheCentrumWiskunde&Informati a(CWI)inAmsterdam.

Chapters2through 4ofthisthesishaveappearedasjournalarti les:

1. Chapter2isbasedonS.DubinkinaandJ.Frank,Statisti alme hani sof

Arakawa's dis retizations, Journal of Computational Physi s 227, pages

12861305,2007.

2. Chapter 3 is based on S. Dubinkina and J. Frank, Statisti al relevan e

of vorti ity onservation with the Hamiltonian parti le-mesh method,

a - eptedforpubli ationinJournal ofComputationalPhysi s,2010.

3. Chapter4 isbasedonS.Dubinkina,J.FrankandB.Leimkuhler,A

ther-mostat losureforpointvorti es,submitted,2009.

Chapter 1,theintrodu tory hapter olle ts mu hof theba kgroundmaterial

neededtoread therest ofthethesis.

SvetlanaBorisovnaDubinkina

(13)

Chapter 1 Introduction

In numeri al simulations of limate variability, the ow of atmosphere and

o eansisknowntobe haoti ,andthereforethequantityofparti ularinterest

isnotasinglesolutiontraje tory,butratheranensembleofpossiblemotionsof

thesystemgiven ertaininformationaboutthebehavior ofthissystem. Then

statisti alme hani sbe omesuseful,sin eitdealswiththeprobabilityofa

sys-tembeinginsomestateunderknown onstraints. Toderiveastatisti altheory

for a system, the system should be volume preserving and possess onserved

quantities. Ergodi ityisusuallytakenas anassumptionto relatedynami s to

statisti al theory, and thevalidityof this assumptionmust be onrmed with

numeri al experiments. Although a ontinuous model has onservation laws,

it isnot learthat after implementinga numeri al method, the orresponding

dis retemodelwillalsohaveanequivalen eofthem. Thisentirelydependson

propertiesof a numeri al method. Thereforethe statisti al resultsof the

dis- retemodelmaydepend(andtheya tuallydo)onthenumeri almethodused.

And in an ample range of available numeri al methods, it therefore be omes

even more ru ial and pronoun ed to understand how the numeri al method

inuen esthestatisti alme hani sof asimulation.

Arakawa'ss hemes for 2D in ompressible ow area onvenient subje t to

investigatetheimportan eof onservation be ausetheysatisfyoneor bothof

a pair of onservation laws, and they are onstru ted from identi al dis rete

operators, dieringonly in the order of appli ation of these. If the operators

would ommuteinthedis rete aseastheydointhe ontinuousthes hemes

would be identi al. Nonetheless we observe ompletely dierent equilibrium

behavior for these methods, see Chapter 2, and for the introdu tion to the

statisti al me hani sseeSe tion1.3.

Symple ti or Poisson dis retizations for Hamiltonian partial dierential

equations are onstru ted in su h a way that the semi-dis rete systems are

again Hamiltonian and possess an equivalen e of onserved quantities of the

original ontinuousHamiltoniansystems. Be auseHamiltonianstru tureplays

animportantroleinstatisti s,wemightexpe tthesemethodstoperformbetter

thanstandardmethodsintermsofstatisti alperforman emeasures. This

mo-tivatedus toderiveinthethesisstatisti altheoriesforaHamiltonian

parti le-meshmethodin the aseofquasigeostrophi owwithtopographi for ing,see

(14)

Se tions 1.1and1.2.

Anotherattra tiveemploymentofstatisti al me hani sis themodelling of

sub-grids ale losureof invis id uidequations, sin esu h ows,whi hmodel

essential omponentsof atmosphereand o ean, support a vorti ity as ade to

ever ner s ales. Any nitedis retization ofthese equationsmustbe ome

un-derresolved,and we an useequilibrium statisti al me hani stheory to

repre-sent the subgrid motions. In Chapter 4, we show for a point vortex ow on

a dis that this an be parameterized byusing a mathemati al thermostat,a

te hniquefrommole ulardynami swhi hallowstomodela systemin onta t

withareservoirofnes alevorti ity. InSe tion1.3,weexplainstatisti al

me- hani sin general, on ept ofmathemati al thermostat, and dene statisti al

temperature.

1.1 Hamiltonian systems and geometric

integra-tion

The ontinuum models studied in this thesis have a ommon mathemati al

stru ture: thatofHamiltoniansystems. Hamiltonianstru tureisalsothepoint

of departureforstatisti al me hani s. WethereforereviewHamiltonian

stru -tureforordinaryandpartialdierentialequationsinthisse tion,andwedis uss

geometri integratorsnumeri almethodsthatpreserveHamiltonianstru ture.

A geometri integrator isa numeri almethod whi h preservessome stru tural

properties of a given problem, e.g. symple ti ity, rst integrals, symmetries.

Hamiltonian systems in turn are ri h in su h properties. Thereforethis

moti-vatesinterestingeometri integrationofHamiltoniansystems.

1.1.1 Hamiltonian systems

Inthefollowingse tionwe onsider Hamilton'sprin iplefornite-dimensional

me hani alsystemsfollowingthederivationsdes ribedin[3 ,42 ℄. Inthe ontext

ofgeophysi aluiddynami ssee also[57 ,73,74,82℄.

The canonical equation

Consider a me hani al system with

d

degrees of freedom. Its state an be spe iedbythegeneralized oordinates

q

= (q

1 , . . . , q

d

)

T

asfun tionsoftime

t

. Thedynami softhesystemaredeterminedbytheLagrangian

L

,thedieren e between the kineti and potential energies. The Lagrangian

L = L(q, ˙q)

is a fun tion of the oordinates

q

and thevelo ities

˙q

, where the dots on

q

stand for dierentiation with respe t to time. The evolution of the system may be

determinedfrom Hamilton'sprin iple[42℄

δ

Z

t

0

(15)

forarbitraryvariations

δq

k

thatvanishat

t = 0

and

t = t

0

. Fromthisvariational problem(1.1)followsLagrange'sequationsofthesystem

d

dt

µ ∂L

∂ ˙q

¶

=

∂L

∂q

.

(1.2)

Denethe onjugatemomenta as

p

k

=

∂L

∂ ˙

q

k

(q, ˙q),

k = 1, . . . , d

(1.3)

andthe

Hamiltonian

H

_{≡ p}

T

_˙q

− L(q, ˙q)

asafun tionof

q

and

p

. ThisistheLegendretransformationof

L

withrespe t to

˙q

. Undertheassumptionthat

˙q

maybeexpressedasa fun tionof

q

and

p

via (1.3),we an hangevariablesto

(q, p)

.

Hamilton'sprin ipleintermsofvariables

(q, p)

leadsto

Hamilton’s canonical

equations

˙

q

k

=

∂H

∂p

k

(q, p),

p

˙

k

=

−

∂H

∂q

k

(q, p),

k = 1, . . . , d.

(1.4)

The anoni alequation(1.4) ontinuestoholdif

L(q, ˙q, t)

,andhen e

H(q, p, t)

, ontainsanexpli ittime-dependen e.

Itisusefulto ombineallthedependentvariablesina

2d

-dimensionalve tor

y

= (q, p)

. Then(1.4)takesa simpleform

˙y = J

∇

y

H,

(1.5)

where

J = (J

ij

)

isthe

2d

× 2d

skew-symmetri matrix

J =

Ã

0

1 −1 0

!

.

(1.6)

Here

1

and

0

representtheunitandzero

d

× d

matri es,respe tively.

The Poisson bracket

An alternativeform for a Hamiltonian systemis the

Poisson bracket

, a skew-symmetri , bilinear form. For two fun tions

F (y)

and

G(y)

this is dened as

{F, G} = ∇

y

F (y)

T

J

∇

y

G(y).

(1.7)

The Poisson bra ket of a fun tion

F (y)

with the Hamiltonian

H(y)

is a on-venientnotationforexpressingthetime-derivativeof

F

alonga solutiontothe Hamiltoniansystem(1.5)

d

dt

F (y(t)) =

{F, H} = ∇

y

F

T

_J

∇

y

H =

∇

y

F

T

dy

dt

.

(16)

Givena Poissonbra ketandaHamiltonian,thedynami alequations(1.5) an

beba kedout

˙y

k

=

{y

k

, H

}, k = 1, . . . , 2d.

GeneralizedHamiltonian systemsareobtained byallowing

J

to bea skew-symmetri operator, dependent on the oordinates

y

, i.e.

J = J(y)

. The Poisson bra ket still takesthe form(1.7), but isrequiredto satisfyskew

sym-metry

{F, G} = −{G, F }

andthe

Jacobi identity

{E, {F, G}} + {F, {G, E}} + {G, {E, F }} = 0, ∀E, F, G.

Theserequirementstranslateintoanalogouspropertiesforthematrix

represen-tationof

J(y)

;respe tively,

J(y)

T

₌

−J(y)

and

J

im

∂J

jk

∂y

m

+ J

jm

∂J

ki

∂y

m

+ J

km

∂J

ij

∂y

m

= 0,

∀i, j, k,

whererepeatedindi esaresummedfrom

1

to

d

. TheJa obiidentityistrivially satised by a onstant

J

, so onlyskew-symmetry isneeded in this ase. Veri- ation or onstru tionof non onstantPoisson bra kets an bea ompli ated

task[49, 64 ℄.

The spe i form (1.6) of

J

orresponding to (1.4) is alled the

canonical

form of a Hamiltoniansystem. Otherwise,the systemis

noncanonical

. If

J

is dependenton

y

,thesystemis alleda

Poisson system

.

First integrals

For an ordinary dierential equation

˙y = f (y)

, a non onstant fun tion

I(y)

is a rstintegral if

I(y(t)) =

onst. alonganysolution. Inother words, ea h solutionis onstrainedtoalevelset of

I

.

d

dt

I(y(t)) =

∇

y

I(y)

T

_{f (y) = 0,}

∀y.

IntermsofthePoissonbra ket,a rstintegralis afun tionwhosePoisson

bra ketwith

H

vanishes

{I, H} = 0.

Asystemmayhavemorethenonerstintegral.Inthat asetheinitial ondition

denes an interse tion of the rst integrals and the solution evolves on this

interse tion. The set ofrst integrals foliate thephasespa e onstrainingthe

dynami s, restri ting the behavoir of a given traje tory. When one derivesa

statisti altheory,theaimisto onstru ttheleastbiaseddistributionthat still

ree ts given information about the system. This given information may be

(17)

Duetotheskewsymmetryofthebra ket,theHamiltonianis onservedfor

any(autonomous)Hamiltoniansystem

dH

dt

=

{H, H} = ∇

y

H

T

_J

∇

y

H = 0.

Inphysi alappli ationsthis often orrespondsto onservation oftotalenergy.

Other rst integrals in Hamiltonian systems arise due to ontinuous

sym-metriesintheHamiltonian,viaNoether'stheorem[64 ℄. Thatis,ifone annd

a one parameter, ontinuous hange of variables that leaves the Hamiltonian

invariant,thenNoether's theoremidentiesa orrespondingrstintegral.

Fornon anoni alHamiltoniansystems,therearesometimesfun tions

C(y)

, alled

Casimirs

, that vanish identi allyin thePoisson bra ket with any other fun tion

{F, C} = 0, ∀F = F (y)

⇔

J

∇

y

C = 0.

Casimirsareobviouslyrstintegrals,relatedto singularityof

J

. Possessionof CasimirsisapropertyofthePoissonbra ket,theydonotdependonaparti ular

hoi eoftheHamiltonian. For anoni alsystemstherearenoCasimirs.

Symplectic structure

Consider anordinarydierentialequation

˙y = f (y)

withthephasespa e

R

2d

.

Denetheowovertime

t

asamapping

φ

t

thatadvan esthesolutionbytime

t

, i.e.

φ

t

(y

0 _{) = y(t, y}

0 ₎

,where

y(t, y

0 ₎

isthesolutionofthesystem orresponding

totheinitialvalue

y(0) = y

0

. Asmoothmap

φ

t

onthephasespa e

R

2d

is alled

a

symplectic map

withrespe ttothe( onstantandinvertible)stru turematrix

J

ifitsJa obian

φ

′

t

(y)

satises

φ

′

t

(y)

T

J

−1

φ

′

t

(y) = J

−1

,

(1.8)

forall

y

in thedomain ofdenitionof

φ

t

.

Theorem (Poincar´

e, 1899).

Theowmap

φ

t

ofaHamiltoniansystem(1.5) issymple ti .

Symple ti mappingsarevolumepreserving. Takingthedeterminantofboth

sidesof(1.8)wehave

|φ

′

t

(y)

|

2 

J

−1



=



J

−1



.

Therefore

|φ

′

t

(y)

|

is either

+1

or

−1

. Sin e at

t = 0

|φ

′

0 (y)

| = 1

, a ontinuity argument showsthat

|φ

′

t

(y)

| = 1

forany

t

,i.e. thevolume ispreservedunder thesymple ti mapping.

Theowmap

φ

t

(y)

ofaPoissonsystemisa

Poisson map

,whi h satises

φ

′

t

(y)J(y)φ

′

t

(y)

T

= J(φ

t

(y)).

(1.9)

(18)

1.1.2 Geometric integration

Geometri integrators arenumeri al methods that preservegeometri

proper-ties oftheowofadierentialequation,e.g. symple ti orPoisson integrators

forHamiltoniansystems,methodspreservingrstintegrals,et . Preservationof

rstintegralsplaysanimportantroleinstatisti alme hani s,be auseoneofthe

aimsofstatisti alme hani sisto onstru taproperaverageofdesired

quanti-tiesofasystemundersu h onstraints. Symple ti itygivesvolumepreservation

oftheHamiltonianphaseow,andinaddition,a numeri alsolutiontraje tory

onashorttimes alestays losetotheexa tsolutionoftheHamiltoniansystem.

On ergodi time intervals,the solution driftsover multiple traje tories, while

preservingthephasedensityofsolutions.

Inthefollowingse tionwe onsiderpro esseswhi hevolve ontinuouslyin

time only, that is havingno ontinuous spatial stru ture. Therefore they an

be des ribed bysystems of ordinary dierential equations. In this se tion we

followthederivationsdes ribedin [35,45 ,79 ℄.

Symplectic integrators

Sin esymple ti ityisa hara teristi propertyofHamiltoniansystems

a ord-ingtotheoremofPoin aré,itistemptingtosear hfornumeri almethodsthat

sharethis property.

Aone-step numeri almethod

Φ

h

: y

n+1

_{= Φ}

h

(y

n

)

,

t

n+1

= t

n

+ h

is alled

symplectic

ifthemap

y

n

_{→ y}

n+1

issymple ti wheneverthemethodisapplied

to asmoothHamiltoniansystem.

Here aresome examples ofsymple ti methods: thesymple ti Euler rule,

the impli it midpoint rule, the Störmer-Verlet s heme, the Gauss ollo ation

methodsandotherRunge-Kuttamethodssatisfyingthesymple ti ity ondition

onthe oe ients[77 ℄.

Somesymple ti integratorspreserverstintegralsofasystem. Forexample,

impli itmidpointpreservesanyquadrati rstintegraloftheform

I =

1

2 y

T

_Ay+

b

T

_y

for a onstant symmetri matrix

A

and an ordinary dierential equation

˙y = f (y)

. Thesymple ti Euler methodand thegeneralized leapfrogmethod preserve any rst integral of the form

I = q

T

_Ap

for the anoni al variables

(q, p)

ofa Hamiltoniansystem.

Symple ti methodsdonot,ingeneral, onservetheHamiltonian

H

exa tly. Inspiteofthistheydo onserveitapproximately. Thisgood onservationof

H

insymple ti integrationisrelatedtotheexisten eandexa t onservationofa

perturbedHamiltonianasshownbyba kwarderroranalysis.

Backward error analysis

The origin of ba kward error analysis dates ba k to [87℄ in numeri al linear

algebra. For thestudy of numeri aldierential equations,its importan e was

(19)

numeri al one-step map is interpreted as the ow of a modied dierential

equation,whi his onstru tedasanasymptoti series.

Consider an ordinary dierential equation

˙y = f (y)

with exa t mapping

φ

t,f

,andaone-stepnumeri almethod

Φ

h

(y)

:

y

n+1

_{= Φ}

h

(y

n

)

,

t

n+1

= t

n

+h

. A

forward, or common, error analysis

onsistsofthestudyoftheerror

φ

h,f

(y

n

₎

−

Φ

h

(y

n

)

, i.e. the dieren e between the exa t solution and the numeri ally omputed approximation. The idea of

backward error analysis

is to derive a

modified differential equation

˙y = ˜

f

i

(y, h),

(1.10)

with

˜

f

i

(y, h) = f (y) + hf

1 (y) + h

2 f

2 (y) +

· · · + h

i

f

i

(y),

(1.11)

where the expressionfor

f

˜

i

(y, h)

depends onthe numeri almethod

Φ

h

of the original system

˙y = f (y)

byusing theTaylor expansion of theow map

φ

t, ˜

f

i

with

t = h

at

y

= y

n

and omparingitsTaylorexpansiontermswiththeterms

of thenumeri almethod

Φ

h

inpowersof

h

. Ingeneralfora one-stepmethod, theseries(1.11)divergesas

i

→ ∞

butmaybeoptimallytrun atedasafun tion of

h

.

Consider now a Hamiltoniansystem (1.5)with a smoothHamiltonian

H :

R

2d

→ R

. If a symple ti method

Φ

h

(y)

is applied to it, then the modied ve tor eld

f

˜

i

(h)

in the equation (1.10) is also Hamiltonian with a modied Hamiltonianfun tion

H

˜

i

(h)

. Morepre isely,

˜

f

i

(y, h) = J

∇

y

H

˜

i

(y, h),

andthemodiedHamiltonian

H

˜

i

(y, h)

is losetotheoriginallygiven Hamilto-nian

H(y)

,i.e.

| ˜

H

i

(y, h)

− H(y)| = O(h

p

),

where

p

≥ 1

istheorderofthesymple ti method

Φ

h

[4 ℄.

For the Hamiltonian system(1.5), themodied Hamiltonian

H

˜

i

(y, h)

an bewrittenas

˜

H

i

(y, h) = H(y) + hH

1 (y) + h

2 H

2 (y) +

· · · + h

i

H

i

(y),

where theset

{H

i

}

depends onthesymple ti method

Φ

h

as it was des ribed above,andwhereea h

H

i

isa Hamiltonianofsome Hamiltoniansystem.

Ba kwarderror analysisshowstheadvantages of usingsymple ti

integra-torsforHamiltoniansystems. Itwasshownin[35 ,45 ,79 ℄thatforasymple ti

integratorappliedtoanautonomousHamiltoniansystem,modiedautonomous

Hamiltonianproblemsexistsothatthenumeri alsolutionoftheoriginal

prob-lem is the exa t solutionof the modied problem. On the ontrary, when a

nonsymple ti integratorisusedthemodiedsystemisnotingeneral

Hamilto-niananymore. ForgeneralHamiltonians,in[33 ℄itwasproventhatasymple ti

method

Φ

h

annot exa tly onserveenergy. Nevertheless,thesymple ti inte-grators do a good job in preserving Hamiltonians approximately. Consider a

(20)

Hamiltonian system with analyti

H : K

→ R

(where

K

⊂ R

2d

is an open

subset), and apply a symple ti numeri almethod

Φ

h

(y)

withstep size

h

. If thenumeri alsolutionstaysina ompa tsubsetof

K

,thenthereexistsa

ϑ > 0

su hthat

˜

H

i

(y

n

, h)

=

H

˜

i

(y

0 , h) + O(e

−ϑ/2h

),

H(y

n

, h)

= H(y

0 , h) + O(h

p

)

over exponentiallylong time intervals

nh

≤ e

ϑ/2h

. Intypi alappli ations, the

Hamiltonianos illatesarounditsinitialvaluewithboundedamplitude

O(h

p

₎

.

In[70 ℄,itisproventhatwhenevertheowmapofagivendierential

equa-tionpossesssome geometri propertiessu hasexisten eofrstintegrals, time

reversibility, preservation of volume, symple ti ness, and the numeri al

dis- retizationpreservesthesepropertiesexa tly,thentheowmapofthemodied

dierentialequationwillalsosatisfythesegeometri properties.

Poisson integrators

Poissonintegrators generalizeofsymple ti integratorsto Poissonsystems.

Anumeri alone-step method

y

n+1

_{= Φ}

h

(y

n

)

,

t

n+1

= t

n

+ h

, isa

Poisson

integrator

for the stru ture matrix

J(y)

, if the transformation

y

n

→ y

n+1

respe ts theCasimirsand ifitis a Poisson map (1.9)whenever themethod is

applied tothePoisson system.

A ordingtoba kwarderroranalysis,ifaPoissonintegrator

Φ

h

isappliedto thePoissonsystem,thenthemodiedequationislo allyaPoissonsystem[35 ℄.

There is no general te hniquefor onstru ting Runge-Kutta type Poisson

methods. Themostgenerallyappli ablealternativeissplittingmethods,e.g.[50 ℄.

1.2 Hamiltonian fluid dynamics

Thestudyofthedynami sofuidsisoneofthemostattra tiveareasinapplied

mathemati s. Thefa tthatuiddynami sisanattra tiveresear hareaisdue

tomanyreasons. Perhapsthemostimportantoneistheintrodu tionofe ient

highresolutionnumeri alsimulationsintouiddynami sasaresear htool. The

signi an e of thistoolis espe ially pronoun edin ase of omplexbehaviour

ofa system.

1.2.1 Hamiltonian PDEs

The dynami s of uids has both propagation in time and a spatial stru ture

and,hen e, annotbedes ribedbyordinarydierentialequations(ODEs)

any-more but by partial dierentialequations (PDEs). Many PDEs that arise in

physi s an beviewedas innite-dimensionalHamiltoniansystems. (Problems

des ribed by ordinary dierential equations are nite-dimensional.) The

(21)

PDEsismu hlessexploredthanthatofODEs,sin ethesolutionbehaviourof

PDEsismu hmore omplex. Nevertheless,thereexistsanumberofpapersfor

HamiltonianPDEs[12 ,14 ,17 ,24 ,26 ,27 ,31 , 32 ,51 ,83 ℄.

Thenite-dimensionalHamiltoniansystem onsistsofatriple(

K,

{·, ·}, H

), where the phasespa e

K

⊂ R

2d

is anopen subset,

H : K

→ R

is the Hamil-tonian fun tion, and

{·, ·}

is a Poisson bra ket with stru ture matrix

J(y)

, see (1.7). When the phase spa e is innite-dimensional, we write the triple

as

(

K, {·, ·}, H)

,andthePoissonoperator as

J

,a ordingto [64 ℄. Typi ally

K

onsistsofsetsofsmoothfun tionsonanite-dimensionalspa e

Y

. Anelement in

K

isdenoted by

u(y)

,

y

∈ Y

. TheHamiltonian

H : K → R

is a fun tional onthisspa e,andthebra ket anbewrittenas

{F, G}[u] =

Z

Y

δ

_F

δu

J (u)

δ

_G

δu

dy,

where

δ

F/δu

isthevariationalderivativedenedby

lim

ǫ→0

F[u + ǫδu] − F[u]

ǫ

≡

Z

Y

δ

_F

δu

δu dy

forappropriate

δu

.

J (u)

is,ingeneral,adierentialoperator, alledthe

Poisson

operator

.

Motivated by thesu ess of symple ti integrators, a reasonable approa h

to HamiltonianPDEsisto tryto dis retizein spa ewhilepreservingthe

sym-ple ti or Poisson stru ture. For anoni al stru ture, it is a simple matter to

dis retizetheHamiltonian fun tionalwith anydesiredquadrature. Theresult

is a Hamiltonian ODE to whi h symple ti integrators may be applied. For

Poisson systems, it is a signi ant hallenge to derive a dis rete bra ket that

preservestheJa obiidentity. IfthereareCasimirs,there shouldbesome

rem-nantofthese. Foruids,thereisaninnitefamily,andonlyanitenumberof

independent integrals an survive, ifthey onstrain thedis rete, nite

dimen-sional phasespa e.

Onthe other hand,sin e Lagrangian uid dynami s is anoni al, one an

approximatethePDEsolutionwithasetofmovingparti lesintera tingthrough

anappropriatepotentialenergyfun tion,andaHamiltoniansemi-dis retization

will be obtained for any quadrature of the Hamiltonian. The set of

nite-dimensionalHamiltonianODEsisthenintegratedintimeusingasuitable

sym-ple ti orPoissonintegrator.

Unfortunately, inthe aseof aPoisson PDE,unlike anoni alHamiltonian

PDEs, it is not possible to establish a ommon generi approa h. For ea h

parti ular problemone hastodevelopa proper wayof redu ingthePDEtoa

systemofODEs. Wewill onsiderthequasigeostrophi potentialvorti ity

equa-tion (a Hamiltonian PDE with Poisson stru ture), for whi h we will des ribe

(22)

1.2.2 Eulerian and Lagrangian descriptions

Geophysi al uiddynami sisthestudyofuidmotionintheatmosphereand

theo ean. Therststepinthisstudyisa hoi eofframework: Eulerianor

La-grangian. The

Eulerian description

is ommonlyusedinliterature,e.g.[52 ,67 ℄, ittreats motionoftheuidasaeld inwhi h thevelo ityistobedetermined

at allpositions and times. The

Lagrangian description

regards the uid as a ontinuouseldof parti les,whosepositionsareto bedetermined[74℄.

IntheEuleriandes ription,theindependentvariablesarethespa e

oordi-nates

x

= (x, y, z)

∈ D

and the time

t

. Thedependent variables in lude the velo ity

v(x, t)

andthemassdensity

ρ(x, t)

Thenon anoni alPoisson bra ket foridealuidowinEulerian variables,

a ording to[58 ℄,is

{F, G} = −

Z

D

dx

· δF

δρ

∇

x

· δG

δv

−

δG

δρ

∇

x

· δF

δv

+

µ ∇

x

× v

ρ

¶

· µ δG

_δv

×

δF

_δv

¶¸

withHamiltonian

H[v, ρ] =

Z

D

dx

· v · v

2ρ

+ ρE(ρ)

¸

.

Here

x

is anindependentvariable.

ThePoisson bra ket has aninnite lassof potentialvorti ityCasimirs of

theform

C[ρ] =

Z

D

dx ρf

µ ∇

x

× v

ρ

¶

foranarbitraryfun tion

f

.

In the Lagrangian des ription, ea h uid parti le is assigned a label

a

=

(a, b, c)

∈ A

. Forexample,thelabelsmaybedenedasthepositionsofparti les at the initial time. The independent variables are set of

a

, whi h are xed for ea h parti le, and the time

t

. The dependent variables are the position oordinates

x(a, t)

. Thevelo ityofa parti leisgivenby

v

=

µ ∂x

∂t

,

∂y

∂t

,

∂z

∂t

¶

.

Themassdensity

ρ

isdenedvia Ja obianmatrix

|

∂

x

∂

a |

as

ρ = ρ

0 



∂x

∂a



−1

,

where

ρ

0 = ρ

0 (a)

doesnotdepend ontime

t

.

Taking derivativesoftheexpressionaboveleadsto the ontinuityequation

dρ

(23)

ThePoissonbra ketforidealuidowin Lagrangianvariablesis anoni al

{F, G} =

Z

A

da

· δF

δx

· δG

δv

−

δG

δx

· δF

δv

¸

withHamiltonian

H[x, v] =

Z

A

da

"

v

_{· v}

2ρ

0 + ρ

0 E

Ã

ρ

0 



∂x

∂a



−1

!#

.

Here

x

is adependentvariable.

The quasigeostrophic model

The two-dimensional quasigeostrophi potential vorti ity (QG) equation [48 ,

67 ,74 ℄des ribesdivergen e-freeowovertopographyby

d

dt

q = 0,

∆ψ(x, t) = q(x, t)

− h(x),

(1.12)

where

q

isthepotentialvorti ity(PV)eld,

ψ

isthestreamfun tion,and

h

is thetopographyoftheearth. TheLapla ian operatoris denotedby

∆

andthe material derivativeby

d

dt

=

∂

∂t

+ u

· ∇

. Here, thedivergen e-freevelo ityeld

u

isrelatedtothestreamfun tionby

u

=

∇

⊥

_ψ

,where

∇

⊥

_{= (}

−

∂y

∂

,

∂

∂x

)

T

. We

onsider theQG equationona doublyperiodi domain

x

= (x, y)

∈ D ≡ [0, 2π) × [0, 2π).

Dene the operator

J (q, ψ) = q

x

ψ

y

− q

y

ψ

x

. The QG model des ribes a HamiltonianPDEwithPoissonstru ture[57 ℄,

{F, G} =

Z

D

q

_J

µ δF

δq

,

δ

G

δq

¶

dx,

implyingthe onservationoftheHamiltonianortotalkineti energy

H = E = −

1 ₂

Z

D

ψ

_{· (q − h) dx}

as wellastheinnite lassofCasimirfun tionals

C[f] =

Z

D

f (q) dx

for any fun tion

f

for whi h the integral exists. Of spe i interest are the momentsofPV

C

r

=

Z

D

q

r

dx,

r = 0, 1, 2, . . . ,

(1.13)

(24)

PreservationoftheCasimirfun tionalsfollowsfromarea-preservationunder

the divergen e-free velo ity eld [55 ℄. Dene a fun tion

G(σ, t)

denoting the measure ofthatpartofthedomain

D

forwhi h thevorti ityislessthan

σ

:

G(σ, t) = meas

_{{x ∈ D | q(x, t) < σ}.}

We note that due to the divergen e-free adve tion of

q

, this fun tionis inde-pendentoftime,

∂G

∂t

= 0

. Dierentiatingwithrespe t to

σ

,thefun tion

g(σ) =

∂G

∂σ

(1.14)

ispreserved. Forthe aseofapie ewiseuniformPVeld,

q(x, t)

∈ {σ

1 , . . . , σ

L

}

, thisquantity

g

ℓ

= G(σ

ℓ+1

)

− G(σ

ℓ

)

isthemeasure ofthevorti itylevelset

σ

ℓ

.

1.2.3 Numerical methods

In the following se tion we des ribe several numeri al methods to solve the

quasigeostrophi model(1.12).

The Zeitlin method

Normalspe tralmethods fortheQG equationpreservetheenergyand

enstro-phy at most. However, the Zeitlin method [88 ℄ is a spe tral method whi h

preservesa Poisson stru ture, theHamiltonian and

2M

Casimirsin a

(2M +

1)

× (2M + 1)

modetrun ation.

Therstequationin(1.12)istransformedthroughtwo-dimensionalFourier

seriestakingtheform ofaninnitesystemofODEs

dˆ

q

k

dt

=

∞

X

k

′

1 ,k

′

2 =−∞

k

′

6=0

k

_{× k}

′

|k

′

|

2 q

ˆ

k+k

′

(ˆ

q

−k

′

− ˆh

−k

′

).

(1.15)

Here

q

ˆ

k

denotes the spe tral oe ient asso iated with the two-dimensional wave ve tor

k

, whose omponents are integers. The skew-symmetri s alar produ t

k

× k

′

is

k

1 k

′

2 − k

2 k

′

1

, andthe norm

|k|

is

pk

2

1 + k

2

. Sin e

q

is real,

ˆ

q

∗

k

= ˆ

q

−k

.

Zeitlin proposed the sine-bra ket trun ation of the equations. The

nite-dimensional setofequationsfortheFourier oe ientsisthengivenby

dˆ

q

k

dt

=

M

X

k

′

1 ,k

′

2 =−M

k

′

6=0

1 ǫ

sin(ǫk

× k

′

)

|k

′

|

2 q

ˆ

k+k

′

(ˆ

q

−k

′

− ˆh

−k

′

),

ǫ =

2π

2M + 1

,

(1.16)

where all indi es are redu ed modulo

2M + 1

to the periodi latti e

−M ≤

k

1 , k

2 ≤ M

. The summation o urs on the

(2M + 1)

× (2M + 1)

domain of the Fourier oe ients. For

M

→ ∞

and given

k

and

k

′

,

ǫ

−1

_sin(ǫk

× k

′

) =

k

_{× k}

′

+ O(ǫ

2 ₎

(25)

Thistrun ationpossessesaHamiltonianstru turewithsymple ti operator

J

kk

′

=

−ǫ

−1

sin(ǫk

× k

′

)ˆ

q

_k+k

′

,

andHamiltonian

H = E =

1

2 X

k6=0

|k|

2 | ˆ

ψ

k

|

2 =

1

2 X

k6=0

|ˆq

k

− ˆh

k

|

2 |k|

2 .

(1.17)

Thesymple ti matrixisskew-symmetri andsatisestheJa obiidentity. The

sine-bra ket trun ation (1.16)preservesthe Hamiltonian(1.17) and

2M

inde-pendentCasimirinvariants orrespondingtotherst

2M

momentsofpotential vorti ity. If in the Zeitlin method thePoisson dis retization is integrated

us-ingthePoissonsplittingofM La hlan[50 ℄,thenthesequantities arepreserved

by the splitting (the energy is only preserved approximately, in the sense of

ba kwarderroranalysis[35 ℄).

Unfortunately, theZeitlin method islimitedto 2Din ompressible owson

periodi geometry.

Arakawa’s scheme

TheZeitlinmethodistheonlyknowndis retizationwithPoissonstru turefor

Eulerian uid models. For more general uid problems ( ompressible,

non-periodi boundary onditions, et .) no Poisson dis retizations are available.

How anwepreserveatleastsomequantities?Awellknowns hemeisArakawa's

s heme[2℄whi hpreserveslinearandquadrati invariants.

For astartwerewritetherstequationin (1.12)as

q

t

=

J (q, ψ),

(1.18)

where theoperator

J

isdened by

J (q, ψ) = q

x

ψ

y

− q

y

ψ

x

.

Arakawa'sidea onsistsofseveralsteps. Firstofall,heuses entraldieren es

for

x

-and

y

-derivatives. Thenherewritesthe ontinuous

J

inthreeequivalent formsbasedonthefa tthatthederivativeswithrespe tto

x

and

y

ommutein the ontinuous ase,namely

J (q, ψ) = ∂

x

(qψ

y

)

− ∂

y

(qψ

x

) = ∂

y

(q

x

ψ)

− ∂

x

(q

y

ψ)

. After dis retizing these three equivalent forms of

J

and taking theiraverage, one gets four dis rete non-equivalentright-hand sidesof (1.18), therefore four

dis retizations. Non-equivalen eofthedis reteright-handsidesisexplainedby

the fa t that theprodu t rule does nothold anymore in the dis rete ase. It

an beshownthat one dis retization does not onserveanythingand, in fa t,

is unstable;se ond one onservesonlyenergy; thirdoneonly enstrophy; and

fourthone,whi hisanaverageofthepreviousthree, onservesbothenergyand

enstrophy. Itisworthmentioningthatallthesedis retizationsarealsovolume

preserving in the sense of the Liouville property, see Se tion (1.3). This is a

ne essary ingredientforastatisti altheory.

(26)

The Hamiltonian particle-mesh method

The Hamiltonian parti le-mesh (HPM) method approximates the solution of

an ideal uid ow using a set of moving parti les that intera t through an

appropriate potential energy fun tion. The HPM method was originally

de-signed for rotating shallow water owwith periodi boundary onditions [25 ℄

andextendedtootherphysi alsettingsrotatingtwo-layershallowwatermodel

withrigid-lid onstraint,barotropi model,non-hydrostati verti alsli emodel,

in [16 ,17,28,83℄.

HPM isbased onthe Lagrangian formulation of uiddynami s, anda set

of moving point masses ombined with a xed Eulerian grid. The spatially

trun atedequationsare anoni alHamiltonianandsatisfyaKelvin ir ulation

theorem[25 ℄. Implementationofasplittingmethodintimetothesemi-dis rete

system of Hamiltonian ODEs gives symple ti ity. Hen e the HPM method

is symple ti , and one an onstru t a ontinuum velo ity eld in whi h the

dis rete parti le velo ities areembeddedfor alltime. Energy is preserved

ap-proximatelyinthesenseofba kwarderroranalysis,i.e. goodlong-timeenergy

onservation. Convergen eofthemethodwas onsideredin[56 ℄.

Theideaofxingthepotentialvorti ityratherthanthemasstoa parti le

wasoriginallyproposedinthis ontextfortwo-dimensionaladve tionunder

in- ompressibleoweldsin[17 ℄. Theresultisaregularizedpointvortexmethod.

Therethepotentialvorti ityofaquasigeostrophi modelissimplyadve tedin

a divergen e-free velo ity eld, obtained by re onstru ting the PV eld on a

uniformgrid. Theparti lemotion anbeembeddedin anarea preservingow

ontheuidlabelspa e. Hen e,theCasimirs(1.13)are trivially onserved ifa

valueofpotentialvorti ity issimplyassignedtoea hparti leon e andforall.

The semi-dis rete system is still Hamiltonian, and an be integrated in time

usinga symple ti integrator. Sin etheHPMmethod issymple ti ,thephase

owisvolumepreservinginthesenseoftheLiouvilleproperty,seeSe tion(1.3),

whi hisa ne essaryingredientforastatisti altheory.

TheHPMmethodforthequasigeostrophi modelisexplainedinmoredetail

in Chapter3.

1.3 Statistical mechanics of fluids

Statisti alme hani sisapowerfultoolforunderstanding omplexphysi al

sys-tems, e.g. [23,30 ℄. There aredierent purposesof statisti al me hani s. One

may onsider a systemfor a very long time. Then the quantity of statisti al

interest is the average behaviour of a system rather than the behaviour of a

system at a ertain time. Anexample is a tra outside of your window. A

statisti al quantityof interestmaybetheaveragespeedof ea h ar. Another

purposeofstatisti al me hani sisto designmethodsto handlesystems whi h

arein ompletelyknown. For example,ifwedonotknowtheinitial onditions

of a system,and wewantto know itsmostlikelybehaviour. Thisinvolvesan

(27)

Statisti alme hani sisbroadlyemployedinnumeri alsimulationsof

mole -ular dynami s and limate variability. Here either a systemis very big (

10

23

mole ules)orlongtimesimulationsare haoti ,orboth. Moreover,dierent

nu-meri al dis retizationsofadynami alsystemhavedierentdis retedynami s.

Thereforethestatisti alresults analsobedistin t.

Themainpremiseoftheequilibriumstatisti altheoriesofeither ontinuous

or dis rete dynami s is ergodi ity relating the equilibrium distribution to the

dynami s. Ne essaryingredientsare onservationlawsandvolumepreservation

ofthephaseow. Therefore,toderiveastatisti altheoryofanumeri almethod

thenumeri almethodhastopossess onservationlaws,and thedis retephase

owshould bevolumepreserving. Whenoneintegratesa Hamiltoniansystem

with a symple ti integrator, there is an automati onservation of dis rete

analoguesoftheexa t onstantsofmotionand,be auseofsymple ti ity,there

isvolumepreservation. Ergodi ityisdi ulttoshowfornontrivialsystems,so

thisisusuallytakenasanassumptionor'approximation'whi hmustbeveried

withnumeri s.

Theaimofthefollowingse tionistoexplainstatisti alme hani s,togivea

denitionofdierentensemblesandpurposeofea h,andtodis ussthe on ept

ofergodi ity, entropyandtemperaturein statisti alsense.

Statistical equilibrium

Consider anordinarydierentialequation

˙y = f (y)

withthephasespa e

K

⊂

R

2d

,thenaprobabilitydistributionfun tion

ρ(y, t)

,

ρ : K

×R → R

,forexample, over a set of un ertain initial onditions,is transported bythe ow a ording

to

∂

∂t

ρ +

∇

y

· ρf = 0.

Now onsider aHamiltoniansystem(1.5).

Theorem (Liouville, 1838).

The phaseow ofa Hamiltonian system(1.5) isvolumepreserving[3℄.

For onstant

J

this followsfrom theskew-symmetry of

J

, or equivalently, thedivergen e-freenatureofthe anoni alphaseow

y

. Onesaysthattheow hasthe

Liouville property

,ifthisowsatisesLiouville'stheorem,i.e. theow is volume preserving, i.e. the ow is divergen e-free. The Liouville property

is a ne essary ingredient fora statisti al theory, sin efrom itfollowsthat the

probability measure is transported under the divergen e-free ow. Therefore

before deriving a statisti al theory, one has to prove that the phase ow is

divergen e-free.

Thetransportequationinthe aseofa Hamiltoniansystemsimpliesto

∂

∂t

ρ + J

∇

y

H

· ∇

y

ρ = 0,

(28)

A steadystateof theLiouville equation forthe Hamiltonianow,

∂ρ

∂t

= 0

,

isthen

J

∇

y

H

· ∇

y

ρ = 0,

whi hisoftenreferedas an

equilibrium

probabilitydistributionfun tion. Note that anyfun tion

ρ(y) = ρ(H(y))

oftheHamiltonianisanequilibrium proba-bilitydistribution,andtondtheproperprobability,whi h orrespondstothe

dynami sdes ribedbytheODE,isthetopi ofergodi theory,seee.g.[68 ,85 ℄.

To bemorepre ise, letusdene thetimeaverageofafun tion

F (y(t))

as

F

≡ lim

T →∞

1 T

Z

T

0 F (y(t))dt,

providedthatthelimitexists.

Theensembleaverageof

F

isdened as

hF i ≡

Z

K

F (y)ρ(y) dy

_≡

Z

K

F ν(dy)

fora propermeasure

ν

su h that

ν > 0

and

R

K

ν(dy) = 1

.

Ergodicity

implies thatthelong time average isequivalentto theensemble average

F =

hF i

(1.19)

withtheprobabilitymeasure

ν

ora reasonableapproximationtoit. Giventhe probability measure

ν

it is a hallenging task to prove theergodi ity, and we will takeitasanassumption.

1.3.1 Statistical ensembles

Let usintrodu eanimportantideaof mi rostateand ma rostate ofa system.

Consideramodelinwhi huidmotionisdes ribedbyasetofmovingparti les.

Thena

microstate

ofsu hsystem anbedes ribedbypositionsoftheseparti les, anda

macrostate

is, forexample,theobservableenergywhi hmay orrespond to alargenumberofmi rostates.

Belowwewillexplainthisinmoredetailfollowingthederivationsdes ribed

in [10 ℄.

The microcanonical ensemble

Consideradis retespa e

K

withasinglema rostategivenbytheenergy

H = E

. Let

y

∈ K

. Thesubset

D(E) =

{y ∈ K : H(y) = E}

onsistsofdis retestates

y

withthesameenergy

E

.

Dene

Ω(E)

to be the total number of states

y

∈ K

with the energy

E

. The mi ro anoni alensembleis the set of all

y

having

H(y) = E

. Assuming

(29)

all su h states are equally likelywe dene the

microcanonical density

for the dis retespa e

Prob

_{{y | H = E} =}







0,

H(y)

6= E

1/Ω(E),

H(y) = E.

Themi ro anoni alentropyisdened as

S(E) = ln Ω(E).

ConsidernowajointHamiltoniansystem

AB

,whi h onsistsoftwosystems

A

and

B

,withtotalenergy

E

su hthat,the ouplingelementallowsex hange of the energy between systems

A

and

B

and adds neither new states to the systemnornewtermsintotheHamiltonian.

A mi rostateof system

A

is dened as

y

A

. If system

A

isin a mi rostate

y

_A

with orrespondingenergy

E

A

,thensin ethereistotalenergy onservation, system

B

shouldbeinastatewiththeenergy

E

− E

A

. Theprobabilityofthe mi rostate

y

A

is

Prob

_{y

A

| H = E} =

Ω

B

(E

− H

A

(y

A

))

N (E)

,

(1.20)

where

N (E)

isanormalization onstantsu h that

X

y

A

∈K

A

Prob

_{y

A

| H = E} = 1

⇔

N (E) =

X

y

A

∈K

A

Ω

B

(E

− H

A

(y

A

)).

This is the mi ro anoni al probability of a mi rostate. The mi ro anoni al

probabilityofthema rostate

H

A

= E

A

istheenergysplit

Prob

{H

A

= E

A

| H = E} =

Ω

A

(E

A

)Ω

B

(E

− E

A

)

N (E)

withthenormalization onstant

N (E)

.

Themostprobable ma rostate anbefoundbymaximizing thenumber of

stateswiththeenergysplitoverallpossiblestatesofsystem

A

withtheenergy

E

A

max

E

A

[Ω

A

(E

A

)Ω

B

(E

− E

A

)] ,

whi hin termsofentropiesgives

max

E

A

[S

A

(E

A

) + S

B

(E

− E

A

)] .

Fromthelastexpressionitis learthatthemostprobablestateisthemaximizer

ofthetotalentropy. Whenthejointsystem

AB

isverylarge,thenitispossible tomaketransitionfromthedis retephasespa etoa ontinuousone,andthen

(30)

the maximizer of the total entropy an be found analyti ally: the maximum

o urs atsome

E

∗

A

where

S

A

′

(E

∗

A

) = S

B

′

(E

− E

∗

A

),

and the prime denotes thederivative with respe tto theargument. This has

motivatedthedenition ofthe

microcanonical statistical temperature

S

′

(E) =

1 T

.

Note that

T

A

=

T

B

at the mostlikelyma rostate. Inthe literatureone often usesthe

inverse statistical temperature

β = 1/

T

.

Now we show how to derive the mi ro anoni al density for a ontinuous

phasespa e. Therestofthemi ro anoni alstatisti altheoryfora ontinuous

phasespa efollowsautomati allywithsumsrepla edbyintegrals.

Iftheenergyistheonly onservedquantityofthesystem,thenwe onsider

thesubspa e

D(E, dE) =

{y ∈ A : H(y) ∈ [E, E + dE]},

(1.21)

with orrespondingdensity

ρ(y) =







0,

H(X) /

_{∈ [E, E + dE]}

1/vol

{D}, H(X) ∈ [E, E + dE].

(1.22)

The density (1.22) is a stationary density, sin e for xed

(E, dE)

it depends onlyonautonomous

H

. Ifwetakethelimit

dE

→ 0

,thedensity

ρ

ispresented onlyonthesurfa e

H = E

. Thenthemi ro anoni aldensityforthe ontinuous phasespa e anbewrittendownin termsofDira deltafun tions

ρ(y) =

1 Ω(E)

δ(H

− E)

with

Ω(E) =

Z

K

δ(H

− E) dy,

where

Ω(E)

isthemeasureofthesurfa e

H = E

.

We want to underlinethat the mi ro anoni alstatisti al me hani s is

de-rivedassuming onservationofsome quantities,e.g. energy. Thereforeall

pos-sible stateswith these onstantquantities,i.e. the ma rostate, determinethe

mi ro anoni alensemble.

The canonical ensemble

Consider again thejoint system

AB

, but withthe size of thesystem

B

mu h largerthenthesizeofthesystem

A

,thesizeofthesystem

A

mayormaynotbe large omparedtounity. Inthis asethesystem

B

is alledanenergyreservoir forthesystem

A

.

Nowwewouldliketoderiveanequivalen eof(1.20)forthedes ribedsystem.

First,letuswritedown(1.20)in termsofentropy:

(31)

Sin ethefun tion

S

isa slowlyvaryingfun tionin anamplerangeofpossible mi rostates ofthelargesystem

B

,in ontrastto

Ω

whi h isnot, we an write downtheTaylorexpansionof

S

B

(E

− E

A

)

,andwetrun atetheseriesafterthe rstterm

S

B

(E

− E

A

)

∼ S

B

(E) + S

B

′

(E)(E

− E

A

− E).

Absorbing

S

B

(E)

into the onstantof normalizationandnoti ing that

S

′

B

(E)

istheinversetemperature

β

B

,weobtain

Prob

{y

A

| H = E} ∼ exp(−β

B

E

A

).

Thismotivatesthedenition ofthe

canonical probability density

ofasystemin onta twithanenergyreservoirin aseofa dis retephasespa e

Prob

_{{y | β} =}

1 N (β)

exp(

−βH(y))

and

N (β) =

X

y∈K

exp(

_{−βH(y)). (1.23)}

For a ontinuousphasespa e

K

, the sumin (1.23) isrepla ed bytheintegral over

K

.

Sampling

Nowthatwehavedenedthestatisti alensemblesandtheprobabilitydensities

asso iatedwiththem,aquestionariseshowtoensuresamplingofthese

ensem-bles. If we onsider a system with onserved energy, and having ergodi ity,

we an samplea mi ro anoni aldistribution of onstant energybysimulating

the dynami sof thesystemfor a longtime. Then theensemble originatedby

the time series of the dynami s is equivalent to the mi ro anoni al ensemble

asso iatedwiththe onstantenergy.

Thereareseveralapproa hestoensuresamplingofa anoni aldistribution.

Theseapproa hesworkinsu hawaythatiteithermodiesthedynami al

sys-temorintrodu esasto hasti perturbation. Themostknown lassi almethod

for sampling a anoni al distribution is the Metropolis algorithm [53℄. It is

basedonarandom hoi eofa stateanda eptan eofthisstatedependingon

theprobabilitywhi hshouldbesampled. Apopularmethodologyinmole ular

dynami s to samplea anoni al distribution isa mathemati al thermostata

tooltomodelthesysteminthermalequilibriumwithareservoir. The

thermo-statis responsiblefortheenergyex hange betweenthesystemandtheenergy

reservoirsu hthatthesystemstaysata giventemperature,whi hfor es

sam-pling of the anoni al equilibrium distribution. Here are several thermostat

te hniques. The lassi alonesareLangevindynami s[80℄,whi hisasto hasti

thermostat,anddeterministi thermostatssu hastheNosémethod[61 ,62 ℄and

theNosé-Hoover method [37 ,62 ℄. InLangevindynami s the ombinationofa

damping for e anda sto hasti termmaintainsthesystem ata given

temper-ature. Be ause of thepresen e of damping and theintrodu tion of a random

for ing, the dynami s are not any longer Hamiltonian. The Nosé and

Nosé-Hoover methods preserve the Hamiltonian stru ture and a hieve sampling by

(32)

1.3.2 Information theory

Wehave onsideredthemaximumentropyprin iplebasedonthemaximization

of the number of states of a system and repla ing this number by the total

entropy. There is an alternative approa h to the maximum entropy prin iple

based on Shannon entropy of information theory, whi h is not derived from

physi alprin iples.

Consideraninnitedimensionalphasespa e

K

(it ouldalsobenite)with an element of it denoted by

y

. Then the

Shannon entropy

, or information entropy,is

S[ρ] =

−

Z

K

ρ ln ρ dy.

(1.24)

The on ept of Shannon entropyplays the entral role in information theory,

sometimesreferredasmeasureofun ertainty[81 ℄. Theprobabilitydensity

fun -tion

ρ

is hosentomaximize

S

under onstraints orrespondingtoobservations on the system. These minimal assumptions lead to the distribution of least

bias, i.e. the distribution whi h is most generaland still explains the

obser-vations. Weshow that themi ro anoni aland anoni aldistributions arethe

maximizersofShannonentropyundersuitable onstraints.

Supposethatwehaveonlyone onstraintontheprobabilitydensity,namely

thenormalization onstraint

Z

K

ρ dy

_{− 1 = 0.}

(1.25)

To ensurethis onstraintwe needto in lude it in thea tion prin iple via

La-grange multiplier

θ

δ

· −

Z

K

ρ ln ρ dy

_{− θ}

µZ

K

ρ dy

_{− 1}

¶¸

= 0

forarbitraryvariations

δρ

and

δθ

. Aftertakingvariationswithrespe tto

ρ

we have

ln ρ =

−1 − θ

with

θ

determinedby thenormalization onstraint(1.25). Therefore

ρ =

1 vol

{K}

,

and

S

∗

_{= ln vol}

{K}

is themaximized entropy. Thus,withnofurther assumptions,theleastbiased

distribution isuniform.

Considerthephasespa e(1.21)withthemi ro anoni alprobabilitydensity

(1.22). Then it an be shown that this density is the maximizer of Shannon

entropyunderboththenormalization onstraint(1.25)andthe onstraintthat

theenergy anonlytakevalues

H

∈ [E, E +dE]

. The orrespondingmaximized Shannon entropyis

S

∗

_{= ln vol}

_{{D(E, dE)}}

.

The anoni aldistribution (1.23) is the maximizer of Shannon entropy as

well, but under other onstraints,namely, under the normalization onstraint

(1.25)and the onstraintofobservedmeanenergy

U

Z

K

(33)

forsome xed

U

.

Thea tionprin iple statesthat

δ

· −

Z

K

ρ ln ρ dy

_{− θ}

µZ

K

ρ dy

_{− 1}

¶

− β

µZ

K

Hρ dy

_{− U}

¶¸

= 0

forarbitrary

δρ

,

δθ

and

δβ

. Here

θ

and

β

areLagrangemultipliers orresponding to thenormalizationandthexedmeanenergy onstraintsrespe tively.

Taking variationswithrespe tto

ρ

weobtain

ρ =

1 exp(1 + θ)

exp(

−βH).

Sin e

θ

here orrespondsto the normalization onstraint, after renamingitas

N

,itbe omes learthat thisisthe anoni aldistribution (1.23).

Thereforethemi ro anoni aland anoni aldistributionsarethemaximizers

of the Shannon entropy under suitable onstraints. This an be extended to

a generalsystemwith more onstraintsmore information about thesystem.

For example, if one onsiders a numeri al method with some quantities, say,

preservationof energy, or anyother information,then thisinformation anbe

usedto onstru ttheleastbiaseddensity onsistentwiththeobservations.

The ee tiveness of a density derived this way depends on the detail to

whi h known information about the system is in luded. For example, the

energy-enstrophystatisti altheoryforthequasigeostrophi model(1.12)based

onpreservationofenergyandenstrophyisamodelandisin omplete[1℄,sin e

ittakesintoa ountonlylinear( ir ulation)andquadrati (enstrophyand

en-ergy) invariants,while weknow that thequasigeostrophi model preservesan

innitenumberofCasimirs(1.13).

Anothermeasureofinformational ontentis

relative entropy

,alsoknownas theKullba k-Leibler'sdistan eor thedivergen e [41℄

S[ρ, Π] =

₋

Z

K

ρ ln

³

ρ

Π

´

dy,

(1.27)

where

Π

is a probability density over

K

representing an externalbias due to someadditionalinformation. Forexample,inthequasigeostrophi model(1.12)

there areaninnitenumber ofCasimirs(1.13). Thismeansthat formallyone

has to onsider an innite number of onstraints on the entropy. To nd a

solution to all these onstraints might be a di ult or even impossible task.

Insteadofthis,

Π

anbe hosensu hthatitree tstheCasimirs,andtherefore itgivesanexternalbiasonthespatialdistributionofPV.

It an be shownthat

S

is non-positive, and

S

is zeroonly if

ρ

≡ Π

every-where. Thisexplainsthetermdistan efor(1.27). Thedensity

Π

isoften alled a

prior distribution

andin atypi alappli ation

Π

isgivenand

S

ismaximized overpossible hoi esof

ρ

.

Statistical mechanics and numerical modelling of geophysical fluid dynamics - Thesis

UvA-DARE is a service provided by the library of the University of Amsterdam (http

s

://dare.uva.nl)