7OyFEb-lecture8

7OyFEb-lecture8 - Symplectic Integration Introduction...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Symplectic Integration Introduction Realistic Objectives for Molecular Dynamics Simulations In general, the aim of a MD simulation is to identify qualitative and statistical properties of molecular motions rather than to reproduce a precise trajectory. Large system size and small time steps lead to compounding of roundoff error. Chaotic dynamics cause nearby trajectories to diverge rapidly. Initial conditions are poorly known, even when experimental structural data is available. Model misspecification: Numerous approximations: Newtonian mechanics, solvent models, periodic BCs. Force field parameters are estimated. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 1 / 104 Symplectic Integration Introduction Chaotic Dynamics Local instabilities in MD simulations lead to rapid but bounded divergence of trajectories. Figures show the RMSD between simulations of a water tetramer. Divergence is initially exponential. Saturation occurs within 10 ps. The timestep has little effect. From Schlick (2006). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 2 / 104 Symplectic Integration Introduction Implications for Numerical Analysis The choice of the numerical integration scheme used in a MD simulation is influenced by several considerations: Simulation of large macromolecules on nanosecond or longer timescales imposes severe demands on time and memory. Stability is at least as important as order of accuracy. Geometric invariants of the system dynamics can sometimes be exploited. Flexible implementation for use with constrained or stochastic dynamics, or with multiple-timestep methods. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 3 / 104 Symplectic Integration Introduction The Strmer-Verlet Algorithm o One of the most popular integrators used in MD is the Strmer-Verlet o algorithm (position Verlet): Xn+1/2 = Xn + Vn+1 Xn+1 t n V 2 = Vn + F(Xn+1/2 ) t n+1 = Xn+1/2 + V 2 The position Verlet algorithm has the following properties: explicit; second-order accurate; symplectic and time-reversible. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 4 / 104 Symplectic Integration Hamiltonian Mechanics Hamiltonian Form of Newtonian Mechanics Recall that Newton's equation of motion can be expressed as a Hamiltonian system: p=- where q is the coordinate vector; p = Mq is the momentum vector; the Hamiltonian H gives the total energy of the system: 1 H(p, q) = pT M-1 p + U(q). 2 H , q q= H p Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 5 / 104 Symplectic Integration Hamiltonian Mechanics Hamiltonian Systems More generally, we say that a system of ODE's on an open domain Rd Rd is an autonomous Hamiltonian system if there is a smooth function H : R such that p=- Here d is the degrees of freedom of the system. q Rd and p Rd are the generalized coordinates and generalized momenta. One can also consider time-dependent Hamiltonians H(p, q, t). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 6 / 104 H , q q= H . p Symplectic Integration Hamiltonian Mechanics The Flow Map Given any Hamiltonian system, for each point (p0 , q0 ) there is a T > 0 such that the initial value problem H , q p(0) = p0 , p=- H p q(0) = q0 q= has a solution t = (t, p0 , q0 ) for all t [0, T ]. If T = for all points (p0 , q0 ) , then the flow map t : is a diffeomorphism of and the family of flow maps (t ; t 0) is a one-parameter semigroup of diffeomorphisms on : t s = t+s . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 7 / 104 Symplectic Integration Hamiltonian Mechanics Conservation of Energy Hamiltonian systems have several important invariants. One of these is the Hamiltonian itself: H(p(t), q(t)) = t = H H p+ q p q H H H H - + = 0. p q q p In mechanics, this corresponds to conservation of the total energy. Dynamically, it means that the trajectories lie on the level sets of H. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 8 / 104 Symplectic Integration Hamiltonian Mechanics Liouville's Theorem Phase space volume is also conserved by Hamiltonian flows. Liouville's Theorem: If D is a measurable subset of , then for all t 0 V (t) dx = dx = V (0). t (D) D Proof: If D has a piecewise smooth boundary D, then the divergence theorem tells us that (t) = V n(y) F(y)dy = F(x)dx D D where F is the Hamiltonian vector field and n(y) is the outward unit normal vector at y D. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 9 / 104 Symplectic Integration Hamiltonian Mechanics However, since F= , p q H H 2H 2H - , =- + =0 q p pq qp we see that the Hamiltonian vector field is divergence-free and v (t) = 0. This implies that V (t) = V (0) for all t 0. To deal with measurable D, observe that t t (D) = dx = x dx t (D) D defines a measure on the Borel -algebra () and that () is generated by the open balls contained in (which have smooth boundaries). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 10 / 104 Symplectic Integration Differentiable Manifolds Topological Manifolds Definition: An n-dimensional topological manifold M is a second-countable Hausdorff space that is locally homeomorphic to Rn , i.e., for every x M there is a pair (U, ) such that U M is an open neighborhood of x; : U Rn is a 1-1 map which is continuous and has a continuous inverse -1 defined on the open set V = (U). (U, ) is said to be a chart on M with coordinate neighborhood U and coordinates (q). Informally, M is a space that locally `looks like' Rn . Any open subset of Rn is an n-dim manifold. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 11 / 104 Symplectic Integration Differentiable Manifolds Coordinate Transformations and Atlases In general, a manifold will be endowed with many overlapping coordinate neighborhoods. An atlas for an n-dim manifold M is a collection of charts (U , ) for which the coordinate neighborhoods form an open cover of M. If U U = , then the coordinate transformations = -1 : (U U ) (U U ) are homeomorphisms. = -1 : (U U ) (U U ) Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 12 / 104 Symplectic Integration Differentiable Manifolds Smooth Atlases We can also insist that the coordinate transformations have additional regularity. Two charts (U , ) and (U , ) are said to be smoothly compatible if the coordinate transformations and are smooth maps. An atlas A is said to be a smooth atlas if every coordinate transformation is a smooth map. A smooth atlas A is said to be maximal if A contains every chart (U, ) that is smoothly compatible with A. Exercise: Show that each smooth atlas A is contained in a unique maximal smooth atlas. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 13 / 104 Symplectic Integration Differentiable Manifolds Differentiable Manifolds Definition: An n-dimensional differentiable manifold is an n-dimensional topological manifold M together with a smooth maximal atlas M. M is said to be a differentiable structure on M. By the exercise, we can specify a differentiable structure M by choosing any smooth atlas A M. Easy example: Any open set U Rn can be made into an n-dimensional differentiable manifold by equipping U with the atlas A = {(U, Id|U )}, where Id|U : U U is the identity map. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 14 / 104 Symplectic Integration Differentiable Manifolds The n-Sphere Example: The n-sphere Sn = {x Rn+1 : |x| = 1} can be equipped with a differentiable structure by using stereographic projection from each of two antipodal points: N : UN S : US = Sn \{N = (1, 0, , 0)} Rn ; = Sn \{S = (-1, 0, , 0)} Rn ; x (x1 , , xn ) 1 - x0 (x1 , , xn ) x . 1 + x0 Indeed, since the coordinate transformations, NS (y) = SN (y) = y , ||y||2 are both smooth functions on N (UN US ) = S (UN US ) = Rn \{0}, it follows that the two charts are smoothly compatible. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 15 / 104 Symplectic Integration Differentiable Manifolds Submanifolds Definition: A subset M of an n dimensional differentiable manifold is an m-dimensional submanifold if for every p M, there is a chart (U, ) for N such that p U and (U M) = (U) (Rm {0}) Rn . In this case, (U M, |UM ) is a chart for M and, endowed with this differentiable structure, M is an m-dimensional differentiable manifold. Whitney Embedding Theorem: Every n-dimensional differentiable manifold is diffeomorphic to an n-dimensional differentiable submanifold of R2n . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 16 / 104 Symplectic Integration Differentiable Manifolds Differentiable Maps A map F : M N between m- and n-dimensional manifolds M and N is said to be differentiable (smooth) at p M if F is continuous at p and if there are charts (U, ) at p and (V , ) at F (p) such that F -1 : (U F -1 (V )) Rm Rn is differentiable (smooth) at (p). F is said to be differentiable on an open subset U M if F is differentiable at every p U. F is said to be a diffeomorphism between M and N if F is a homeomorphism and both F and F -1 are differentiable. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 17 / 104 Symplectic Integration Differentiable Manifolds Claim: Differentiability is well-defined. Proof: Suppose that (U1 , 1 ) and (U2 , 2 ) are charts at p M, that (V1 , 1 ) and (V2 , 2 ) are charts at f (p) N, and that 1 F -1 is 1 smooth at p. Since each pair of charts is compatible, the coordinate -1 transformations 21 = 1 -1 and 12 = 2 1 are local 2 diffeomorphisms and consequently 2 F -1 = 12 1 F -1 21 2 1 is smooth at p. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 18 / 104 Symplectic Integration Differentiable Manifolds Regular and Critical Values Definition: Suppose that M and N are m- and n-dimensional differentiable manifolds and that f : M N is differentiable. The rank of f at p M is the rank of the Jacobian map D( f -1 )((p)) : Rm Rn , where (U, ) and (V , ) are charts at p and f (p), respectively. p M is said to be a critical point for f and f (p) is called a critical value if the rank of f at p is less than n. y N is said to be a regular value for f if f has rank n at every point in f -1 (y ) M. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 19 / 104 Symplectic Integration Differentiable Manifolds Submanifolds as Fibers of Smooth Maps Proposition: If M and N are m- and n-dimensional differentiable manifolds and y N is a regular value of a smooth map f : M N, then f -1 (y ) is an (m - n)-dimensional submanifold of M. Example: If f : Rn+1 R is defined by f (x) = |x|2 , then f is a smooth map with derivative f (x) = 2x, which has rank 1 at every x = 0. It follows that S n = f -1 (1) is an n-dimensional submanifold of Rn+1 . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 20 / 104 Symplectic Integration Differentiable Manifolds Tangent Vectors and Tangent Spaces Let M be an n-dim manifold and let p (M) be the collection of smooth curves : R M with (0) = p M. We can define an equivalence relation on p (M) by stipulating that 1 2 if for any a chart (U, ) at p we have ( 1 ) (0) = ( 2 ) (0). A tangent vector to M at p is an equivalence class [1 ]p = {2 p (M) : 2 1 }. The set of such tangent vectors is denoted Tp (M) and is called the tangent space to M at p. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 21 / 104 Symplectic Integration Differentiable Manifolds Vector Operations in Tangent Spaces Tp (M) has a natural vector space structure that makes it isomorphic to Rn . If (U, ) is a chart at p with (p) = 0, then there is a one-to-one correspondence between Tp (M) and Rn defined by (p ) = ( ) (0). We can then use -1 to define vector operations in Tp (M): [1 ]p + [2 ]p -1 ([1 ]p ) + ([2 ]p ) (, R). Exercise: Show that the vector space structure on Tp (M) does not depend on the choice of (U, ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 22 / 104 Symplectic Integration Differentiable Manifolds Local Coordinates The coordinate function can also be used to assign a basis to Tp (M). Suppose that e1 , , en is a basis for Rn and that xi : (-1, 1) Rn is defined by xi (t) = tei . Then {x1 , , xn } is a basis for Tp (M) where xi = [-1 xi ]p . Caveat: Unlike the vector space structure, this basis does depend on the chart. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 23 / 104 Symplectic Integration Differentiable Manifolds The Tangent Bundle Definition: If M is an n-dimensional differentiable manifold, then the tangent bundle of M is the 2n-dimensional manifold TM = pM Tp (M) where: For each chart (U, ) on M, we define a chart (T (U), ) on TM with T (U) = pU Tp M and (v) = ((p), (v)), v Tp (M). In particular, every point v TM has a neighborhood T (U) that is diffeomorphic to U Rn . Examples: T Rn R2n ; TS 1 S 1 R1 . = = Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 24 / 104 Symplectic Integration Differentiable Manifolds Vector Fields on Manifolds Definition: A smooth vector field on a differentiable manifold M is a smooth mapping X : M TM with the property that X (p) Tp (M) for every p M. Example: Suppose that A is the matrix 0 -1 A= . 1 0 If we regard S1 as the embedded unit circle in R2 , then X : S1 T S1 T R2 : X (v) = (Av)v is a non-vanishing smooth vector field on S1 . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 25 / 104 Symplectic Integration Differentiable Manifolds Derivatives Suppose that f : M N is a map between manifolds. If f is differentiable at p M, then the derivative of f at p is the map f (p) : Tp M Tf (p) N defined by f (p)(p ) = [f ]f (p) . If f is differentiable on all of M, then the derivative of f is the map f : TM TN defined by f (p ) = f (p)(p ). Example: If M and N are open subsets of Rm and Rn , then f (p) = Df (p), where Df (p) is the Jacobian of f . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 26 / 104 Symplectic Integration Differentiable Manifolds Reminder: Duality in Linear Algebra If V is an n-dimensional (real) vector space, then the (algebraic) dual is the n-dimensional vector space V containing the linear functionals w : V R. If A : V W is a linear map, then the dual map A : W V is defined by A (w ) = w A V Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 27 / 104 Symplectic Integration Differentiable Manifolds Cotangent Spaces Definition: A cotangent vector at p M is a linear functional on Tp (M). The collection of such vectors is a vector space Tp (M) called the cotangent space to M at p. If (x1 , , xn ) are local coordinates on M in a neighborhood of p, then the dual basis for Tp (M) is given by (dx 1 , , dx n ), where dx i (xj ) = ij . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 28 / 104 Symplectic Integration Differentiable Manifolds Differentials Example: If f : M R is smooth, then f (p) : Tp (M) Tf (p) (R) R = is a linear functional on Tp (M) called the differential of f at p and denotes df (p). Given local coordinates, we can write df (p) = n f (p) dx i xi i=1 which shows that df (p)(v) is the directional derivative of f at p in the direction v Tp (M). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 29 / 104 Symplectic Integration Differentiable Manifolds The Cotangent Bundle As with the tangent spaces, the union of the cotangent spaces on an n-dimensional manifold is a 2n-dimensional manifold called the cotangent bundle: T M = pM Tp (M). If (U, ) is a chart at p, then we take (T (U), ) to be a chart for T M at v Tp (M) with T (U) = pU Tp (M) and (v ) = ((p), (-1 ) v ) U Rn . where we identify the vector spaces T(p) (Rn ) Rn . = Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 30 / 104 Symplectic Integration Exterior Algebra Alternating Forms Definition: An alternating form of degree k (k-form for short) on a vector space V is an antisymmetric multilinear functional : : V V R Multilinearity means that is linear in each coordinate: (v1 , , vi + wi , , vn ) = (v1 , , vi , , vn ) + (v1 , , wi , , vn ) Antisymmetry means that (v(1) , , v(n) ) = (-1)|| (v1 , , vn ) where is a permutation of {1, , n} with sign ||. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 31 / 104 Symplectic Integration Exterior Algebra Examples: A 1-form is just an element of V . If S is a skew-symmetric n n matrix, i.e., S T = -S, then T T T (v1 , v2 ) = v1 Sv2 = v2 S T v1 = -v2 Sv1 = -(v2 , v1 ) is a 2-form. If e1 , , en is a basis for V = Rn and P is the orthogonal projection onto the span of ei1 , , eik , then (v1 , , vk ) = det(Pv1 , , Pvk ) is a k-form. The number (v1 , , vk ) is just the oriented volume of the k-paralleliped formed by the vectors Pv1 , , Pvk . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 32 / 104 Symplectic Integration Exterior Algebra Alternating Forms and Linear Independence Proposition: If is an alternating k-form on V and v1 , , vk is a collection of linearly dependent vectors, then (v1 , , vk ) = 0 Proof: We first observe that the antisymmetry of implies that (v1 , , w, , w, , vk ) = -(v1 , , w, , w, , vk ) = 0, i.e., annihilates any k-tuple that contains a repeated vector. Now, w.l.o.g., suppose that we can write vk = k-1 ai vi . Then i=1 (v1 , , vk ) = as claimed. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 33 / 104 k-1 i=1 ai (v1 , , vk-1 , vi ) = 0, Symplectic Integration Exterior Algebra Example: Alternations If : V V R is a multilinear map, then the alternation of is the k-form defined by Alt()(v1 , , vk ) = where Sk is the symmetric group of permutations of {1, , k}; () {1} is the sign of the permutation. 1 ()(v(1) , , v(k) ), k! Sk Since ( ij ) = -() for any transposition ij , it follows that Alt() is antisymmetric. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 34 / 104 Symplectic Integration Exterior Algebra Alternation is a Projection If is a k-form, then Alt()(v1 , , vk ) = = 1 ()(v(1) , , v(k) ) k! 1 ()2 (v1 , , vk ) k! Sk Sk = (v1 , , vk ), and so Alt() = . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 35 / 104 Symplectic Integration Exterior Algebra Exterior Powers The set of k-forms over V is a vector space denoted k (V ) under pointwise addition and scalar multiplication: (1 + 2 )(v) = 1 (v) + 2 (v) ()(v) = (v). We can then form the graded vector space (V ) = n k=0 k (V ) where by convention we take the space of 0-forms to be 0 = R. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 36 / 104 Symplectic Integration Exterior Algebra The Exterior Product We can also define a multiplication on (V ). Definition: The wedge product (the exterior product) of two forms k (V ) and l (V ) is the k + l form k+l (V ) defined as ( )(v1 , , vk+l ) = where the tensor product is the multilinear functional ( )(v1 , , vk+l ) = (v1 , , vk )(vk+1 , , vk+l ). k +l Alt( ) k Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 37 / 104 Symplectic Integration Exterior Algebra Alternative Expression for the Exterior Product The exterior product can also be expressed more parsimoniously as ( )(v1 , , vk+l ) = ()(v(1) , , v(k) )(v(k+1) , , v(k+l) ) Sh(k,l) where Sh(k, l) Sk+l is the set of (k, l)-shuffles: Sh(k, l) if (1) < (2) < < (k); (k + 1) < (k + 2) < < (k + l). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 38 / 104 Symplectic Integration Exterior Algebra Example: Exterior Products of 1-Forms Suppose that , are 1-forms over V . Then their exterior product is the determinantal 2-form (v1 ) (v2 ) . ( )(v1 , v2 ) = (v1 )(v2 ) - (v2 )(v1 ) = (v1 ) (v2 ) In particular, = -( ) for any , 1 (V ). = 0 for any 1 (V ). More generally, if 1 , , k are 1-forms, then their exterior product is the determinantal k-form (1 k )(v1 , , vk ) = det i (vj ) . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 39 / 104 Symplectic Integration Exterior Algebra Basis Sets Suppose that (e1 , , en ) is a basis for V and that (e , , e ) is the n 1 dual basis for V . Then the collection of k-fold exterior products ei ei , 1 k 1 i1 < < ik n is a basis for k (V ). Indeed, any k-form is uniquely determined by its values i1 , ,ik on the k-tuples (ei1 , , eik ) and so = It follows that The k'th exterior power k (V ) has dimension The exterior algebra Jay Taylor (ASU) i1 <<ik i1 , ,ik ei ei . 1 k n k (V ) has dimension 2n = ; n n k=0 k . 40 / 104 APM 530 - Lecture 8 Fall 2010 Symplectic Integration Exterior Algebra Interpretation The basis forms ei ei have the following interpretation. If 1 k (v1 , , vk ) is a k-tuple of vectors, then is the oriented volume of the projection of the k-parallelepiped given by (v1 , , vk ) onto the k-dimensional plane with coordinates (ei1 , , eik ). In particular, if n (V ) is an n-form, then = a e1 en (ei ei )(v1 , , vk ) = det ei (vj ) = det vj,ik 1 k k for some scalar a and so (v1 , , vn ) = a det(vj ) is just a scalar multiple of the volume of the parallelepiped formed by the n vectors. is called a volume form. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 41 / 104 Symplectic Integration Differential Forms Differential Forms If M is an n-dimensional differentiable manifold, then we write k (M) = k (Tp (M)) p and k (M) = pM k (M) p for the k'th exterior power of the cotangent space at p M and the union n over M of these spaces. k (M) is an n k -dimensional manifold called the bundle of k-forms on M: Each chart (U, ) on M induces a chart (k (U), ) on k (M) where n : k (U) = pU k (M) k ((U)) (U) R(k ) . = p Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 42 / 104 Symplectic Integration Differential Forms Differential Forms A smooth map : M k (M) is called a differential k-form on M. If x = (x1 , , xn ) are local coordinates in some open U M, then can be expressed as = i1 , ,ik (x)(dxi1 dxik ). i1 <<ik where each function i1 , ,ik : U R is smooth. We write k (M) for the space of all k-forms on M. If 1 , 2 k (M), then 1 + 2 k (M). If is a k-form on M and f : M R is a smooth map, then f is also a k-form on M. k (M) is a module over the ring of smooth real-valued functions on M. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 43 / 104 Symplectic Integration Differential Forms Examples Every smooth function f : M R is a 0-form. n f dxi . xi i=1 Given such f , df is a differential 1-form on M: df = where (x1 , , xn ) are local coordinates. Given k (M) and l (M), we can define the exterior product k+l (M) by pointwise multiplication: ( )p = p p . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 44 / 104 Symplectic Integration Differential Forms The Pull-back of a Differential Form Definition: If f : M N is a smooth map at p M, then the derivative f (p) induces an algebra homomorphism f from k(p) (N) to k (M) p f defined by f ()(p) v1 , , vk = (f (p)) f (v1 ), , f (vk ) . The k-form f is called the pull-back of by f . To say that f is an algebra homomorphism means that f (1 1 + 2 2 ) = 1 f (1 ) + 2 f (2 ); f (1 2 ) = f (1 ) f (2 ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 45 / 104 Symplectic Integration Differential Forms The Exterior Derivative on Rn We first consider differential forms defined on an open subset U Rn . In the following, we will use I to denote multi-indices of length k and write dx I = dx i1 dx ik . Definition: If = I aI dx I is a k-form on U, then we define the exterior derivative of to be the (k + 1)-form d = daI dx I I where daI 1 (U) is the differential of aI . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 46 / 104 Symplectic Integration Differential Forms The Exterior Derivative This definition induces a sequence of maps - k-1 (U) - k (U) - k+1 (U) - with the following properties: Linearity: d is R-linear on k (U): d(1 + 2 ) = d1 + d2 ; Product rule: If 1 is a k-form, then d(1 2 ) = d1 2 + (-1)k 1 d2 ; Cochain property: d 2 = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 47 / 104 d d d d Symplectic Integration Differential Forms Proof of the Product Rule If f , g smooth functions on U, then d(fg ) = (df )g + f (dg ). Thus, if are 1 = I aI dx I k (U) and 2 = J bJ dx J l (U), then d(1 2 ) = I ,J I ,J I ,J d(aI bJ ) dx I dx J = = bJ daI dx I dx J + I J daI dx bJ dx + (-1)k k I ,J aI dbJ dx I dx J I ,J aI dx I dbJ dx J = d1 2 + (-1) 1 d2 , where the (-1)k comes from the interchange of the 1-form dbJ with the k-form dx I . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 48 / 104 Symplectic Integration Differential Forms Proof of the Cochain Property First note that d(dx i ) = 0 since da = 0 when a is a constant function on U. Consequently, if a is a smooth function on U, then a a d(da) = d dx j = d dx j xj xj j j 2a = dx i dx j xi xj i,j = 0 since dx i dx j = -dx j dx i when i = j and dx i dx i = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 49 / 104 Symplectic Integration Differential Forms Next, let = aI dx I be a k-form on U. Then d(d) = d daI dx I I = = 0 I I d(daI ) dx I - daI d(dx I ) since d(daI ) = d(dx I ) = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 50 / 104 Symplectic Integration Differential Forms The Exterior Derivative on Manifolds Definition: If is a k-form on a manifold M, then we define d to be the (k + 1) form such that for each chart (U, ) on M |U = d(-1 ) () where : k ((U)) k (U) is the isomorphism of exterior powers = induced by the isomorphism : Tp (M) T(p) ((U)). = Remarks: It can be shown that is well-defined, i.e., if (V , ) is a second chart and U V = , then d(-1 ) () = d( -1 ) () . The properties of d shown for open subsets of Rn are also true on manifolds. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 51 / 104 Symplectic Integration Differential Forms Integration on Manifolds Integrals of differentiable forms are defined using certain geometric objects called cells or singular polyhedra rather than through measures. Definition: A k-cell in an n-dimensional differentiable manifold M is a triple = (D, f , O) where D Rk is a convex polyhedron; f : D M is a smooth map; O is an orientation on Rk . Examples: A 1-cell is an oriented curve in M. A 2-cell is an oriented surface in M. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 52 / 104 Symplectic Integration Differential Forms Integration of Differential Forms Definition: The integral of a k-form over a k-cell contained in M is defined as = f = a(x)dx1 dxk , D D where f = a(x)dx 1 dx k is the pullback of under f . Remark: This is an oriented integral. In particular, if we let - denote the k-cell with the reverse orientation, then = - . - Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 53 / 104 Symplectic Integration Symplectic structures Symplectic Structures on Manifolds Definition: A symplectic structure on a 2n-dimensional manifold M is a differential 2-form satisfying: is closed: d = 0; is non-degenerate: for every p M and v Tp (M) with v = 0 there exists w Tp (M) such that (v, w) = 0. The pair (M, ) is called a symplectic manifold and is called a symplectic form. Remark: There can be topological obstructions to the existence of symplectic structures, e.g., S2 is the only sphere that admits a symplectic structure. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 54 / 104 Symplectic Integration Symplectic structures Symplectic Structures on Euclidean Phase Spaces Example: Suppose that U is an open subset of Rn and let (q, p) be coordinates for the 2n-dimensional cotangent bundle M = T U U Rn . = Then = dp dq = is a symplectic form on T U. n i=1 dpi dqi = n i=1 d(pi dqi ) Proof: That is closed follows from the calculation d = n i=1 d 2 (pi dqi ) = 0 and the cochain property d 2 = 0 that holds for any differential form . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 55 / 104 Symplectic Integration Symplectic structures To show that is non-degenerate, let (q1 , , qn , p1 , , pn ) be a basis for Tq,p (M), where dqi (qj ) = ij , dpi (pj ) = ij , Then (dpl dql )(pi , qj ) = dpl (pi )dql (qj ) - dpl (qi )dql (pj ) = li lj (dpl dql )(qj , pi ) = -(dpl dql )(pi , qj ) = -li lj . dqi (pj ) = 0 dpi (qj ) = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 56 / 104 Symplectic Integration Symplectic structures These computations show that the matrix representation of with respect to the basis (q1 , , qn , p1 , , pn ) is (v, w) = (Jv)T w = vT JT w where J is the 2n 2n structure matrix 0 In J= . -In 0 In particular, given any v = 0 and w = Jv, we have (v, w) = vT v = ||v||2 > 0, which shows that is non-degenerate. Remark: We will refer to this symplectic structure = dp dq as the canonical symplectic structure on the phase space T U. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 57 / 104 Symplectic Integration Symplectic structures Symplectic Maps Definition: Let (M1 , 1 ) and (M2 , 2 ) be symplectic manifolds of dimension 2n. A smooth map : M1 M2 is said to be symplectic if 2 = 1 , i.e., at every p M1 1 (p)(v, w) = ( 2 )(p)(v, w) = 2 ((p)) v, w , Remarks: A symplectomorphism is a diffemorphism that is also symplectic. Darboux's theorem shows that every 2n-dimensional symplectic manifold is locally symplectomorphic to R2n endowed with the canonical symplectic form. v, w Tp (M). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 58 / 104 Symplectic Integration Symplectic structures Surface Area Invariants of Symplectic Maps Suppose that : (M, ) (M, ) is a symplectic map and that : D M is a 2-chain in M. Then = is also a 2-chain in M and ~ = = = ~ = ~ D D D . Interpretation: A symplectic map preserves the oriented area embedded surface . of any Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 59 / 104 Symplectic Integration Symplectic structures The Surface Area Invariant of the Canonical Form Suppose that : D M is a 2-chain in M = T U. If i denotes the projection from M onto the (qi , pi )-coordinate plane, denoted (Qi , Pi ), then i (dPi dQi ) = dpi dqi and so (dpi dqi ) = = = (dpi dqi ) = D D D i (dPi dQi ) (i ) (dPi dQi ) (dPi dQi ) = Ai i () , i where Ai denotes the oriented area in the (Qi , Pi )-plane. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 60 / 104 Symplectic Integration Symplectic structures If we let = dp dq be the canonical symplectic form on M = T U, then it follows that = = = n i=1 n i=1 n i=1 (dpi dqi ) (dpi dqi ) dPi dQi = n i=1 i Ai i () and so the canonical surface area invariant is just the sum of the oriented areas of the projections of into each of the coordinate planes (qi , pi ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 61 / 104 Symplectic Integration Symplectic structures Symplectic Maps of Euclidean Phase Space Example: If is the canonical symplectic form on the phase space M = T U and : M M is a smooth map, then is symplectic if and only if for every z = (q, p) M and all v, w Tz M we have vT J T w = (z) v, w = (z)(v, w) = (z)( v, w) = ( (z)v, (z)w) T = (z)v J T (z)w) = vT [ (z)]T JT (z) w, [ (z)]T J T (z) = J T . i.e., is symplectic if and only if for every z M, we have Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 62 / 104 Symplectic Integration Symplectic structures Alternative Test of Symplecticness Symplecticness can also be checked directly using differential forms. If = dp dq is the canonical symplectic form on T (U) and : M M can be expressed as ^ q = 1 (q, p) then is symplectic if and only if ^ ^ d p d q = dp dq, where 1 1 ^ d q = q dq + p dp 2 2 ^ d p = q dq + p dp. and ^ p = 2 (q, p), and Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 63 / 104 Symplectic Integration Symplectic structures Canonical Hamiltonian Systems Suppose that H(z) is a Hamiltonian on the phase space T (U) U Rn = and that the conjugate variables q and p evolve according to the Hamiltonian system: q = p H(q, p) p = -q H(q, p). Then, writing z = (q, p) and letting J be the canonical symplectic structure matrix introduced previously, this system of equations can be rewritten as z = Jz H(z). with solution z(t; z0 ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 64 / 104 Symplectic Integration Symplectic structures Variational Equations for a Hamiltonian System Suppose that the flow map (t, z0 ) = z(t; z0 ) is a diffeomorphism for every t. If F(t) = z0 t (z0 ) is the Jacobian of t , then d F(t) = dt t (z0 ) = t (z0 ) t z0 z0 t = [Jz H(t (z0 ))] z0 = JHz z (z(t)) t (z0 ) z0 = JHz z (z(t))F(t). The resulting matrix ODE's are the variational equations along z(t). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 65 / 104 Symplectic Integration Symplectic structures Hamiltonian Flows are Canonically Symplectic Theorem: The flow map t is canonically symplectic. Proof: We need to verify that the Jacobian F(t) = z0 t (z0 ) satisfies the condition F(t)T JT F(t) = JT . Since F(0) = I2n , this statement is true for t = 0. Thus it suffices to show that d K(t) = 0; dt K(t) F(t)T JT F(t). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 66 / 104 Symplectic Integration Symplectic structures However, by using the variational equations for F(t), we find K = FT JT F + [F]T JT F T T = F J JHz z (z(t))F + FT Hz z (z(t))JT JT F = FT Hz z (z(t))F - FT Hz z (z(t))F = 0, where we have used the fact that JT J = I2n and JT JT = -I2n . Remark: It follows that the Hamiltonian flow preserves not only the phase space volume, but also the oriented surface area of the projections of any 2-cell onto the (qi , pi )-coordinate planes. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 67 / 104 Symplectic Integration Symplectic Integration One-Step Integrators Suppose that t (z0 ) is the flow map of a Hamiltonian system on M = T (U). We will say that a numerical solution (zn ; n 0) for this system is a one-step method with step size h if there is a family of maps h : M M with h (0, h0 ) such that n+1 zn+1 = h (zn ) = h (z0 ) (n+1)h (z0 ). Recall that: h is consistent if limh0 h-1 ||h (z0 ) - h (z0 )|| = 0; h is of order r if limh0 h-r ||h (z0 ) - h (z0 )|| = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 68 / 104 Symplectic Integration Symplectic Integration Symplectic Integrators Definition: The one-step method determined by h is said to be a symplectic integrator if for each h (0, h0 ), the map h : M M is canonically symplectic: T h (z) J T h (z) = J T . z z Motivation: Because symplectic maps preserve the oriented area of surfaces projected onto conjugate coordinate planes, numerical solutions produced using a symplectic integrator will share many of the same integral invariants as the actual solution to the Hamiltonian system. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 69 / 104 Symplectic Integration Symplectic Integration Example: Asymmetrical Euler Method Given the Hamiltonian H, consider the following asymmetrical Euler scheme: qn+1 = qn + hp H(qn , pn+1 ) pn+1 = pn - hq H(qn , pn+1 ). Since this method is implicit, we first need to show that for sufficiently small h > 0 there is a smooth map h such that (qn+1 , pn+1 ) = h (qn , pn ) is the unique solution to the scheme in some open neighborhood of (qn , pn ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 70 / 104 Symplectic Integration Symplectic Integration To verify the existence and uniqueness of such a solution, let F : Rn Rn R Rn Rn be the function F (qn , pn , h, u) = u - pn - hq H(qn , u), and observe that F (qn , pn , 0, pn ) = 0 u F (qn , pn , 0, pn ) = In . Since F has full rank at (qn , pn , 0), it follows from the implicit function theorem that there is a smooth function pn+1 = u(qn , pn , h) with pn = u(qn , pn , 0) such that for sufficiently small h > 0 F qn , pn , h, u(qn , pn , h) = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 71 / 104 Symplectic Integration Symplectic Integration We can then define h locally by setting: q h (q, p) = q + hp H(q, u(q, p, h)) p h (q, p) = u(q, p, h) which is smooth because u and H are both smooth mappings. Thus the asymmetrical Euler method is well-defined. Our next aim is to show that this method is canonically symplectic: dpn+1 dqn+1 = dpn dqn where we recall the notation dp n+1 dq n+1 = n i=1 (dpin+1 dqin+1 ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 72 / 104 Symplectic Integration Symplectic Integration We first observe that d q H(qn , pn+1 ) d Hq1 (qn , pn+1 ), , Hqn (qn , pn+1 ) n n = (Hq1 qj dqjn + Hq1 pj dpjn+1 ), , (Hqn qj dqjn + Hqn pj dpjn+1 ) n n = Hq1 qj dqjn , , Hqn qj dqjn + j=1 j=1 j=1 j=1 j=1 j=1 = Hqq dq + Hqp dp n n Hq1 pj dpjn+1 , , Hqn pj dpjn+1 n n+1 . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 73 / 104 Symplectic Integration Symplectic Integration Similarly, d p H(qn , pn+1 ) Hpq dqn + Hpp dpn+1 , T where Hpq = Hqp . Using these two results along with the implicit definitions of qn+1 and qn+1 shows that dqn+1 = d qn + hp H(qn , pn+1 ) = dqn + h Hqp dqn + Hpp dpn+1 dpn+1 = d pn - hq H(qn , pn+1 ) = dpn - h Hqq dqn + Hpq dpn+1 . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 74 / 104 Symplectic Integration Symplectic Integration The next result can be used to simplify the product dpn+1 dqn+1 : Lemma: If dx and dy are k-vectors of 1-forms and A is a k-by-k matrix, then dx (Ady) = (AT dx) dy where the exterior product of two vectors of 1-forms is the sum of the component-wise exterior products. In particular, we have: dpn+1 Hpp dpn+1 = Hpp dpn+1 dpn+1 dpn+1 Hqp dqn = Hpq dpn+1 dqn Jay Taylor (ASU) APM 530 - Lecture 8 dqn Hqq dqn = Hqq dqn dqn = 0 = -dpn+1 Hpp dpn+1 = 0 Fall 2010 75 / 104 Symplectic Integration Symplectic Integration We can then compute dpn+1 dqn+1 = = dpn+1 dqn + h dpn+1 Hqp dqn + Hpp dpn+1 = dpn dqn - h Hqq dqn + Hpq dpn+1 dqn + h dpn+1 Hqp dqn + Hpp dpn+1 = dpn dqn , = dpn dqn - h Hpq dpn+1 dqn + h Hpq dpn+1 dqn which shows that the asymmetric Euler method is symplectic. A direct calculation also shows that this method is first-order accurate. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 76 / 104 Symplectic Integration Symplectic Integration Symplectic Runge-Kutta Methods The general form of a Runge-Kutta method for a Hamiltonian system is Qi Pi = q +h = pn + h n s j=1 s j=1 s i=1 s i=1 aij Fj aij Gj bi Fi , bi Gi , i = 1, , s, i = 1, , s, qn+1 = qn + h pn+1 = pn + h where s is the number of stages and (Qi , Pi ) are the internal stage variables. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 77 / 104 Symplectic Integration Symplectic Integration Remarks: The functions Fi and Gi are given by Fi Gi = p H(Qi , Pi ) = -q H(Qi , Pi ). s i=1 This method is consistent if and only if bi = 1. This method is implicit except when the matrix (aij ) is strictly lower triangular: aij = 0 for j i. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 78 / 104 Symplectic Integration Symplectic Integration Theorem: The Runge-Kutta method is symplectic if and only if bi aij + bj aji - bi bj = 0, i, j = 1, , s. Proof: It will be convenient to introduce tensor notation: if T = (tij ) is an m n matrix and S = (sij ) is a p q matrix, then T S is the mp nq matrix defined by t11 S t12 S t1n S t21 S t22 S t2n S TS= . . . . . . . . . . tm1 S tm2 S tmn S Mixed-product identity: If the products RT and QS are defined, then (R Q)(T S) = (RT QS). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 79 / 104 Symplectic Integration Symplectic Integration The RK method can be rewritten using tensor notation as Q = e qn + h(A I)F qn+1 = qn + h(bT I)F where A = (aij ), b = (bi ) are the step coefficients; e = (1, 1, , 1)T Rs ; I is the n n identity matrix; Q = (Q1 , , Qs )T Rsn , P = (P1 , , Ps )T Rsn ; F = (F1 , , Fs )T Rsn , G = (G1 , , Gs )T Rsn . P = e pn + h(A I)G pn+1 = pn + h(bT I)G, Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 80 / 104 Symplectic Integration Symplectic Integration Likewise, we can also express the differential forms of the RK transformation using tensor notation: dQ = e dqn + h(A I)dF dqn+1 = dqn + h(bT I)dF where dF = FQ dQ + FP dP and dG = GQ dQ + GP dP. dP = e dpn + h(A I)dG dpn+1 = dpn + h(bT I)dG, Furthermore, because F and G are determined by a Hamiltonian system, we also have FQ = -GT , P Jay Taylor (ASU) FP = FT , P APM 530 - Lecture 8 GQ = GT . Q Fall 2010 81 / 104 Symplectic Integration Symplectic Integration Taking exterior products gives dpn+1 dqn+1 = = dpn + h(bT I)dG dqn + h(bT I)dF + h2 (bT I)dG (bT I)dF + h2 dG (bbT I)dF, = dpn dqn + h dpn (bT I)dF + h (bT I)dG dqn = dpn dqn + h dpn (bT I)dF - h dqn (bT I)dG where we have used the identity (C D)T = CT DT . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 82 / 104 Symplectic Integration Symplectic Integration If we let B be the n n diagonal matrix with Bii = bi , then dQ (B I)dG = = e dqn + h(A I)dF (B I)dG and dP (B I)dF = = e dpn + h(A I)dG (B I)dF = dqn (eT B I)dG + h (BT A I)dF dG = dqn (bT I)dG - h dG (BA I)dF = dpn (eT B I)dF + h dG (AT B I)dF = dpn (bT I)dG + h dG (AT B I)dF. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 83 / 104 Symplectic Integration Symplectic Integration Summarizing, we have derived the following identities: dpn+1 dqn+1 = dpn dqn + h dpn (bT I)dF - dQ (B I)dG = dqn (bT I)dG - h dG (BA I)dF h dqn (bT I)dG + h2 dG (bbT I)dF dP (B I)dF = dpn (bT I)dF + h dG (AT B I)dF, which can be combined to give dpn+1 dqn+1 = h2 dG bbT - BA - AT B I dF + +h dP (B I)dF - h dQ (B I)dG. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 84 / 104 Symplectic Integration Symplectic Integration To make further progress, we must exploit the fact that F and G are determined by a Hamiltonian vector field: dP (B I)dF - dQ (B I)dG = = dP (B I) (FQ dQ + FP dP) = dP (B I)FP dP - dQ (B I)GQ dQ - dQ (B I) (GQ dQ + GP dP) = dP (B I)FQ dQ + dP (BT I)GT dQ P T = dP (B I) FQ + GP dQ = 0 since FQ = -GT . P Jay Taylor (ASU) APM 530 - Lecture 8 + dP (B I)FQ dQ - dQ (B I)GP dP Fall 2010 85 / 104 Symplectic Integration Symplectic Integration It follows that dpn+1 dqn+1 - dpn dqn = = h2 dG bbT - BA - AT B I dF bbT - BA - AT B = 0 holds, i.e., if and only bi bj - bi aij - bj aji = 0 for all i, j = 1, , s. which vanishes if and only if the matrix identity Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 86 / 104 Symplectic Integration Symplectic Integration Corollary: Symplectic Runge-Kutta methods are necessarily implicit. Proof: Taking j = i in the symplecticness condition gives 0 = bi2 - 2bi aii = bi (bi - 2aii ) and so aii = 0 whenever bi = 0. Since at least one of the bi 's is non-negative, it follows that (aij ) is not strictly lower triangular. One of the shortcomings of the RK methods considered so far are that we have adopted the same numerical scheme for the position and momentum coordinates. Methods that are both explicit and symplectic can be formulated by applying a different RK scheme to each set of coordinates. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 87 / 104 Symplectic Integration Symplectic Integration Partitioned Runge-Kutta Methods The general form of a partitioned Runge-Kutta method for a Hamiltonian system is: Qi Pi = qn + h = pn + h s j=1 s j=1 s i=1 s i=1 aij Fj ~ij Gj a bi Fi , ~ bi Gi . i = 1, , s, i = 1, , s, qn+1 = qn + h pn+1 = pn + h Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 88 / 104 Symplectic Integration Symplectic Integration Theorem: The partitioned Runge-Kutta method is symplectic if and only if ~ ~ bi ~ij + bj aji - bi bj = 0, i, j = 1, , s. a Remarks: This can be proven using differential forms as before. Consistency requires s i=1 s i=1 bi = ~ bi = 1. Taking i = j now gives ~ ~ bi ~ii + bi aii - bi bi = 0 a ~ which can be satisfied if bi bi = 0. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 89 / 104 Symplectic Integration Symplectic Integration Conservation of Quadratic Invariants Theorem: A symplectic Runge-Kutta scheme preserves any quadratic invariant of the form I = qT Wp + dT q + dT p. 1 2 Sketch of Proof: For simplicity, we will assume that d1 = d2 = 0. Because the method is symplectic, we know that BA + AT B - bbT = 0, where A = (aij ) and B is the diagonal matrix with Bii = bi . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 90 / 104 Symplectic Integration Symplectic Integration As above, we can write Q = e qn + h(A I)F qn+1 = qn + h(bT I)F which shows that (qn+1 )T Wpn+1 = (qn )T Wpn + h (qn )T (bT W)G and QT (B W)G = (qn )T (bT W)G + hFT (AT B W)G FT (B W)P = FT (b W)pn + hFT (BAT W)G. APM 530 - Lecture 8 P = e pn + h(A I)G pn+1 = pn + h(bT I)G, + h FT (b W)pn + h2 FT (bbT W)G Jay Taylor (ASU) Fall 2010 91 / 104 Symplectic Integration Symplectic Integration Combining these expressions gives (qn+1 )T Wpn+1 = (qn )T Wpn + hFT bbT - BA - AT B W G = (qn )T Wpn + h QT (B W)G + h FT (B W)P since symplecticness implies that the second term vanishes. However, if I is invariant for the Hamiltonian flow, then 0 = I (q, p) = qT Wp + qT Wp = (p H(q, p))T Wp + qT W(-q H(q, p)). + h QT (B W)G + h FT (B W)P Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 92 / 104 Symplectic Integration Symplectic Integration In particular, by taking q = Qi and p = Pi and recalling that Fi Gi this last identity gives FT WPi + QT WGi = 0 i i which then implies that QT (B W)G + FT (B W)P = 0. = p H(Qi , Pi ) = -q H(Qi , Pi ), Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 93 / 104 Symplectic Integration Symplectic Integration This allows us to conclude that (qn+1 )T Wpn+1 = (qn )T Wpn and so I is conserved by the RK map. A similar calculation shows that symplectic partitioned RK methods conserve quadratic invariants. If H is translation invariant, then the total linear momentum is a conserved quantity, i.e., I = dT p is invariant when d = (1, , 1). If H is rotation invariant, then angular momentum is conserved and so I = (qi pi ) a i is invariant for any a R3 . Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 94 / 104 Symplectic Integration Symplectic Integration The position Verlet method is symplectic Recall that the position Verlet scheme (expressed in (q, p)-coordinates) is qn+1/2 = qn + pn+1 qn+1 h -1 n M p 2 = pn + h F(qn+1/2 ) h = qn+1/2 + M-1 pn+1 2 h -1 n n = q + M p + pn+1 2 Remark: This scheme is explicit. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 95 / 104 Symplectic Integration Symplectic Integration The position Verlet method can be expressed as a partitioned Runge-Kutta method with two stages: Q1 = P1 = qn+1 = pn+1 = giving ~ ~ b1 = b2 = 1/2; b1 = 0, b2 = 1 a11 = a12 = 0; a21 = 1/2, a22 = 0 ~11 = ~12 = ~21 = 0; ~22 = 1. a a a a Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 96 / 104 1 -1 q + 0, Q2 = q + h M P1 + 0 M-1 P2 2 pn + 0, P2 = pn + h 0 + F(Q2 ) 1 -1 1 -1 n q +h M P1 + M P2 2 2 pn + h 0 F(Q1 ) + F(Q2 ) , n n Symplectic Integration Symplectic Integration ~ ~ There are four conditions to be checked: bi ~ij + bj aji - bi bj = 0 a (1, 1) : (1, 2) : (2, 1) : (2, 2) : 1 2 1 2 1 2 1 2 1 0=0 2 1 1 0+1 - 1=0 2 2 1 0+00- 0=0 2 1 1 + 1 0 - 1 = 0, 2 0+00- and since these are all satisfied, it follows that the position Verlet method is symplectic. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 97 / 104 Symplectic Integration Symplectic Integration The Velocity Verlet Method The velocity Verlet method updates the velocity at the midpoint of each time interval: h pn+1/2 = pn + F(qn ) 2 qn+1 = qn + h M-1 pn+1/2 h pn+1 = pn+1/2 + F(qn+1 ) 2 h n = p + F(qn ) + F(qn+1 ) 2 A similar calculation shows that this method is symplectic. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 98 / 104 Symplectic Integration Symplectic Integration Splitting and Composition In some cases, symplectic methods can be constructed by composing the flow maps for separable components of a Hamiltonian system. Suppose that H = H1 + + Hk is a splitting of the Hamiltonian and that t,i is the flow map corresponding to Hi : d i,t = Jz Hi . dt Then h = 1,h k,h is a consistent symplectic method for the full Hamiltonian system. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 99 / 104 Symplectic Integration Symplectic Integration Example: Suppose that the Hamiltonian is separable: H(q, p) = T (p) + V (q) and set H1 (q, p) = T (p) and H2 (q, p) = V (q). Then the Hamiltonian system corresponding to H2 can be written as q = 0 p = -q V (q) which is completely integrable with linear flow map q V ,t (q, p) = . p - tq V (q) Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 100 / 104 Symplectic Integration Symplectic Integration Similarly, the system corresponding to H1 q = p T (p) p = 0 is completely integrable with linear flow map q + tp T (p) T ,t (q, p) = . p The composition of these two maps is h (q, p) = T ,h V ,h (q, p) q = T ,h p - hq V (q) q + hp T p - hq V (q) = . p - hq V (q) Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 101 / 104 Symplectic Integration Symplectic Integration Since T ,h and V ,h are both symplectic maps, it follows that h gives rise to a symplectic numerical scheme for H which can be written as qn+1 = qn + hp T (pn+1 ) pn+1 = pn - hq V (qn ). This is the asymmetrical Euler method specialized to a mechanical Hamiltonian. Of course, we will also obtain a symplectic method by composing these ^ maps in the opposite order: h = V ,h T ,h . This leads to an alternative asymmetrical Euler scheme: qn+1 = qn + hp T (pn ) pn+1 = pn - hq V (qn+1 ). Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 102 / 104 Symplectic Integration Symplectic Integration Example: Suppose that H = T + V is separable and consider the splitting H = H1 + H2 + H3 with 1 H1 = V (q), 2 H2 = T (p), 1 H1 = V (q). 2 The composition method associated with this splitting is h = V ,h/2 T ,h V ,h/2 which leads to the recursion h pn+1/2 = pn - q V (qn ) 2 n+1 n q = q + h p T (pn+1/2 ) h pn+1 = pn+1/2 - q V (qn ) 2 which we recognize as the velocity Verlet method. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 103 / 104 Symplectic Integration References References Arnold, V.I. (1978) Mathematical Methods of Classical Mechanics., trans. by K. Vogtman and A. Weinstein. Springer-Verlag. Barden, D. and Thomas, C. (2003) An Introduction to Differential Manifolds. Imperial College Press. Leimkuhler, B. and Reich, S. (2004) Simulating Hamiltonian Dynamics. Cambridge Monographs on Applied and Computational Mathematics, v. 14. Sanz-Serna, J.M. (1991) Symplectic integrators for Hamiltonian problems: an overview. Acta Numerica 1: 243-286. Schlick, T. (2006) Molecular Modeling and Simulation. Springer. Jay Taylor (ASU) APM 530 - Lecture 8 Fall 2010 104 / 104 ...
View Full Document

This note was uploaded on 03/11/2012 for the course APM 530 taught by Professor Staff during the Fall '10 term at ASU.

Ask a homework question - tutors are online