Adaptivity and symmetry for ODEs and PDEs Chris Budd Basic Philosophy ….. • ODES and PDEs develop structures on many time and length.

Download Report

Transcript Adaptivity and symmetry for ODEs and PDEs Chris Budd Basic Philosophy ….. • ODES and PDEs develop structures on many time and length.

Adaptivity and symmetry for ODEs and PDEs
Chris Budd
Basic Philosophy …..
• ODES and PDEs develop structures
on many time and length scales
• Structures may be uncoupled (eg. Gravity waves and slow
weather evolution) and need multi-scale methods
• Or they may be coupled, typically through (scaling)
symmetries and can be resolved using adaptive methods
Talk will look at
• variable step size adaptive methods for ODES
• scale invariant adaptive methods for PDES
The need for adaptivity: 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
Kepler orbits
Forward
Euler
Symplectic
Euler
Stormer
Verlet
FE
Global error
SV
H error
Larger error at close approaches
t
Kepler’s third law is not respected
Main error
Adaptive time steps are highly desirable for accuracy
and symmetry
But …
Adaptivity can destroy the symplectic shadowing structure
[Calvo+Sanz-Serna]
Adaptive methods may not be efficient as a splitting method
AIM: To construct efficient, adaptive, symplectic
methods EASY which respect symmetries
H error
t
Hamiltonian ODE system:
dq
 H p,
dt
dp
 H q
dt
The Sundman transform introduces a continuous
adaptive time step.
IDEA: Introduce a fictive computational time

dt
 g ( p, q )
d
g ( p, q)
SMALL if solution requires small time-steps
Rescaled system for p,q and t
dq
 g ( p, q ) H p ,
d
dp
  g ( p, q ) H q
d
dt
 g ( p, q )
d
Can make Hamiltonian via the Poincare Transform
New variables
Hamiltonian
pt  t, qt  H0  H ( p(0), q(0))
K ( p, pt , q, qt )  g ( p, q)(H ( p, q)  qt )
Now solve using a Symplectric ODE solver
Choice of the scaling function g(q)
Performance of the method is highly dependent on the
choice of the scaling function g.
Approach: insist that the performance of the numerical
method when using the computational variable should be
independent of the scale of the solution and that the
method should respect the symmetries of the ODE
du
i
The differential equation system
 f i (u1 ,...,u N )
dt
Is invariant under scaling if it is unchanged by the symmetry
i
t  t , u   u,   0
eg. Kepler’s third law relating planetary orbits
It generically admits particular self-similar solutions
satisfying
ui (t )  i ui (t )
Theorem [B, Leimkuhler,Piggott] If the scaling function
satisfies the functional equation
1
N
g ( u1,..., uN )  g (u1,...,uN )
Then
Two different solutions of the original ODE mapped
onto each other by the scaling transformation are the
same solution of the rescaled system scale invariant
A discretisation of the rescaled system admits a
discrete self-similar solution which uniformly
approximates the true self-similar solution for all time
Example: Kepler problem in radial coordinates
A planet moving with angular momentum
with

radial coordinate r = q and with dr/dt = p satisfies a
Hamiltonian ODE with Hamiltonian
p 1 
H
  2
2 q q
2
If
 0
symmetry
t  t , q  2 / 3q,
p  1/ 3 p
Numerical scheme is scale-invariant if
g ( q)  g (q)  g (q)  q
2/ 3
3/ 2
If 0    1 there are periodic solutions with close
approaches
Hard to integrate with a non-adaptive scheme
q
t
Consider calculating them using the scaling
g (q)  q 




0
1
 3/ 2
2
No scaling
Levi-Civita scaling
Scale-invariant
Constant angle change
H Error
Surprisingly
sharp!!!

 opt
3 1
 
2 2n
Method order
Scale invariant methods for PDES
These methods extend naturally to PDES with
scaling and other symmetries
ut  f (u,  xu,  xu,...),
t  t , ( x, y )   ( x, y ), u   u
Examples
ut  u xx  u 3 ,
ut  u xxxx  u 3 ,
Parabolic blow-up
High-order blow-up
iut  u  u u  0,
2
ut  u  .(uv), vt  v  u  v,
NLS
Chemotaxis
ut  (u mu x ) x ,
ut  (qum ) x  v, u (0, t )  0, u ( s, t )  vt.
PME
Rainfall
Need to continuously adapt in time and space
Introduce spatial analogue of the fictive time
Adapt spatially by mapping a uniform mesh from a
computational domain  C into a physical domain  P
C ( , )
 P ( x, y)
Use a strategy for computing the mesh which takes
symmetries into account
Q( , , t )
Introduce a mesh potential
( x(t ), y(t ),..)   Q  (Q , Q ,..)
Geometric scaling
 Q
 ( x, y )
H (Q) 
 det
( , )
 Q
Q 

Q 

Q

Q Q  Q
1D
Control scaling via a measure
2
2D
M (u, ux , u y ,..)
Evolve mesh by solving a MK based PDE
 I    Qt   M (Q) H (Q) 
1/ d
(PMA)
Spatial smoothing
(Invert operator
using a spectral
method)
Averaged
measure
Ensures right-handside scales like P in
d-dimensions to give
global existence
Parabolic Monge-Ampere equation PMA
Because PMA is based on a geometric approach,
it has natural symmetries
1. System is invariant under translations and rotations
2. For appropriate choices of M the system is invariant
under scaling symmetries
t  Tt, ( x, y)  L( x, y), u  Uu
( x, y)  L( x, y)  Q  LQ
L
Qt  Qt ,
T
( M ( LQ) H ( LQ))1/ d  L( M ( LQ) H (Q))1/ d
PMA is scale invariant provided that
M ( LQ)1/ d  M (Uu( L( x, y),Tt ))1/ d  T 1M (u( x, y, t ))1/ d
Example: Parabolic blow-up in d dimensions
ut  uxx  u yy  u3
u   t  t , ( x, y)  ( x , y )
*
*
(t  t )   (t  t ), ( x, y)   ( x, y), u  
1/ 2
1/ 2
Scale:
T  (t  t ), U  T
1 / 2
, L T
1/ 2
log(T )
u
1/ 2
M (T 1/ 2u)1/ d  T 1M (u)1/ d  M (u)  u 2d

Regularise: M ( X , Y , t )  u ( X , Y , t ) 2 d  u ( X , Y , t ) 2 d d d X
Basic approach
• Discretise PDE and PMA in the computational domain
• Solve the coupled mesh and PDE system either
(i) As one large system (stiff!)
or
(ii) By alternating between PDE and mesh
Method admits exact discrete self-similar solutions
ut   X u  u3 solve PMA simultaneously with the PDE
10
10^5
Solution:
Y
Mesh:
X
Solution in the computational domain
10^5
2
1
Same approach works well for the Chemotaxis eqns,
Nonlinear Schrodinger eqn, Higher order PDEs
Now extending it to CFD problems: Eady, Bousinessq