globals [x-max y-max diam filename
g b k dt xzero yzero]
turtles-own
[
x y vx vy ax ay xtemp ytemp vxtemp vytemp
kx1 kx2 kx3 ky1 ky2 ky3 jx1 jx2 jx3 jy1 jy2 jy3
]
to setup ; Setup the model for a run, build a plot.
;; (for this model to work with NetLogo's new plotting features,
;; __clear-all-and-reset-ticks should be replaced with clear-all at
;; the beginning of your setup procedure and reset-ticks at the end
;; of the procedure.) ..
__clear-all-and-reset-ticks
setup-globals
setup-patches
setup-turtles
setup-slinky
set xzero 0
set yzero y-max - 3
ask turtle 0 [ setxy 0 yzero ]
setup-plot
end
to go
set k spring-force
set b mydamping
set g mygravity
;; Now go use Runge-Kutta
ifelse (do-drop)
[ drop-slinky ]
[ ifelse (drive)
[driven ]
[fixed-top]
]
do-plots
if filename = 0 [setup] ; an attempt to work even tho user forgets setup
; this allows us to grab a node and drag it . . .
let t 0
if mouse-down? [
set t closest-xy mouse-xcor mouse-ycor turtles ;
while [mouse-down?] [
ask t [
set x mouse-xcor
set y mouse-ycor
setxy x y
]
]
]
check-movie
; if (ticks = 1000) [stop] ; stop after a while, if you want . . .
tick
end
to setup-globals ; separate procedure for setting globals
set diam 1.5
set x-max max-pxcor - (diam / 2) + 1 ; 0.5
set y-max max-pycor - (diam / 2) + 1 ; 0.5
set filename "" ; change to collect images (or just use command center after setup)
set-default-shape turtles "circle"
set dt 0.01 ; step size for Runge Kutta
end
to setup-patches
; make the ground green
ask patches [ if pycor < (5 + min-pycor) [set pcolor green ]]
end
to setup-turtles
crt number-of-nodes [
ifelse (who = 0)
[set color blue]
[ifelse (who = number-of-nodes - 1)
[set color magenta]
[set color yellow]
]
set label-color black
set size diam
; set label who
set vx 0
set vy 0
set ax 0
set ay 0
set x 0
set y (y-max - who - 5) ; start them out stacked . . .
setxy x y
]
end
to setup-slinky ; Build the slinky
ask turtles with [(who > 0)] [
let mywho who
ask self [
connect-edge turtle (mywho - 1)
]
]
end
; find new velocity and position of nodes using 4th order Runge Kutta.
to update-coordinates
let kx4 vx
let ky4 vy
;
; here is the first place where we put in the differential equation
; we are solving
; ax, ay is the acceleration from the springs (neighbors),
; - b * vx, -b * vy is damping (friction) -- a simple model
; where damping is linear in velocity
; -g is the force of gravity (constant, - y direction . . .)
;
let jx4 (ax - b * vx )
let jy4 (ay - b * vy - g)
set vx vxtemp + (dt / 6)*(jx1 + (2 * jx2) + (2 * jx3) + jx4)
set vy vytemp + (dt / 6)*(jy1 + (2 * jy2) + (2 * jy3) + jy4)
set x xtemp + (dt / 6)*(kx1 + (2 * kx2) + (2 * kx3) + kx4)
set y ytemp + (dt / 6)*(ky1 + (2 * ky2) + (2 * ky3) + ky4)
set ax 0
set ay 0
; If a node found below the green ground we give it positive velocity.
; This corresponds to a perfectly elastic collision with the ground.
;
if y < 5 + min-pycor [ set vy abs vy ]
;; Here we hide nodes that have moved off screen. They are still part of the simulation, and
;; if they return on screen they become visible again.
ifelse (y > max-pycor) or (y < min-pycor)
[ set hidden? true
setxy x max-pycor]
[set hidden? false
setxy x y ]
end
to oscillate
;
; this is for a driven system . . . either the top (blue) node
; moves in a simple vertical oscillation,
; or the top node moves in a circle
;
ifelse (two-d-drive)
[set x xzero + drive-amplitude * cos (ticks * drive-rate / ( 2 * pi))]
[set x xzero]
set y yzero + drive-amplitude * sin (ticks * drive-rate / ( 2 * pi))
setxy x y
end
to fixed-top
; this is the Runge Kutta solver with the top node fixed
;
ask turtles with [who != 0] [find-acceleration ]
ask turtles with [who != 0][find-k1]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][find-k2]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][find-k3]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][update-coordinates]
end
to driven
; this is the Runge Kutta solver with driven top node
;
ask turtle 0 [ oscillate ]
ask turtles with [who != 0] [find-acceleration ]
ask turtles with [who != 0][find-k1]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][find-k2]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][find-k3]
ask turtles with [who != 0][find-acceleration ]
ask turtles with [who != 0][update-coordinates]
end
to drop-slinky
; this is the Runge Kutta solver with all nodes
; free to move
;
ask turtles [find-acceleration ]
ask turtles [find-k1]
ask turtles [find-acceleration ]
ask turtles [find-k2]
ask turtles [find-acceleration ]
ask turtles [find-k3]
ask turtles [find-acceleration ]
ask turtles [update-coordinates]
end
; The following three procedures are steps in the Runge Kutta process.
; note that our differential equation shows up in each of these . . .
;
to find-k1
set xtemp x
set ytemp y
set vxtemp vx
set vytemp vy
set kx1 vx
set ky1 vy
set jx1 (ax - b * vx )
set jy1 (ay - b * vy - g)
set vx vxtemp + jx1 * (dt / 2)
set vy vytemp + jy1 * (dt / 2)
set x xtemp + kx1 * (dt / 2)
set y ytemp + ky1 * (dt / 2)
set ax 0
set ay 0
end
to find-k2
set kx2 vx
set ky2 vy
set jx2 (ax - b * vx )
set jy2 (ay - b * vy - g)
set vx vxtemp + jx2 * (dt / 2)
set vy vytemp + jy2 * (dt / 2)
set x xtemp + kx2 * (dt / 2)
set y ytemp + ky2 * (dt / 2)
set ax 0
set ay 0
end
to find-k3
set kx3 vx
set ky3 vy
set jx3 (ax - b * vx )
set jy3 (ay - b * vy - g)
set vx vxtemp + jx3 * dt
set vy vytemp + jy3 * dt
set x xtemp + kx3 * dt
set y ytemp + ky3 * dt
set ax 0
set ay 0
end
; Using Newton's second law F=ma, find the acceleration of each node
; based on the interaction with link neighbors.
; The mass is set to 1.
; Note that when spring-length = 0, this is a simple spring with force linearly
; dependent on length.
; With positive spring-length, the spring pulls (or pushes) toward
; a rest length of spring-length . . .
;
to find-acceleration
ask my-links [
let fx 0
let fy 0
;let f (-1) * k * (link-length)
let f k * (spring-length - link-length)
let myheading link-heading
ifelse (myself = [end2] of self)
[set fx f * sin myheading
set fy f * cos myheading
]
[set fx (-1) * f * sin myheading
set fy (-1) * f * cos myheading
]
ask myself [set ax (ax + fx)]
ask myself [set ay (ay + fy)]
]
end
to setup-plot
end
;; plot y locations
to do-plots
set-current-plot-pen "top"
let ytop 0
ask turtle 0 [ set ytop y ]
plot ytop
set-current-plot-pen "bottom"
ask turtle (number-of-nodes - 1) [ set ytop y ]
plot ytop
set-current-plot-pen "c-of-mass"
plot (sum [y] of turtles) / number-of-nodes
end
to check-movie ; if filename is non-empty, make another snapshot
if length filename > 0 [
if (1 = (ticks mod 10))
[ export-view word word word filename (round log ticks 10) ticks ".png" ]
]
end
;;;; Edge & Node utilities
to connect-edge [other-node] ; node proc: attach an edge between self and other
ask other-node [ create-link-with myself
[
set color gray
]
]
end
;;;; Utilities
to-report closest-xy [myx myy agent-set] ; Return closest agent to x, y
report min-one-of agent-set [distancexy-nowrap myx myy]
end
@#$#@#$#@
GRAPHICS-WINDOW
204
10
583
410
20
-1
9.0
1
10
1
1
1
0
1
0
1
-20
20
0
40
0
0
1
ticks
30.0
BUTTON
7
10
62
43
NIL
setup
NIL
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
BUTTON
65
10
120
43
NIL
go
T
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
BUTTON
124
10
179
43
step
go
NIL
1
T
OBSERVER
NIL
NIL
NIL
NIL
1
SLIDER
7
94
185
127
number-of-nodes
number-of-nodes
1
20
10
1
1
NIL
HORIZONTAL
SLIDER
6
150
182
183
spring-force
spring-force
0
30
15
5
1
NIL
HORIZONTAL
SLIDER
8
193
184
226
spring-length
spring-length
0
10
1.75
0.25
1
NIL
HORIZONTAL
PLOT
596
35
796
227
mass-plot
ticks
y pos
0.0
10.0
0.0
10.0
true
false
"" ""
PENS
"top" 1.0 0 -13345367 true "" ""
"bottom" 1.0 0 -5825686 true "" ""
"c-of-mass" 1.0 0 -1184463 true "" ""
SWITCH
7
52
97
85
do-drop
do-drop
1
1
-1000
SLIDER
8
262
180
295
mygravity
mygravity
0
10
2.4
0.2
1
NIL
HORIZONTAL
SLIDER
9
304
181
337
mydamping
mydamping
0
2
0.3
0.1
1
NIL
HORIZONTAL
SWITCH
101
53
191
86
drive
drive
1
1
-1000
SLIDER
598
240
770
273
drive-rate
drive-rate
0
20
10
1
1
NIL
HORIZONTAL
SLIDER
598
282
770
315
drive-amplitude
drive-amplitude
0
3
2
0.1
1
NIL
HORIZONTAL
SWITCH
601
335
712
368
two-d-drive
two-d-drive
1
1
-1000
@#$#@#$#@
## WHAT IS IT?
This is a dynamic model of a "slinky". The "slinky" is a chain of nodes connected by springs, and pulled down by gravity. There is also "damping" in the system (a "friction" force, opposite to (and linear in) velocity).
This model was partly inspired by this video: http://youtu.be/wGIZKETKKdw .
(but also by the desire to explore dynamical systems and numerical solvers . . . )
Some aspects of this model (especially the basic form of the Runge Kutta solver) came from the Bouncing Ball netlogo model: http://ccl.northwestern.edu/netlogo/models/community/Bouncing_Ball
## HOW IT WORKS
This model uses a 4th order Runge Kutta numerical solver for movement of the nodes. There are several modes:
Default: One node (the top) is fixed, the others are free to move. The moving nodes are pulled down by a fixed gravitational force. They also experience forces caused by "springs" attaching them to other nodes. These "springs" have a rest-length, and exert a (linear) force on their ends, trying to return to their rest-length.
Do-drop: This releases the top node, so they all fall (if gravity is on).
Drive: In this mode, the top node is "driven" -- either oscillating up and down, or rotating in a circle.
There is also a "ground" at the bottom. If nodes hit the ground, they bounce back up (elastic collision). If nodes go off screen, they become invisible, but stay active in the model. The model wraps horizontally, but not vertically . . .
## HOW TO USE IT
Click on setup to initialize the slinky. Click on go to start the model. While the model is running, you can change the various parameters. You also can grab a node with the mouse and move it around.
The plot shows the y coordinates of the top node, the bottom node, and the center-of-mass of the system.
You can slow things down to see better what is happening by using the Interface "speed" slider.
## THINGS TO NOTICE
This is a fairly simple model.
## THINGS TO TRY
Grab nodes with the mouse while the model is running. Change the various parameter settings while the model is running.
Movies: (Use this with local models, not via the web) You can make a sequence of images for creating a movie of the model. To do this, type 'set filename "movie"' in the command center immediately after clicking "setup". When the model sees the filename has been set to a non-null value while running ("go"), it will take a picture every step of the model. The filenames will look like movie0001.png movie0002.png ... etc. You can stop the process by clicking on "go" when you have enough images. Alternatively, the model can be stopped by "set stoptick 100" or some similar value, again using the command center.
## EXTENDING THE MODEL
Go for it! :-)
I'm hoping to add a "pendulum" mode, but not yet . . .
## NETLOGO FEATURES
Nothing particularly special -- uses "links" . . .
## RELATED MODELS
As noted above, the "Bouncing Ball" http://ccl.northwestern.edu/netlogo/models/community/Bouncing_Ball
## CREDITS AND REFERENCES
There are a variety of (somewhat) related physics applets available here: http://www.myphysicslab.com/
@#$#@#$#@
default
true
0
Polygon -7500403 true true 150 5 40 250 150 205 260 250
circle
false
0
Circle -7500403 true true 35 35 230
line
true
0
Line -7500403 true 150 0 150 301
@#$#@#$#@
NetLogo 5.0.3
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
@#$#@#$#@
default
0.0
-0.2 0 1.0 0.0
0.0 1 1.0 0.0
0.2 0 1.0 0.0
link direction
true
0
Line -7500403 true 150 150 90 180
Line -7500403 true 150 150 210 180
@#$#@#$#@
0
@#$#@#$#@