Ordinary Differential Equations Ordinary Differential Equations (ODEs)(ODEs)NRiC Chapter 16.ODEs involve derivatives wrt one independent variable, e.g. time t.ODEs can always be reduced to a set of first-order equations (involving only first derivatives).e.g.is equivalent to the set =  =   =  −  ODE Intro, Cont'dODE Intro, Cont'dUsually new variables just derivatives of old, but sometimes need additional factors of t to avoid pathologies.General problem is solving set of 1st order ODEs:But also need boundary conditions: algebraic conditions on values of yi at discrete time(s) t...=         =    fi' are known functionsODE Boundary Conditions (BC)ODE Boundary Conditions (BC)Two categories of BC:Initial Value Problem (IVP): all yi are given at some starting point ts, and solution is needed from ts to tf.Two-point Boundary Value Problem (BVP): yi are specified at two or more t, e.g. some at ts, some at tf (only one BC needed for each y).Generally, IVP much easier to solve than 2-pt BVP, so consider this first.Finite DifferencesFinite DifferencesHow do you represent derivatives with discrete number system?Basic idea: replace dy/dt with finite differences y/t. Then:How do you use this to solve ODEs?   Euler's MethodEuler's MethodWrite y/t = f'(t,y) ⇒ y = t f'(t,y).Start with known values yn at tn (initial values).Then yn+1 at tn+1 = tn + h is: yn+1 = yn + h f'(tn,yn)h is called the step size.Integration is not symmetric: derivative evaluated only at start of step ⇒ error term O(h2), from Taylor series. So, Euler's method is first order.Euler's Method, Cont'dEuler's Method, Cont'dRunge-Kutta MethodsRunge-Kutta MethodsWe can do better by symmetrizing derivative:Take a trial step to midpoint, evaluate yn+1/2 and tn+1/2.Use these to evaluate derivative f'(tn+1/2,yn+1/2).Then use this to go back and take a full step.Thus: yn+1 = yn + h f'(tn + ½ h,yn + ½ h f'(tn,yn)) + O(h3)Can show that O(h2) terms "cancel," so leading error term is O(h3) ⇒ 2 nd-order Runge-Kutta.Runga-Kutta, Cont'dRunga-Kutta, Cont'dFourth-Order Runge-KuttaFourth-Order Runge-KuttaActually, there are many ways to evaluate f' at midpoints, which add higher order error terms with different coefficients. Can add these together in ways such that higher order error terms cancel. e.g. can build fourth-order Runge-Kutta (RK4): k1 = h f'(tn,yn) k2 = h f'(tn + h/2,yn + k1/2) k3 = h f'(tn + h/2,yn + k2/2) k4 = h f'(tn + h,yn + k3)Then yn+1 = yn + k1/6 + k2/3 + k3/3 + k4/6 + O(h5)Fourth-Order Runge-Kutta, Cont'dFourth-Order Runge-Kutta, Cont'dFourth-Order Runge-Kutta, Cont'dFourth-Order Runge-Kutta, Cont'dDisadvantage of RK4: requires f' to be evaluated 4 times per step.But, can still be cost effective if larger steps OK.RK4 is workhorse method. Higher-order RK4 takes too much effort for increased accuracy.Other methods (e.g. Bulirsch-Stoer, NRiC §16.4) are more accurate for smooth functions.But RK4 often "good enough".The Leapfrog IntegratorThe Leapfrog IntegratorExtremely useful for second-order DEs in which d2x/dt2 = f(x), e.g. SHM, N-body, etc.Suppose x is position, so d2x/dt2 is acceleration.Procedure: define v = dx/dt at the midpoints of the steps, i.e. velocities staggered wrt positions.Define vi+1/2 = v(t + ½ t), i = 0, 1, 2, ...Then advance xi to xi+1 and vi+1/2 to vi+3/2: xi+1 = xi + vi+1/2 tvi+3/2 = vi+1/2 + f(xi+1) tLeapfrog, Cont'dLeapfrog, Cont'dLeapfrog, Cont'dLeapfrog, Cont'dComplication: need to "jumpstart" and "resync"...vi+1/2 = vi + ½ t f(xi) [opening "kick": Euler] xi+1 = xi + vi+1/2 t ["drift"] vi+1 = vi+1/2 + ½ t f(xi+1) [closing "kick": resync] Note vi+3/2 = vi+1 + ½ t f(xi+1) = vi+1/2 + t f(xi+1).Also have "drift-kick-drift" (DKD) scheme.Like midpoint method, Leapfrog is second order.So why is Leapfrog so great?...Leapfrog, Cont'dLeapfrog, Cont'dAnswer: Leapfrog is time reversible.Suppose we step back from (ti+1,xi+1,vi+3/2) to (ti,xi,vi+1/2). Applying the algorithm: vi+1/2 = vi+3/2 + f(xi+1)(- t) xi = xi+1 + vi+1/2(-t) These are precisely the steps (in reverse) that we took to advance the system in the first place!Hence if we Leapfrog forward in time, then reverse to t = 0, we're back to where we started, precisely.Leapfrog, Cont'dLeapfrog, Cont'dLeapfrog is time reversible because of the symmetric way in which it is defined, unlike the other schemes.In Euler, forward and backward steps do not cancel since they use different derivatives at different times.In Midpoint, the estimate of the derivative is based on an extrapolation from the left-hand side of the interval. On time reversal, the estimate would be based on the right-hand side, not the same.Similarly, RK4 is not time reversible.Leapfrog, Cont'dLeapfrog, Cont'dTime reversibility is important because it guarantees conservation of energy, angular momentum, etc. (in many cases).Suppose the integrator makes an error  after one orbital period. Now reverse. Is the error -? No! The time-reversed orbit is a solution of the original ODE (with v replaced with -v), so the energy error is still +. But we've returned to our starting point, so we know the final energy error is zero. Hence  = 0!Leapfrog is only second order, but very stable.Adaptive Stepsize ControlAdaptive Stepsize ControlUp to now, have assumed stepsize h is constant.Clearly prefer choosing h small when f ' is large, and h large when f ' is small.The tradeoff is extra trial steps to determine optimum h, but may achieve factor of 10 to 100 increase in stepsize, so it's often worth it.NRiC provides a routine odeint() for RK4 with adapative stepsize control. Complicated to