Summerschool in Aveiro (Sept. 2018), Ernst Hairer Part I. Geometric - PowerPoint PPT Presentation

Summerschool in Aveiro (Sept. 2018), Ernst Hairer Part I. Geometric numerical integration ◮ Hamiltonian systems, symplectic mappings, geometric integrators, St¨ ormer–Verlet, composition and splitting, variational integrator ◮ Backward error analysis, modified Hamiltonian, long-time energy conservation, application to charged particle dynamics Part II. Differential equations with multiple time-scales ◮ Highly oscillatory problems, Fermi–Pasta–Ulam-type problems, trigonometric integrators, adiabatic invariants ◮ Modulated Fourier expansion, near-preservation of energy and of adiabatic invariants, application to wave equations Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 1 / 29

Lecture 3. Highly oscillatory systems Numerical energy conservation 1 Molecular dynamics model Oscillatory energy – FPU-type problem Conservation of oscillatory energy (exact solution) Numerical integrators 2 Step size restriction for St¨ ormer–Verlet Exponential integrators IMEX (implicit–explicit) integrators Conservation of energies (numerical solution) Several high frequencies – resonances 3 FPU type problem Example with resonant frequencies Numerical energy conservation Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 2 / 29

Example 1. Molecular dynamics model Consider the N -body problem i − 1 � � N � � q = −∇ U ( q ) , ¨ U ( q ) = V ij � q i − q j � , i =2 j =1 1 V ( r ) = r − 12 − 2 r − 6 where V ij ( r ) = V ( r ) 0 Lennard –Jones −1 1 2 Initial values: positions q : randomly perturbed values on a grid (see the figure) momenta p = ˙ q : zero Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 3 / 29

Numerical simulation St¨ ormer–Verlet method p n − h p n +1 / 2 = 2 ∇ U ( q n ) q n +1 = q n + h p n +1 / 2 p n +1 / 2 − h p n +1 = 2 ∇ U ( q n +1 ) with two different step sizes h = 0 . 08 t = h = 0 . 04 0.00 Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 4 / 29

Let us look at the Hamiltonian We plot the total energy q ) = 1 q T ˙ H ( q , ˙ 2 ˙ q + U ( q ) as a function of time for step sizes: h = 0 . 08 (black), h = 0 . 04 , 0 . 02 , 0 . 01 , 0 . 005 (blue) 10 −1 10 −2 10 −3 10 −4 10 −5 0 20 40 1000 1020 1040 Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 5 / 29

Is the step size small (compared to the frequency)? q = − ω 2 q the St¨ For the harmonic oscillator ¨ ormer–Verlet scheme leads to the recursion � 2( h ω ) 2 � 1 − 1 q n +1 − 2 q n + q n − 1 = 0 , which is stable for 0 < h ω < 2. For the N -body problem ¨ q = −∇ U ( q ) we linearise the vector field, we approximate ω 2 by the largest eigenvalue of the Hessian ∇ 2 U ( q ), and we approximate ω by the square root of the largest eigenvalues. 25 square root of the largest eigenvalues of ∇ 2 U ( q ) 20 46.2 46.3 46.4 The largest ω is between 20 and 25. Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 6 / 29

In the situation of our experiment we have ω max ≈ 25 . Therefore step sizes h = 0 . 08 , 0 . 04 , 0 . 02 , 0 . 01 , 0 . 005 correspond to h ω max ≈ 2 , 1 , 0 . 5 , 0 . 25 , 0 . 125 Comments divergence for h = 0 . 08, because h ω ≈ 2 is too close to the linear stability limit backward error analysis cannot be applied for the considered step sizes, because h ω is not small Conclusion it is an open problem to explain the long-time behaviour of the previous experiment Let us turn to simpler problems which still have high oscillations. Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 7 / 29

Example 2. Chain of alternating soft and stiff springs Q 2 m − 1 Q 2 m Q 1 Q 2 · · · stiff soft 2 m m m � � � Q ) = 1 k + 1 j ( Q 2 j − Q 2 j − 1 ) 2 + H ( Q , ˙ ˙ Q 2 ω 2 ( Q 2 j +1 − Q 2 j ) 4 2 4 k =1 j =1 j =0 √ q j = ( Q 2 j − Q 2 j − 1 ) / 2 compression/extension of a stiff spring √ q j + m = ( Q 2 j + Q 2 j − 1 ) / 2 its mean position q j + ω 2 ¨ j q j = −∇ j U ( q ) j = 1 , . . . , m q j + m = −∇ ¨ j + m U ( q ) Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 8 / 29

Oscillatory energies Oscillatory energies of the stiff springs q j + ω 2 ¨ j q j = −∇ j U ( q ) � � q j ) = 1 q j + m = −∇ ¨ j + m U ( q ) q 2 j + ω 2 j q 2 I j ( q j , ˙ ˙ j 2 Total oscillatory energy I ( q , ˙ q ) = I 1 ( q 1 , ˙ q 1 ) + . . . + I m ( q m , ˙ q m ) I 1 I 1 I 3 I 2 0 0 100 200 Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 9 / 29

Highly oscillatory differential equations Oscillatory Hamiltonian system j = 0 , 1 q j + ω 2 ¨ j q j = − ∇ j U ( q ) � q 1 | 2 � q ) = 1 + 1 Total energy q 0 | 2 + | ˙ 2 ω 2 | q 1 | 2 + U ( q ) H ( q , ˙ | ˙ 2 � q 1 | 2 + ω 2 | q 1 | 2 � q 1 ) = 1 Oscillatory energy I ( q 1 , ˙ | ˙ 2 ω 0 = 0 and ω 1 = ω ≥ 1 arbitrarily large there exist δ > 0 and a set K such that the potential U has bounded derivatives of all orders in a δ -neighborhood of K × { 0 } . the initial values satisfy q 1 (0) = O ( ω − 1 ) , q 0 (0) = O (1) , q 0 (0) = O (1) , ˙ q 1 (0) = O (1) , ˙ � � such that I q 1 (0) , ˙ q 1 (0) ≤ E , where the bound E is independent of ω . Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 10 / 29

Preservation of the oscillatory energy � � � � We know that H q ( t ) , ˙ q ( t ) − H q (0) , ˙ q (0) = 0 for all t . Theorem (Benettin, Galgani & Giorgilli, 1987) For an arbitrarily fixed integer N ≥ 1 there exist C > 0 and ω ∗ > 0 such that the following holds for ω ≥ ω ∗ : along the exact solution the total oscillatory energy deviates from its starting value by no more than � � � � �� � ≤ C ω − 1 � I 0 ≤ t ≤ ω N , q 1 ( t ) , ˙ q 1 ( t ) − I q 1 (0) , ˙ q 1 (0) for provided that q 0 ( t ) stays in the set K for such long times. The threshold ω ∗ and the constant C depend on N, on the energy bound E and on bounds of derivatives of the potential U. This result explains the excellent long-time preservation of the total oscillatory energy in Example 2 (for the exact solution). Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 11 / 29

Lecture 3. Highly oscillatory systems Numerical energy conservation 1 Molecular dynamics model Oscillatory energy – FPU-type problem Conservation of oscillatory energy (exact solution) Numerical integrators 2 Step size restriction for St¨ ormer–Verlet Exponential integrators IMEX (implicit–explicit) integrators Conservation of energies (numerical solution) Several high frequencies – resonances 3 FPU type problem Example with resonant frequencies Numerical energy conservation Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 12 / 29

Step size restriction for St¨ ormer–Verlet Harmonic oscillator q + ω 2 q = 0 ¨ St¨ ormer–Verlet q n +1 − 2 q n + q n − 1 +( h ω ) 2 q n = 0 The characteristic equation of this difference relation is � � 1 − ( h ω ) 2 ζ 2 − 2 ζ + 1 = 0 2 and we have stability (bounded solutions) if and only if � � � 1 − ( h ω ) 2 � � � < 1 i.e. 0 < h ω < 2 2 Aim. Consider numerical integrators without such a severe step size restriction. Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 13 / 29

Exponential integrators q + Ω 2 q = −∇ U ( q ) we consider For a system ¨ q n +1 − 2 cos( h Ω) q n + q n − 1 = − h 2 Ψ ∇ U (Φ q n ) 2 h sinc ( h Ω) p n = q n +1 − q n − 1 where Ψ = ψ ( h Ω) , Φ = φ ( h Ω) and ψ ( x ), φ ( x ) are even functions satisfying ψ (0) = 1, φ (0) = 1. Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 14 / 29

Exponential integrators q + Ω 2 q = −∇ U ( q ) we consider For a system ¨ q n +1 − 2 cos( h Ω) q n + q n − 1 = − h 2 Ψ ∇ U (Φ q n ) 2 h sinc ( h Ω) p n = q n +1 − q n − 1 where Ψ = ψ ( h Ω) , Φ = φ ( h Ω) and ψ ( x ), φ ( x ) are even functions satisfying ψ (0) = 1, φ (0) = 1. One-step formulation cos( h Ω) q n + h sinc ( h Ω) p n + 1 2 h 2 Ψ g n q n +1 = � � − Ω sin( h Ω) q n + cos( h Ω) p n + 1 p n +1 = 2 h Ψ 0 g n + Ψ 1 g n +1 where g n = −∇ U ( φ q n ), Ψ 0 = ψ 0 ( h Ω), Ψ 1 = ψ 1 ( h Ω), and the functions ψ 0 ( x ) and ψ 1 ( x ) are given by ψ ( x ) = sinc ( x ) ψ 1 ( x ) , ψ 0 ( x ) = cos( x ) ψ 1 ( x ) . Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 14 / 29

Interpretation as splitting method The previous one-step formulation can be split into � p n � � p + � � p − � � p n +1 � n n +1 �→ �→ �→ q n q n q n +1 q n +1 where p n − h p + = 2Ψ 1 ∇ U (Φ q n ) n n +1 − h p − p n +1 = 2Ψ 1 ∇ U (Φ q n +1 ) and the central mapping is the exact flow of the differential equation q + Ω 2 q = 0. ¨ Ernst Hairer (Universit´ e de Gen` eve) Geometric Numerical Integration September 10 -14, 2018 15 / 29