PowerPoint-presentatie

1
A Core Course on Modeling
Week 3 Time for Change

• Enter ACCEL
• basics
• running a script
• Time in ACCEL
• the operator
• examples

2
A Core Course on Modeling
Week 3 Time for Change
ACCEL

• intro and PR
• line-by-line, navigate a script
• assistance
• enter a script as text
• numerical analysis
• optimization
• larger view on simulation output

3
A Core Course on Modeling
Week 3 Time for Change
3
ACCEL
Example of a script
p5 qslider(10,0,20) rpq
paste into IO/edit box click run
Every line introduces a quantity Quantities can
be constants (p) Quantities can be functions r
f(p,q) Quantities can be user-entered (q)
4
A Core Course on Modeling
Week 3 Time for Change
4
ACCEL

• If something STRANGE happens
• don't panic
• goto IO/edit
• ctrl-A (select all)
• ctrl-C (copy all)
• ctrl-V into text editor to save your script
• reload ACCEL
• goto IO/edit
• ctrl-V to load script into ACCEL
• retry

image http//partlycloudyjuly.wordpress.com/2011/
02/25/look-both-ways/
5
A Core Course on Modeling
Week 3 Time for Change
5
ACCEL
Example of a script with comment
sslider(10,0,20) // this is a slider rpq
// this is an expression p5 // this is a
constant qst // this is an
expression tpow(s,3) // this is a standard
// function
to see values click show/hide values to see
dependencies click on pauze, next click on any
quantity
6
A Core Course on Modeling
Week 3 Time for Change
6
ACCEL
Example of a script with visual output
sslider(10,0,100) zplot(gr1,gr2) gr1str,s
gr2str,s 10 str'xmodeintp,ymodeshi
ft,ref1'
plot(graph1,graph2,, graphn) plots n graphs
graphi format,data format string, e.g. '
xmodeintp,ymodeshift,ref1' data one or
more quantities
Black Magic for now
7
A Core Course on Modeling
Week 3 Time for Change
7
time in ACCEL
Remember recursive functions Qcurrent f
(Qprev, Pprev) Simplest example timecurrent
timeprev 1
8
A Core Course on Modeling
Week 3 Time for Change
8
time in ACCEL
pt tt11
ACCEL uses 1 to access previous value ACCEL
uses n to access n-th previous value, ngt0 ACCEL
initializes after modification in script Start
conditions ACCEL initializes all quantities to
0, unless you use (value in case historic
value is lacking)
9
A Core Course on Modeling
Week 3 Time for Change
9
time in ACCEL
pk kif(tgt0,k15,0) tt11
ACCEL can be forced to re-initialize after nr
steps
10
A Core Course on Modeling
Week 3 Time for Change
10
time in ACCEL
zdescartes(gr1,gr2) gr1str1,s gr2str2,
50s-s1 s5025sin(t/10) zdescartes(gr1,g
r2) gr1locationsxmode'intp',ymode'shi
ft',values,edgesthickness2,col_r255
gr2locationsxmode'intp',ymode'shift',v
alue50s-s1,edgesthickness2,col_b255
tt11
Implement first derivative
image http//www.treklens.com/gallery/photo593123
.htm
11
A Core Course on Modeling
Week 3 Time for Change
11
time in ACCEL
zdescartes(gr1,gr2) gr1locationsxmode'i
ntp',ymode'shift',value5010rate,edgesth
ickness2,col_rif(rategt0,0,255),col_gif(rategt0,2
55,0) gr2locationsxmode'intp',ymode'
shift',value50supply,edgesthickness2,col_b
255 rateslider(0,-0.5,0.5) supplysupply1r
ate
Implement integral
image http//www.themarketingexpert.net/2011/12/l
eaky-bucket.html
12
A Core Course on Modeling
Week 3 Time for Change
12
time in ACCEL
zdescartes(gr2,gr1) gr1locationsxmode'i
ntp',ymode'shift',values,edgesthickness2
,col_r255 gr2locationsxmode'intp',ym
ode'shift',valueds,edgesthickness2,col_b25
5 sslider(50,0,100) ds(1-damp)sdampds1
dampslider(0.5,0.01,0.99)
Implement damping
image http//www.acecontrols.co.uk/product-range
13
A Core Course on Modeling
Week 3 Time for Change
13
time in ACCEL
zdescartes(gr1,gr2) gr1locationsxmode'i
ntp',ymode'shift',values,edgesthickness2
,col_r255 gr2locationsxmode'intp',ym
ode'shift',valuesdelay,edgesthickness2,co
l_b255 sslider(50,0,100) delayslider(1,1,10
0)
Implement delay
image http//www.rainbowresource.com/prodlist.php
?subject20category8753
14
A Core Course on Modeling
Week 3 Time for Change
14
time in ACCEL
t ? i? K(t) ? Ki u(t) ? ui u(t) ?
(u(t?)-u(t))/? (ui1-ui)/? v(t)
?(v(t?)-v(t))/? (vi1-vi)/? So ui1 ui?vi
vi1 vi?ai vi?Ki/m
vi? C(urest-ui) /m Recursive
function QcurrF(Qprev,Pprev) ucurrF1(uprev,vpre
v)uprev?vprev vcurrF2(vprev,uprev)
vprev? C(urest-uprev) /m, with
suitable u0 and v0
Equal intervals a mass-spring system with
damping. From physics, we know Kma where
KC(urest-u) and au. Write vu,
then av. Here, KK(t), uu(t). Let us
approximate by sampling. For ? ?0
out
15
A Core Course on Modeling
Week 3 Time for Change
15
time in ACCEL
t ? i? K(t) ? Ki u(t) ? ui u(t) ?
(u(t?)-u(t))/? (ui1-ui)/? v(t)
?(v(t?)-v(t))/? (vi1-vi)/? So ui1 ui?vi
vi1 vi?ai vi?Ki/m
vi? C(urest-ui) /m Recursive
function QcurrF(Qprev,Pprev) ucurrF1(uprev,vpre
v)uprev?vprev vcurrF2(vprev,uprev)
vprev? C(urest-uprev) /m, with
suitable u0 and v0
Equal intervals a mass-spring system with
damping. From physics, we know Kma where
KC(urest-u) and au. Write vu,
then av. Here, KK(t), uu(t). Let us
approximate by sampling. For ? ?0
out
-?v
( -?vi)
( -?vprev)
16
A Core Course on Modeling
Week 3 Time for Change
16
time in ACCEL
Equal intervals WARNING ? ? 0 is not the same
as '?small'. The substitutions for u(t) and
v(t) are approximations only. As we will see
in a later example, aliasing errors are worse
when the sampled signal changes rapidly, compared
to the sampling rate. Concrete unless ? is
much smaller than the period of the oscillation
(?m/C), the approximations are bad with the
risk of instability.
17
A Core Course on Modeling
Week 3 Time for Change
17
time in ACCEL
Infinitesimal intervals the mass-spring system
with damping revisited. From physics, we
know Kma where KC(urest-u) and au.
For some purposes the numerical solution
(sampling) is not acceptable. Interested in
outcome? ? use better numerical
methods Interested in insight? ? try to use
symbolic methods, i.e. analysing differential
equations Only possible in very few special
cases In particular linear (sets of) DEs.
out
18
A Core Course on Modeling
Week 3 Time for Change
18
time in ACCEL
Infinitesimal intervals the mass-spring system
with damping revisited. From physics, we
know Kma where KC(urest-u) and au.
Here, KK(t), uu(t). Let us try a solution of
the form uurestA sin(t/T) u
-AT-2sin(t/T) Substitute back -mAT-2sin(t/T)C(u
rest-urest-Asin(t/T) ), mT-2C, or T ?m/C (
T0)
out
19
A Core Course on Modeling
Week 3 Time for Change
19
time in ACCEL
Infinitesimal intervals the mass-spring system
with damping revisited. From physics, we
know Kma where KC(urest-u) and au.
Here, KK(t), uu(t). try uurestAe?t, then
C??m?20, or ?1,2(- ???(?2-4mC))/2m Let
?0?(4mC). For ? ?0? critical damping for ? lt
?0? oscillations TT0/?(1- (? /?0)2) for ? gt ?0?
super critical damping (further details lectures
dynamical systems)
out
-?v
By substituting u Ae?t into C(urest-u)-?umu
and using uA?e?t uA?2e?t. C(urest-urest-Ae?
t)- ?A?e?tmA?2e?t Must hold for any t, so
indeed C ?? m?20
20
A Core Course on Modeling
Week 3 Time for Change
20
time in ACCEL
A mass-spring system
cslider(5,1,20.0) deltaslider(0.1,0.001,0.2) ms
lider(5,1,10) muslider(0,0,1.5) resetbutton() z
descartes(gr1,gr2) gr1locationsxmode'intp
',ymode'shift',value4010u,edgesthicknes
s2,col_r255 gr2locationsxmode'intp',y
mode'shift',value4510u_exact,edgesthickne
ss2,col_g255 F-cu1-muv1 aF/m discrmu
mu-4mc lambda_imgcond(discrgt0,0,sqrt(-discr)/(2
m)) lambda_realcond(discrgt0,(-musqrt(discr))/(2
m),-mu/(2m)) u_exactcond(discrgt0,exp(lambda_rea
lt),exp(lambda_realt)cos(lambda_imgt)) tif(re
set,0,t1delta) uif(reset,1,u1deltav) vif(
reset,0,v1deltaa)
21
A Core Course on Modeling
Week 3 Time for Change
21
time in ACCEL
y(t) vy t gt2/2 x(t) vx t y(t) 0 , so t
2 vy / g Then x 2 vx vy / g ? 2
(cos ? sin ? ) / g sin 2 ?
/g Maximize x ? 45o But what in case of
damping?
y
?
x
http//thehumanmarvels.com/1634/zazel-the-human-ca
nnonball/talents
A cannon shot simulation
22
A Core Course on Modeling
Week 3 Time for Change
22
time in ACCEL
airDampslider(0.005,0,0.03) dtslider(0.1,0.02,40
.0) phislider(0,0,1.5) plotBalldescartes(locat
ionsfill'interior',xr.x,yr.y,rad2.5)pxr.
x g-0.003 groundDamp0.02 v0.55 boom(r1.vxlt0.
000001) newVxr1.vx(1-velDampdt) newVyr1.vy
(g-velDampr1.vy)dt newXmin(r1.xdtr1.vx
,100) newYmax(r1.ydtr1.vy,0) rif(boom,star
tR,simR) simR'x'newX,'y'newY,'vx'newVx,'vy'n
ewVy startR'x'0,'y'20,'vx'vcos(phi),'vy'v
sin(phi) velDampif(r1.ygt0,airDamp,groundDamp)
A cannon shot simulation
23
A Core Course on Modeling
Week 3 Time for Change
23
time in ACCEL
Lanchester's Law e1 nr troops of army 1 e2 nr
troops of army 2 ?1 relative strength of 1 ?2
relative strength of 2 then de1/dt - ?2 e2
de2/dt - ?1e1 Given e1(t0)e10,
e2(t0)e20, who wins?
http//www.napolun.com/mirror/napoleonistyka.atspa
ce.com/default.htm
A combat simulation
24
A Core Course on Modeling
Week 3 Time for Change
24
time in ACCEL
pcheck(true) e1if(p,e1_0,max(0,e11-lambda_2e2
1)) e2if(p,e2_0,max(0,e21-lambda_1e11)) la
mbda_1slider(0.01,0,0.02) lambda_2slider(0.01,0,
0.02) e1_0slider(100,0,500) e2_0slider(100,0,500
) resdescartes(graph1,graph2) graph1locations
xmode'intp',ymode'shift',valuee1/5,rad
1,fcol_r255,fill'interior' graph2locations
xmode'intp',ymode'shift',valuee2/5,rad1,
fcol_g255,fill'interior'
A combat simulation
25
A Core Course on Modeling
Week 3 Time for Change
25
time in ACCEL
Lotka Volterra model br birth rate of rabit
(autonomous) bf birth rate of fox, feeding on
rabit dr rabits killed by fox df death rate
of fox (autonomous) ?r / ? t r (br dr f) ?f /
? t f (bf r df )
http//q-rai.deviantart.com/
A predator - prey simulation
26
A Core Course on Modeling
Week 3 Time for Change
26
time in ACCEL
bfslider(1.5,0,5) brslider(3,0,5) dfslider(1,0,
5) drslider(1,0,5) dtslider(0.15,0.01,0.2) reset
button() pdescartes(gboth,gfox,grab) gbothl
ocationsfill'interior',nrLocations50,rad2.5,f
col_r50,xmode'shift',value103f,ymode'sh
ift',value103r gfoxlocationsymode'in
tp',xmode'shift',103f,edgesthickness2,c
ol_r255 grablocationsxmode'intp',ymod
e'shift',103r,edgesthickness2,col_b255
fif((timelt2)reset,150,max(0.01,f1dt(bff1
r1-dff1))) rif((timelt2)reset,150,max(0.0
1,r1dt(brr1-drf1r1))) timetime11
recursive functions
A predator - prey simulation
27
A Core Course on Modeling
Week 3 Time for Change
27
time in ACCEL
v'w
v'r
http//joyreactor.com/tag/billiards
? ?1
?r
?w - ?r ? ?? ? (vw-vr , ?)
vr
vw
A billiards simulation
28
A Core Course on Modeling
Week 3 Time for Change
28
time in ACCEL
v'w
v'r
? ?1
?r
?w - ?r ? ?? ? (vw-vr , ?)
vr
vw
Proof Energy conservation vr2vw2(vr?)2(vw-?)
2, hence vr2vw2(vr??)2(vw-??)2,
hence 0(2?(vr, ?)?2)(-2?(vw, ?)?2), hence ?
(vw-vr, ?), QED
A billiards simulation
29
A Core Course on Modeling
Week 3 Time for Change
29
time in ACCEL
v'w
v'r
? ?1
?r
?w - ?r ? ?? ? (vw-vr , ?)
vr
vw
Proof Energy conservation vr2vw2(vr?)2(vw-?)
2, hence vr2vw2(vr??)2(vw-??)2,
hence 0(2?(vr, ?)?2)(-2?(vw, ?)?2), hence ?
(vw-vr, ?), QED
v'rvr(vw-vr , ?) ? v'wvw(vr-vw , ?) ?
A billiards simulation
30
A Core Course on Modeling
Week 3 Time for Change
30
time in ACCEL
and a further handfull of details
maxX40 minX-40 maxY25 minY-25 halfWay50 //
half of the screen used to correct the queueX
and queueY coordinates helpcheck(false) // plot
the direction of the queue nBalls3 // nr
balls // This is also hard coded in the colors
if the nr balls should change, only the color //
arrays have to be updated rhoBallslider(1.4,0.6,4
) // radius of the balls mBall5 // mass of the
balls timetime11 rollDamp0.994 // the
damping coefficient for the rolling
balls collDamp0.6 // the damping coefficient for
colliding balls eThreshold0.5 // the minimum
kinetic energy for the balls to stay
rolling dT0.8 // to go from physical time to
simulation time PRplot(plotTable,plotShadow,plot
Balls,plotHighLight,plotQueue,plotQueue2) //
one graph consisting of three balls. The x and y
coordinates are // taken from the 'x'- and
'y'-components of the r-vector that has to be
transposed in order to get // the x-s and y-s
in two separate vectors. tVec(i,j)_at_(r1,j)-_at_(r
1,i) // the vector pointing from ball i to
ball j, with coefficients x and
y cpl(i,vSequence(0,nBalls),cplOneBall(i),vAppen
d) // the couplings matrix containing all
info about the relations to ball i, // that
is 'close' to indicate if these two balls are in
a collision-state // if close is true,
'force' is the current reaction force between
them. tMat(i,vSequence(0,nBalls),tMatOneBall(i),
vAppend) // the tMatrix contains, for every
pair, i-j, the t-vector between the
centres tMatOneBall(i)(j,vSequence(0,nBalls),tVe
c(i,j),vAppend) // the tMatOneBall vector
contains, for every ball, // the the
t-vectors between its centre and the other
centres touch(i,j)((vNormEuclid(_at_(_at_(tMat,i),j))lt(
_at_(rho,i)_at_(rho,j))) (i!j)
vDot(_at_(_at_(tMat,i),j),_at_(p1,i)/_at_(m,i)-_at_(p1,j)/_at_(m
,j))gt0) // condition for colliding contact.
Three terms // 1. Are the balls close
enough? // 2. No self-collision? // 3. Is
the relative velocity opposite to the vector
connecting the centres? cplOneBall(i)(j,vSequenc
e(0,nBalls),cond(touch(i,j),'close'true,'force'
(mpt(i,j)/m1m2tt(i,j))_at_(_at_(tMat,i),j),'close'fa
lse,'force'0),vAppend) // the cplOneBall
vector contains, for ball i, the info for the
collisions between this // ball and all other
balls. // It sets the value to the boolean
'close', and if closetrue, the // current
force vector. The derivation of the force vector
// is based on conservation of momentum,
angular momentum and kinetic energy // in a
coordinate-free version. mpt(i,j)-2vDot(_at_(m,j)_at_
(p1,i)-_at_(m,i)_at_(p1,j),_at_(_at_(tMat,i),j)) //
the product (((m1p2-m2p1)t),t) the numerator of
the force vector m1m2tt(i,j)(_at_(m,i)_at_(m,j))vDot(
_at_(_at_(tMat,i),j),_at_(_at_(tMat,i),j)) // the product
(m1m2)(t,t) the denominator of the force
vector f(i,vSequence(0,nBalls),forceOnOneBall(i)
,vAppend) // calculate the forces for all
balls by concatenating forceOnOneBall(i)(j,vSequ
ence(0,nBalls),cond(_at_(_at_(_at_(cpl,i),j),'close'),_at_(_at_(_at_
(cpl,i),j),'force'),0),add) // adding the
forces due to all other of balls
rho(i,vSequence(0,nBalls),rhoBall,vAppend)
// the radii of the balls
m(i,vSequence(0,nBalls),mBall,vAppend) //
the mass of the balls col'red'255,255,255,'gr
n'0,255,255,'blu'0,255,0 // colors of the
balls first one is red, two is white, three is
yellowish rif(gameState'roll',r1dTp/m,cond(
time1,(i,vSequence(0,nBalls),'x'halfWay(rand
om()-0.5),'y'0.5halfWay(random()-0.5),vAppend)
,r1)) // r is obtained by integrating v
at the starting time (1) give the initial
random positions pif(gameState'roll',rollDamp
wallCollide(p1f),'x'0,'y'0,'x'_at_(_at_(r1,1
),'x')-queueX,'y'_at_(_at_(r1,1),'y')-queueY,'x'0,
'y'0) // p is obtained by tame integrating
f at the starting time give the initial momenta.
The only // non-vanishing initial momentum is
the momentum of ball 1 this is obtained from the
queue position minus the centre of ball
1. wallCollide(a)wallCollideLeft(wallCollideRight
(wallCollideBottom(wallCollideTop(a)))) //
the collisions with all walls wallCollideLeft(a)
(i,vDom(a),cond((_at_(_at_(r1,i),'x')gt(minX_at_(rho,i)))
(_at_(_at_(a,i),'x')gt0),_at_(a,i),'x'-collDamp_at_(_at_(a,
i),'x'),'y'collDamp_at_(_at_(a,i),'y')),vAppend) //
collide left wall wallCollideRight(a)(i,vDom(a),
cond((_at_(_at_(r1,i),'x')lt(maxX-_at_(rho,i)))
(_at_(_at_(a,i),'x')lt0),_at_(a,i),'x'-collDamp_at_(_at_(a,i),'
x'),'y'collDamp_at_(_at_(a,i),'y')),vAppend) //
collide right wall wallCollideBottom(a)(i,vDom(a
),cond((_at_(_at_(r1,i),'y')gt(minY_at_(rho,i)))
(_at_(_at_(a,i),'y')gt0),_at_(a,i),'x'collDamp_at_(_at_(a,i),'x
'),'y'-collDamp_at_(_at_(a,i),'y')),vAppend) //
collide bottom wall wallCollideTop(a)(i,vDom(a),
cond((_at_(_at_(r1,i),'y')lt(maxY-_at_(rho,i)))
(_at_(_at_(a,i),'y')lt0),_at_(a,i),'x'collDamp_at_(_at_(a,i),'x
'),'y'-collDamp_at_(_at_(a,i),'y')),vAppend) //
collide top wall eKin(i,vSequence(0,nBalls),vDot
(_at_(p,i),_at_(p,i))/(2_at_(m,i)),add) // the kinetic
energy gameStatecond((eKin1lteThreshold)
!queueHit,'hit',cond(queueHit,'roll','hit')) //
determine which state we are in queueHitcursorB()
cXcursorX() cYcursorY() queueXcX-halfWay queu
eYcY-halfWay // the queue end position is
derived from the location of the
cursor plotQueuecond(gameState'hit','plotType
line,col_bvalue0,col_gvalue0,col_rvalue
0,xmodedata,ref1,ymodedata,ref2',halfWa
yqueueX,_at_(_at_(r1,1),'x'),halfWayqueueY,_at_(_at_(r
1,1),'y'),'plotTypeline,xmodedata,ref1,y
modedata,ref2',0,0,0,0) plotQueue2cond(
gameState'hit' help,'plotTypeline,col_gva
lue100,col_rvalue0,widthvalue0.4,xmode
data,ref1,ymodedata,ref2',halfWay10(_at_(_at_
(r1,1),'x')-queueX)queueX,_at_(_at_(r1,1),'x'),hal
fWay10(_at_(_at_(r1,1),'y')-queueY)queueY,_at_(_at_(r1
,1),'y'),'plotTypeline,xmodedata,ref1,y
modedata,ref2',0,0,0,0) // the queue is
drawn when in hit-mode plotTable'plotTypevbar,y
Basevalue25,widthvalue80,heightvalue50
,col_rvalue0', // plot the
table plotBalls'plotTypebubble,xmodedata,ref
1,ymodedata,ref2,diametermodedata,ref3
,col_rmodedata,ref4,col_gmodedata,ref5,c
ol_bmodedata,ref6',halfWay_at_(vTranspose(r),'x
'),halfWay_at_(vTranspose(r),'y'),2rho,_at_(col,'red')
,_at_(col,'grn'),_at_(col,'blu') // plot the
balls plotShadow'plotTypebubble,xmodedata,re
f1,ymodedata,ref2,diametervalue'2rhoBa
ll',col_rvalue0,col_gvalue0,col_bvalue
0,col_avalue0.5',halfWay-0.15rhoBall_at_(vTra
nspose(r),'x'),halfWay-0.5rhoBall_at_(vTranspose(r)
,'y') // plot the shadows beneath the
balls plotHighLight'plotTypebubble,xmodedata
,ref1,ymodedata,ref2,diametervalue'0.75
rhoBall',col_rvalue255,col_gvalue255,co
l_bvalue255',halfWay0.1rhoBall_at_(vTranspose(
r),'x'),halfWay0.5rhoBall_at_(vTranspose(r),'y')
// plot the highlights on the balls
A billiards simulation