Transcript Geometric Integration of Differential Equations 1. Introduction and ODEs Chris Budd Want to simulate a physical system governed by differential equations Expect the numerical.

Geometric Integration of Differential Equations
1. Introduction and ODEs
Chris Budd
Want to simulate a physical system governed
by differential equations
Expect the numerical approximation to have the
same qualitative features as the underlying solution
Traditional approach
 Carefully approximate the differential operators in the system
 Solve the resulting difference equations
 Monitor and control the local error
Basis of most black box codes and gives excellent results
over moderate computing times
BUT This is a local process and does not pay attention to
the qualitative (global) features of the solution
Geometric Integration
Aims to reproduce qualitative and global features
Some qualitative properties:
Conservation laws
Global quantities: Energy, momentum, angular momentum
Flow invariants: Potential vorticity, Casimir functions
Phase space geometry
Symplectic structure
Symmetries
Galilean
Reversal
Scaling: Nonlinear Schrodinger
Lie Group: SO3 (Rigid body)
Some global features
Asymptotic behaviour
Orderings
Often linked: Noether’s theorem for Lagrangian flows
Example: The Kepler Problem


2
2
d
x
x d
y
y


,


2
2 2
3
/
2 2
2 2
3
/
2
dt
(
x

y
)
dt
(
x

y
)
Conserved quantities:

1
22
H

(
u

v
)

,L

xv

yu
2 2
1
/
2
2
(
x

y
)
Hamiltonian
Angular Momentum
Symmetries: Rotation, Reflexion, Time reversal, Scaling
 

2
/
3

1
/
3
t

t
,(
x
,
y
)

(
x
,
y
),
(
u
,
v
)

(
u
,
v
)
Kepler's Third Law
Geometric Integration
 Aims to preserve a subset of these features
 Take advantage of powerful global error estimates
(shadowing)
 Powerful methods for important physical problems
Examples of GI methods
Symplectic and multi-symplectic
Splitting
Lie Group/Magnus
Discrete Lagrangian
Scale Invariant
Examples of GI applications
Molecular and celestial mechanics
Rigid body mechanics
Weather forecasting
Integrable systems (optics)
Self-similar PDEs
Highly oscillatory problems
Example of the traditional and the GI approach:
Integrating the Harmonic Oscillator
dx dy
y, 
x
dt dt
Qualitative features: Bounded periodic solutions, time reversal symmetry,
1 2 2
H (x y)
2
Conserved
Forward Euler method (non GI)
X

X

hY
,
Y

Y

hX
n

1
n
n
n

1
n
n
12 2
2
H

X

Y
)

H

(
1

h
)
H
n (
n
n
n

1
n
2
Problem: Energy increases, lack of periodicity, lack of symmetry
Backward Euler Method (non GI)
X n1  X n  hYn1, Yn1  Yn  hXn1
2
H

H
/(
1

h
)
n

1
n
Problem: Energy decreases, lack of periodicity, lack of symmetry
Mid-point rule (a GI method)
h
h
X n 1  X n  (Yn  Yn 1 ), Yn 1  Yn  ( X n  X n 1 )
2
2
Hn1  Hn
FE
Mid-Point rule
BE
Mid point rule conserves:
Energy
Symmetry
Backward (Modified) Equation Analysis
Discrete equation has an exact solution
h2
X n  cos(tn ), Yn  sin(tn ),   1   O(h 4 )
12
Solutions are: Bounded, periodic
Phase error proportional to
h 2t
Discrete solution shadows the continuous one
Symplectic Methods
The mid-point rule behaves well because it conserves the
symplectic structure of the system
Classical Hamiltonian ordinary differential equation:
dp
H dq

H

 ,

dt 
q dt
p
u  ( p, q)
0I


du

1



J
H
,J




I
0
dt


Differential equation induces a FLOW
 t (u )
FLOW MAP is symplectic
' J'J
T
Symplecticity places a strong constraint on the flows
1.
2.
3.
4.
5.
Preservation of phase space volume (and wedge product)
Recurrence
No evolution on a low dimensional attractor
KAM behaviour for near integrable systems
Composition of two symplectic flows is a symplectic flow
Numerical method applied with a constant step size h gives a
map 
h
Traditional analysis:
Show that
GI approach:
Show that 
h


O
(h)
h
h
p
is a symplectic map (symplectic method)
Advantage:
Symplectic methods have good ergodic properties
Strong error estimates via backward error analysis
method is exact solution of a perturbed Hamiltonian problem
Symplectic Methods include: Runge-Kutta, Splitting, Variational
Runge-Kutta methods for du/dt = f(u)
There is a large class of implicit symplectic Runge-Kutta methods
Butcher Tableaux
c
A
b
Construct matrix M
M
b

b

b
ij
ia
ij
ja
ji
ib
j
Method is symplectic if M = 0
All linear and quadratic invariants
T
u Lu
conserved
Example: the implicit mid-point rule
u

u

hf
((
u

u
)
/
2
)
n

1
n
n
n

1
Symplectic, implicit, symmetric, unconditionally stable,
Conserves linear and quadratic invariants
All Gauss-Legendre Runge-Kutta methods and associated
collocation methods are symplectic
Splitting and composition methods
Runge-Kutta methods are implicit, but for certain problems we
can construct explicit symplectic methods via splitting
du
f1(u
)f2(u
)
dt
Construct flow maps
1,h and 2,h
for
f1
and
Compose the split maps

T,h 
1,h
2,h
Lie-Trotter splitting




S
,h
1
,h
/2
2
.h
1
,h
/2
Strang splitting
f2
Some important results
If 1, h
T ,h S ,h
and 2 , h are symplectic, so are
The Campbell-Baker-Hausdorff theorem implies that
2




O
(
h
)
h
T
,h


O
(h)
h
S
,h
3
If
H(u) = T(p) + V(q)
The splittings lead directly to two important numerical methods
Symplectic Euler SE
p

p

h

V
(
q
),
q

q

h

T
(
p
)
n

1
n
n
n

1
n
n

1
Symplectic, explicit, non-symmetric, order 1
Stormer-Verlet SV (Leapfrog)
h
h
p

p


V
(
q
),
q

q

h

T
(
p
),
p

p


V
(
q
)
n

1
/
2
n
n
n

1
n
n

1
/
2n

1
n

1
/
2
n

1
2
2
Symplectic, explicit, symmetric, order 2
Unstable for large step size
There are higher order, explicit, splitting methods due to Yoshida,
Blanes.
Apply to the Kepler problem
FE
SE
SV
Global error
FE
H error
SE
Method
FE
SE
SV
Global error
H error
L error
t^2 h
th
t h^2
th
h
h^2
th
0
0
NOTE: Kepler’s third law is NOT conserved by these methods …
see the next talk!
Backward Error Analysis
 h* / h
)
Up to an exponentially small (In h) error E  O(e
the solutions of a symplectic method of order p are the
discrete samples of a solution of a related Hamiltonian
differential equation with Hamiltonian
p1
Hh (u)  H (u)  h H p (u)  h H p1 (u)  ...
p
Can construct the perturbed Hamiltonian explicitly
H error remains bounded for all times
Doesn’t apply if h varies!
Example: A problem in structural mechanics
Discrete Euler Beam
Small h limit
p2
2
H
 1 q
2
Hamiltonian for Symplectic Euler discretisation = original problem
2
p
h pq
2
2
H h  H  hH 1  O(h ) 
 1 q 
 O(h 2 )
2
2 1 q2
h = 0.05
Hh
H  hH1
h = 1.1
h = 2.2
Symmetry Group Methods
Spo
Important class of GI methods are used to solve problems
tty
with Lie Group Symmetries (deep conservation laws)
dog
du
 A(t , u )u, u  G,
dt
G: (matrix) Lie group
Eg. G = SO3 (rotations),
A g
g: Lie algebra
g = so3 (skew symmetry)
Can a numerical method ensure that the solution remains in G?
Rigid body mechanics, weather forecasting, quantum mechanics, Lyapunov
exponents, QR factorisation
Idea: Do all computations in the Lie Algebra (linear space)
And map between this and the Lie Group (nonlinear space)
G
Un
n
numerical method
g
U n1
 n1
Examples of maps from g to G
An
  exp(A)   ,
n!
General g , G
  cay( A)  ( I  A / 2) 1 ( I  A / 2)
g = so3, G = SO3
satisfies the dexpinv equation
Bl l
d
  ad  A(t , e U n )
dt
l!
Integrate the dexpinv equation numerically
Conserve the group structure by making sure that all
numerical approximations to the dexpinv equation always
lie in the Lie algebra
Fine provided method uses linear operations and commutators
Runge-Kutta/Munthe-Kaas (RKMK) methods use this approach
Magnus series methods:
du
 A(t )u
dt
Magnus series:
1 t 1
 (t )   A( )d    [ A( 2 ), A(1 )]d1d 2
0
2 0 0
t
1 t 1  2
    [[ A( 3 ), A( 2 )], A(1 )]d1d 2 d 3  ...
4 0 0 0
Obtain method by series truncation and careful calculation of the
commutators
VERY effective for Highly Oscillatory Problems [Iserles]
Eg. Evolution on the surface of the sphere
0
t
 cos(t ) / 2 



A(t )    t
0
t / 2   so3
 cos(t ) / 2  t / 2

0


I u u
T
invariant
FE
RKMK
RK
Magnus