COMP155 / EMGT155Dec 1 2008Image from http://en.wikipedia.org/wiki/Integral We’ll commonly have a derivative of some function(or instantaneous values of the derivate)and need to construct values for the function Example: We know velocity, need to compute positionneed to compute position Example: We know acceleration, need to compute velocity and position Simulation is updated at regular interval in time dt = time step/slice To approximate integration:compute derivative at time stepscompute derivative at time steps make assumption about derivative value at other times accumulate area to approximate integration ∆t = 1, assume current value held for previous time sliceblue: correct derivative functionover-estimate of positive functionred: approximated derivative functionunder-estimate of positive areaover-estimate of negative areaof positive areaunder-estimate of negative area∆t = 1blue actual velocityred red approx velocitygreen actual positionmagenta computed position∆t = 0.5∆t = 1∆t = 0.5∆t = 1∆t = 0.1∆t = 0.01euler_pos_vel.pyDefinition of a derivative:taftafaft∆−∆+=→∆)()(lim)('0∆t = 1As ∆t gets smaller, we get closer to the true value of the derivative.∆t = 1∆t = 0.1Definition of derivative:Forward Euler Method:ttfttftft∆−∆+=→∆)()(lim)('0∆t = 1Forward Euler Method:Backward Euler Method:(we’re using backward Euler))(')()( tfttfttf⋅∆+=∆+)(')()( ttfttfttf∆+⋅∆+=∆+∆t = 1∆t = 0.1 Taylor Series:as ∆t  0, the higher order terms quickly approach 0, giving Euler’s approximation:...)(''')('')(')()(32+⋅∆+⋅∆+⋅∆+=∆+ tfttfttfttfttfgiving Euler’s approximation: The higher order terms are the error in the approximation )(')()( tfttfttf⋅∆+=∆+ In interactive simulations, we often do not have a function for the derivative, rather we have instantaneous values for the derivativeExample: driving simulationExample: driving simulationacceleration value is set by driver’s joystick Euler’s method can still be applied (twice)to compute velocity and position assumes that current value of acceleration was constant through previous time step Acceleration is set by user Function updates velocity, then positionvoid Vehicle::update(double elapsed_time){{double dt = elapsed_time/1000.0;velocity += acceleration*dt;double vx = velocity*sin(heading);double vy = velocity*cos(heading);double x = position.getX() + vx*dt;double y = position.getY() + vy*dt;position.moveTo(Point(x,y));}Numeric Integration Methods(explicit or forward) Euler Integration2ndorder Runga Kutta Integration (Midpoint Method)4thorder RungaKuttaIntegration4thorder RungaKuttaIntegrationThree slides borrowed from:www.cse.ohio-state.edu/~parent/BookWebMaterial/Chapter07/Powerpoint/Ch7_1_Integration.pptRunge Kutta Integration: 2ndorderAka Midpoint MethodFor unknown function, f(t); known f ’(t)ttftftfiii∆⋅′+=+)(21)()(21)(itf)(itf′ttftftfiii∆⋅′+=++)()()(211)(21+′itfRunge Kutta Integration: 4thorderFor unknown function, f(t); known f ’(t))(tf)(21+′itf)(itf′)(21+′itf1k2k3k)(21+itf)(itf2+i)(1+′itf++++=+ 4321161313161)()( kkkkhtftfii4k Euler: x(t+∆t) = x(t) + f'(t, x)*∆t Runga-Kutta 2: k1 = f'(t, x)*∆tk2 = f'(t+∆t /2, x+k1/2)*∆tBetter methods use more points on f’ to better approximate results.k2 = f'(t+∆t /2, x+k1/2)*∆t x(t+∆t) = x + k2 Runga-Kutta 4: k1 = f'(t, x)*∆t k2 = f'(t+∆t /2, x+k1/2)*∆t k3 = f'(t+∆t /2, x+k2/2)*∆t k4 = f'(t+∆t, x+k3)*∆t x(t+∆t) = x + k1/6 + k2/3 + k3/3 + k4/6approximate results. A hand-coded Euler integration of x’ = x + 3from pylab import *# initial value for xx = 1.0while t < 5.0:# one step of integrationdx = x + 3x = x + dx*dt# store for plottingx = 1.0# initial value for timet = 0.0# time step for integrationdt = 0.1x_data = []t_data = []# store for plottingx_data.append(x)t_data.append(t)# increment timet = t + dtxlabel("t = time")ylabel("df/dt")plot(t_data, x_data, "b")show()hand_coded_euler.py The attached Python class implements Euler, RK2 and RK4 techniques. It operates on a set of first order differential equations.Class constructor receives diffeqsand ∆tClass constructor receives diffeqsand ∆t Calls to integrate function move integration forward by ∆t timeintegrator.py# function to compute f'(t,x)def diff_eqs(t, x):# create return vectordx= [0.0]*num_eqs Solve x’ = x + 3dx= [0.0]*num_eqs# solve diff eqsdx[0] = x[0] + 3# return answerreturn dxfirst_example.py# setup the simulation # initial values for xx = [1.0] Solve x’ = x + 3x = [1.0]# initial value for timet = 1.0# time step for integrationdt = 0.1# set up the integrator integrator = Integrator(num_eqs, dt, diff_eqs)first_example.py# lists to store resultsx_data = []t_data = [] Solve x’ = x + 3# run the simulationwhile t < 5.0:# one step of integration for current timex = integrator.integrateEuler(x, t)# store current values for plottingx_data.append(x)t_data.append(t)# increment timet = t + dtfirst_example.py Result:first_example.pynum_eqs = 2k1 = 0.1 # prey birthk2 = 0.001 # predator/prey interaction x1' = k1*x1 - k2*x1*x2 # prey population change x2' = k3*x1*x2 - k4*x2 # predator population changek2 = 0.001 # predator/prey interactionk3 = 0.001 # predator/prey interactionk4 = 0.3 # predator death# define two differential equationsdef diff_eqs(t, x):dx = [0.0]*num_eqsdx[0] = k1*x[0] - k2*x[0]*x[1]dx[1] = k3*x[0]*x[1] - k4*x[1]return dxpredator_prey_euler.py# setup simulationx = [100.0, 100.0] # initial populationsdt = 0.1t = 0.0 x1' = k1*x1 - k2*x1*x2 # prey population change x2' = k3*x1*x2 - k4*x2 # predator population changet = 0.0integrator = Integrator(num_eqs, dt, diff_eqs)# run simulationwhile t < 1000:x = integrator.integrateEuler(x, t)p1_data.append(x[0])p2_data.append(x[1])t_data.append(t)t = t + dtone call solves both equationscurrent values of each equation are in arraypredator_prey_euler.py x1' = k1*x1 - k2*x1*x2 # prey population change x2' = k3*x1*x2 - k4*x2 # predator population changeblue = population 1blue = population 1red = population 2Euler integrationpredator_prey_euler.py x1' = k1*x1 - k2*x1*x2 # prey population change x2' = k3*x1*x2 - k4*x2 # predator population changeblue = population 1blue = population 1red = population 2Runga-Kutta 2 integrationx =