Unformatted text preview: LeapFrog and Verlet method for ODEs Example: Simple Pendulum ODE of vectors Example: Projectile motion with air drag HW#6 due today Get HW#6 off website Reading for Differential Equations in Appendix B. Consider general ODE with initial condition: We are interested in a solution from t=0 to t=T in n time steps: du = f (u(t)) and u(0) = u0 dt T ∆t = n Notation for solution at time step k: Using foward difference, ODE becomes Finally solution for time k+1 is uk+1 − uk = f (uk ) ∆t uk+1 = uk + ∆tf (uk ) uk ≡ u(tk ) Need value of u at k=0 time step to get this started (initial condition). In many systems, the rate at which a quantity changes depends linearly on the value of that quantity. du (chemical reactions, nuclear fission, = −bu dt heat transfer,…) Initial condition (t=0): u(0) = A
Solution is an exponential: u(t) = Ae −bt Euler is based on forward difference approximation with error of order τ 2 . Big problem is it is NOT energy conserving. Using center difference method, we get the leap-frog method: Error is of order τ 3 , and it IS energy conserving. However it is not self starting (need k=-1 step at start). One solution is to take an initial Euler step, then go with LeapFrog after the first step is computed. Why “leap-frog”? Let’s apply to a standard mechanics problem… uk+1 = uk−1 + 2τ f (uk ) By Newton’s laws and basic mechanics force (F), acceleration (a) velocity (v), position (r): F = ma
dv = a(r(t)) dt dr = v (t) dt Apply Leap-Frog method to solve for position and velocity at time step k based on previous values: vk+1 = vk−1 + 2τ a(rk ) rk+2 = rk + 2τ vk+1
Note how v and r values leap-frog each other. Here we use the central difference formulas for BOTH the first derivative (v) AND second derivative (r). 2 dr dr = a(r) 2 dt dt = v (t) rk+1 − rk−1 So scheme becomes: v k= 2τ rk+1 = 2rk − rk−1 + τ 2 ak
Error for v is of order τ 2 and error for r is of order Calculation of v is not required to compute r. τ 4 Physics – velocity becomes angular velocity and position becomes angular position. Equations of motion are: dφ =ω dt dω g = α(φ) = − sin φ dt l Or numerically using Verlet (just interested in position) φk+1 = 2φk − φk−1 + τ αk
2 Euler Verlet Often in science the quantities we are interested in are vectors which require a magnitude and direction. Typically we express a vector by it’s individual components: = (vx , vy , vz ) v
Computationally, our 1D problems now become 3 1-D problems. We just perform the calculations on each of the components individually. dvx 1 = Fdx (vx ) dt m dvy 1 = Fdy (vy ) − g dt m Restrict to 2-D, so only need x and y components of r and v. Force of gravity is only in y direction, drag force is always opposite to velocity (x and y components constantly change). drx = vx dt dry = vy dt 1 Fd = − Cd ρA| | vv 2 ...
View Full Document