Integration of Ordinary Differential Equations and Chaos in Dynamics

Chapter 16 Integration of Ordinary Differential Equations

Examples of Differential Equations ODE: Chemical reaction dynamics, dC/dt = -C Pressure change with height, dp/dy = -mgp/kT Newton s equation of motion, F = md2r/dt2 PDE: Maxwell s equations for electricity and magnetism Structure and fluid mechanics Schr dinger equation in quantum mechanics

Higher ODE Reduces to 1stOrder 2 d y dx dy dx dy dx + + + = = 2 '''( ) Define dy dx dy dx dy dx ( ) ''( ) y x y x y y = 2( '( )) ', y = ( ) '', we have 0 y x y x y x y + = ( ) q x ( ) r x 2 , y y 1 2 3 = Define , we have z = 1 y 2 dy dx dz dx = z = 2 y 3 = ( ) r x ( ) q x z = 2 3 2( ) y y y y 1 3 2 1 In general, it is sufficient to solve coupled first-order ordinary differential equations of the form dy dx= = ( , , , ), 1,2, , i f x y y i N 1 i N

Initial Value Problem It is convenient to consider the independent variable x as time t. The solution to the equations dy f t y dt = = ( , , , ), 1,2, , i y i N 1 i N is uniquely determined if the initial value at t=0, yi(0), is given. The equation can be written in vector form Y d dt = F ( , ) t Y

Some General Properties of Autonomous Systems F(t,Y) = F(Y), independent of time t The space spanned by Y (a set of all possible Y) is called phase space F forms a vector field (a vector at each point Y) y2 Solution of dY/dt = F(Y) produces a parametric curve Y(t) in phase space. Intersection of trajectories cannot happen, why? F y1

Fixed Points A location in phase space such that F(Y)=0. Attractor, repellor Saddle point or hyperbolic fixed point

Chaos Extremely sensitive to initial conditions [dY(t) = exp( t)dY(0)]. E.g., Lorenz s weather model: y3 = = = + ( ) y y y y y y y 1 1 2 ry y + 2 1 2 1 3 by y y 3 3 1 2 y2 y1

Finite difference + ( ) dx ( ) ( ) y x dy x y x h Forward difference: = + ( ) O h h ( ) dx ( ) y x ( ) dy x y x h h = + Backward difference: ( ) O h + ( ) dx ( ) ( ) dy x y x h y x h = + 2 ( ) O h Central difference: 2 h + Y Y Y ( ) x dx Y ( ) ( ) x d x h Euler Method: = F Y ( , ) x h F = + Y Y or ( , ) h x + 1 n n n n

Euler and Midpoint += + + 2 ( , ) ( ) y y h f x y O h 1 n n n n h = ( , ) 1 2 + k h f x y 1 n n 1 2 = + + ( , ) k h f x h y k 2 1 n n = + 3 ( ) y y k O h + 1 2 n n

4-th Order Runge-Kutta Method = ( , ) k h f x y 1 n n k h = + + ( , ) 1 k h f x y 2 n n 2 h 2 k = + + ( , ) 2 k h f x y 3 n n 2 h y 2 = + + ( , ) k h f x k k 4 3 n n k k k = + + + + + 5 3 ( ) 1 2 4 y y O h + 1 n n 6 3 3 6

rk4( ) Python code for ODE from scipy.integrate import ode solver = ode(f).set_integrator( dopri5 ) solver.set_initial_value(y0,t0) solver.integrate(t+dt)

Some General Concepts Discretized equation, such as yn+1=yn+hf(xn,yn), is consistent, if as h->0, it approaches the original differential equation The error |y(xn+1)-yn+1| =O(hk) in one step from xn to xn+1 is called local truncation error The error |y(x)-yn| for some finite x and initial condition y(0) = y0 is the global error The method is convergent if the global error goes to zero as h -> 0 and n -> .

Adaptive Stepsize Control Estimate local truncation error from difference between one h step and two steps of h/2 Or difference of 4 and 5-th order Runge- Kutta Increase h if error is small than tolerance, decrease h if error is bigger than tolerance. See NR p.721, odeint() for details.

Richardson Extrapolation and Bulirsch-Stoer Method Take a large step size H, consider the answer as an analytic function f(h) of h=H/n. Fit the function by polynomial or rational function interpolation. Choose a method (e.g., midpoint) such that f(h) is even in h. And finally extrapolate to h=0.

Multi-step, Explicit, Implicit, etc Solving equation y =f(x,y) is to compute x y x y f t y t dt = + ( ) ( , ( )) n x n In general, this results in = + ( ) n n n + + + y y h y y y + + 1 0 1 1 2 1 n n where y = ( , f x y ) n n n

Hamiltonian System The system of equations = = H q = + = , ( , ) H p q p T V pq L i i H p = , 1,2, , q i N i i has special properties. It is equivalent to Newton s equation with a potential energy. V q = = , mq F F i i i i

Verlet or Strmer Algorithm Solve 2 d q dt j = ( ) F q j 2 by central difference + ) 2 ( ) h + = + 2 4 ( ( ) ( ( )) ( ) q t q t q t h h F q t O h j j j j Show the error is O(h4)

2-Form and Symplectics The Hamiltonian dynamics, beside having a conserved energy, also has additional conserved quantities ( 2)n,n=1,2,..,N: wedge product: = 2 = = dp dq 0, dx dx dx dx i i j j dx dx j i j j i A canonical transform is a mapping from (p,q) to (P,Q) such that the form of 2 is the same. I.e. j j becomes dp dq dP dQ j j j j

Canonical Transformation Equivalent condition for canonical mapping z to Z is D JD J = Note: dz T dt= J H where p q 0 I Z z I ( ) z i = = = , [ ] , i J D z ij 0 i j 2N means volume element in phase space Hamiltonian dynamics preserves the volume Liouville s theorem.

Example of Symplectic Algorithm Euler method is not symplectic p p hF q q hp + = + 2 = + p = + assuming ( , ) ( ) H p q V q + 1 n n n 2 m / m 1 n n n But the following is p p q q + = = + Let ( ) P p hF q 2 h m h m h m = + = + + ( ) Q q P p q F q = + + hF hp + 1 n n n then 1/ + m 2 h m h m ( ) = + + + 1 n n n '( ) 1 '( ) dP dQ dp hF q dq dp F q dq 2 2 h m dp h m = + + '( ) 1 '( ) F q dq dp dp F q dq = dq

Second-Order Symplectic or Velocity Verlet Combine two half-step size first-order symplectic algorithms, one can obtain: h q q p m h p p F + = + + 2 h = + + F + 1 n n n n 2 m F + 1 1 n n n n 2 Symplectic algorithm preserves the symplectic properties of the Hamiltonian system exactly.

Deriving Symplectic Algorithm The equation of motion of any dynamic variable A=A(p,q) is dA/dt = - {H,A}. The formal solution is A(t) = exp(-t{H,.})A. Since H=T+V, we can perform operator splitting (Trotter-Suzuki), A(t) exp(-t{T,.}) exp(-t{V,.}) A. Since exp(-{G,.}) is symplectic, the operator splitting is symplectic provided each factor is computed exactly.

Problem set 10 (14 Nov 24) 1. Show that the last 2nd order symplectic algorithm is indeed symplectic! 2. Show that the 4-th order Runge-Kutta is equivalent to Simpson rule if y =f(x,y)=f(x) independent of y. 3. Verify that the 4-th order Runge-Kutta formula is indeed accurate to 4-th order [Taylor expanding both side of equation (Fig. 16.1.3)]. Do this with Mathematica.