Programming Solutions for Integration

1
Programming Solutions for Integration

• COMP155 / EMGT155
• Dec 1 2008

2
Integrals compute area
Image from http//en.wikipedia.org/wiki/Integral
3
Calculus in Simulation

• Well 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
  position
• Example We know acceleration, need to
  compute velocity and position

4
Time Slicing

• Simulation is updated at regular interval in time
• dt time step/slice
• To approximate integration
• compute derivative at time steps
• make assumption about derivative value at other
  times
• accumulate area to approximate integration

5
Time slice approximation

• ?t 1, assume current value held for
  previous time slice

blue correct derivative functionred
approximated derivative function
over-estimate of positive area
under-estimate of positive area
under-estimate of negative area
over-estimate of negative area
6
Time slice approximation
?t 1 blue ? actual velocity red ?
approx velocity green ? actual
position magenta ? computed position
7
Smaller ?t ? better approximation
?t 0.5
?t 0.01
euler_pos_vel.py
8
Smaller ?t ? better approximation
Definition of a derivative As ?t gets
smaller, we get closer to the true value of the
derivative.
9
Time slice approximation
Definition of derivative Forward Euler
Method Backward Euler Method (were using
backward Euler)
10
Eulers Method

• Taylor Seriesas ?t ? 0, the higher order
  terms quickly approach 0, giving Eulers
  approximation
• The higher order terms are the error in the
  approximation

11
Instantaneous Derivative Values

• In interactive simulations, we often do not
  have a function for the derivative, rather we
  have instantaneous values for the derivative
• Example driving simulation acceleration
  value is set by drivers joystick
• Eulers method can still be applied (twice)
  to compute velocity and position
• assumes that current value of acceleration was
  constant through previous time step

12
Example Vehicle Dynamics

• Acceleration is set by user
• Function updates velocity, then position

void Vehicleupdate(double elapsed_time)
double dt elapsed_time/1000.0 velocity
accelerationdt double vx velocitysin(headin
g) double vy velocitycos(heading) double
x position.getX() vxdt double y
position.getY() vydt position.moveTo(Point(x
,y))
13
Numeric Integration Methods
(explicit or forward) Euler Integration 2nd order
Runga Kutta Integration (Midpoint Method) 4th
order Runga Kutta Integration
Three slides borrowed from www.cse.ohio-state.edu
/parent/BookWebMaterial/Chapter07/Powerpoint/Ch7_
1_Integration.ppt
14
Runge Kutta Integration 2nd order Aka Midpoint
Method
For unknown function, f(t) known f (t)
15
Runge Kutta Integration 4th order
For unknown function, f(t) known f (t)
16
Numerical Integration Techniques

• Euler
• x(t?t) x(t) f'(t, x)?t
• Runga-Kutta 2
• k1 f'(t, x)?t
• k2 f'(t?t /2, xk1/2)?t
• x(t?t) x k2
• Runga-Kutta 4
• k1 f'(t, x)?t
• k2 f'(t?t /2, xk1/2)?t
• k3 f'(t?t /2, xk2/2)?t
• k4 f'(t?t, xk3)?t
• x(t?t) x k1/6 k2/3 k3/3 k4/6

Better methods use more points on f to better
approximate results.
17
Implementation - Python

• A hand-coded Euler integration of x x 3

while t lt 5.0 one step of integration
dx x 3 x x dxdt store for
plotting x_data.append(x)
t_data.append(t) increment time t t
dt xlabel("t time") ylabel("df/dt") plot(t_data
, x_data, "b") show()
from pylab import initial value for x x
1.0 initial value for time t 0.0 time step
for integration dt 0.1 x_data t_data
hand_coded_euler.py
18
Implementation - Python

• The attached Python class implements Euler,
  RK2 and RK4 techniques.
• It operates on a set of first order
  differential equations.
• Class constructor receives diffeqs and ?t
• Calls to integrate function move integration
  forward by ?t time

integrator.py
19
Example one diffeq

• Solve x x 3

function to compute f'(t,x) def diff_eqs(t,
x) create return vector dx
0.0num_eqs solve diff eqs dx0
x0 3 return answer return dx
first_example.py
20
Example one diffeq

• Solve x x 3

setup the simulation initial values for x x
1.0 initial value for time t 1.0 time
step for integration dt 0.1 set up the
integrator integrator Integrator(num_eqs, dt,
diff_eqs)
first_example.py
21
Example one diffeq

• Solve x x 3

lists to store results x_data t_data
run the simulation while t lt 5.0 one
step of integration for current time x
integrator.integrateEuler(x, t) store
current values for plotting x_data.append(x)
t_data.append(t) increment time t
t dt
first_example.py
22
Example one diffeq

• Result

first_example.py
23
Example two diffeqs

• x1' k1x1 - k2x1x2 prey population change
• x2' k3x1x2 - k4x2 predator population
  change

num_eqs 2 k1 0.1 prey birth k2 0.001
predator/prey interaction k3 0.001
predator/prey interaction k4 0.3 predator
death define two differential equations def
diff_eqs(t, x) dx 0.0num_eqs dx0
k1x0 - k2x0x1 dx1 k3x0x1
- k4x1 return dx
predator_prey_euler.py
24
Example two diffeqs

• x1' k1x1 - k2x1x2 prey population change
• x2' k3x1x2 - k4x2 predator population
  change

setup simulation x 100.0, 100.0 initial
populations dt 0.1 t 0.0 integrator
Integrator(num_eqs, dt, diff_eqs) run
simulation while t lt 1000 x
integrator.integrateEuler(x, t)
p1_data.append(x0) p2_data.append(x1)
t_data.append(t) t t dt
one call solves both equations
current values of each equation are in array
predator_prey_euler.py
25
Example two diffeqs

• x1' k1x1 - k2x1x2 prey population change
• x2' k3x1x2 - k4x2 predator population
  change

blue population 1 red population 2 Euler
integration
predator_prey_euler.py
26
Example two diffeqs

• x1' k1x1 - k2x1x2 prey population change
• x2' k3x1x2 - k4x2 predator population
  change

blue population 1 red population
2 Runga-Kutta 2 integration
x integrator.integrateRK2(x, t)
predator_prey_rk2.py
27
Example two diffeqs

• x1' k1x1 - k2x1x2 prey population change
• x2' k3x1x2 - k4x2 predator population
  change

Euler
RK2
Common problem with Euleraccumulating errors
RK4
28
Aside Chaotic Attractors

• Lorenz system
• x1 a(x2-x1)x2 (1d-x3)x1 x2x3
  x1x2 bx3
• blue ? Eulerred ? RK4

lorenz.py
29
Aside Chaotic Attractors

• Silnikov system
• x3 -ax3 x2 bx1(1 cx1 dx12)x2
  x3 x1 x2
• blue ? Eulerred ? RK4

silnikov.py
30
Higher-order Systems

• Our solution solves systems of first order
  equations
• To deal with higher-order systems, we can
  make substitutions to convert to first-order
  equations

31
Higher-order Systems

• Example projectile y(0) 0, y(0) 10,
  y(t) -9.8 x(0) 0, x(0) 10, x(t)
  0.0

y
y gravity -9.8
y 10
x 10
x
32
Higher-order Systems

• Example projectile
• y(0) 0, y(0) 10, y(t) -9.8
• x(0) 0, x(0) 10, x(t) 0.0
• Too get a system of 1st order equations,
  introduce new variables
• A0 yA1 yA2 xA3 x

• Derivatives of new variables
• A0 y A1A1 y -9.8A2 x A3A3
  x 0.0

33
Higher-order Systems

• projectile

function to compute f'(t,x) def diff_eqs(t,
A) create return vector dA
0.0num_eqs solve diff eqs dA0
A1 dA1 g dA2 A3 dA3
0 return answer return dA
system constants g
-9.8 meters/second m 2 kilograms y0
100 meters number
of equations num_eqs 4
projectile.py
34
Higher-order Systems

• projectile result

projectile.py
35
Harmonic Oscillation

• A spring and block system (no friction) will
  exhibit harmonic oscillation
  mx -kx

36
Harmonic Oscillation

• mx -kx x (-k/m)x
• New variables
• A0 x
• A1 x
• A0 A1
• A1 (-k/m)x (-k/m)A0
• initial values
• A0(0) 5
• A1(0) 0

spring_harmonic.py
37
Harmonic Oscillation

• Pendulum md? mg sin(?) b? k?
  ? (g/d) sin(?) (b/md)? (k/md)?
• m mass of bob g gravity d length of
  string b damping (friction) k spring
  (string)

38
Pendulum

• Pendulum md? mg sin(?) b? k?
  ? (g/d) sin(?) (b/md)? (k/md)?
• new variables A0 ? A1 ? A0
  ? A1 A1 ? (g/d) sin(A0)
  (b/md)A1 (k/md)A0

pendulum.py
39
Pendulum

• Pendulum, initial values A0(0) 1.5
  radians A1(0) 0.0

pendulum.py
40
Pendulum
angle vs. time
positions vs. time blue ? x position green ? y
position
pendulum.py