Transcript 11 Music of the Spheres

Slides for Chapter 11
Note to Instructors
This Keynote document contains the slides for “The Music of the Spheres”, Chapter 11 of
Explorations in Computing: An Introduction to Computer Science.
The book invites students to explore ideas in computer science through interactive tutorials where
they type expressions in Ruby and immediately see the results, either in a terminal window or a 2D
graphics window.
Instructors are strongly encouraged to have a Ruby session running concurrently with Keynote
in order to give live demonstrations of the Ruby code shown on the slides.
License
The slides in this Keynote document are based on copyrighted material from Explorations in
Computing: An Introduction to Computer Science, by John S. Conery.
These slides are provided free of charge to instructors who are using the textbook for their courses.
Instructors may alter the slides for use in their own courses, including but not limited to: adding new
slides, altering the wording or images found on these slides, or deleting slides.
Instructors may distribute printed copies of the slides, either in hard copy or as electronic copies in
PDF form, provided the copyright notice below is reproduced on the first slide.
© 2012 John S. Conery
The Music of the Spheres
Computer simulation and the N-body problem
✦
Running Around in Circles
✦
The Force of Gravity
✦
Force Vectors
✦
N-Body Simulation of the
Solar System
Explorations in Computing
© 2012 John S. Conery
Heliocentric Orbits
✦
Most ancient astronomers and
philosophers thought the planets
revolved around the Earth
❖
✦
Nicolaus Copernicus (1473--1543)
gave a detailed mathematical model
for planets orbiting the sun
❖
✦
geocentric orbits
heliocentric orbits
Johannes Kepler (1571--1630) was
able to come up with simpler equations
for orbits by allowing planets to move
in elliptic rather than circular paths
wikipedia.com
Newton’s Laws of Motion
✦
✦
According to Isaac Newton (1643 -- 1727):
❖
orbits are not simple circles or ellipses
❖
planets, moons, comets all move
along trajectories that are influenced
by gravity
The force acting between two bodies is
❖
proportional to their masses
❖
inversely proportional to distance
between them
b = make_system(:random, 7, 3)
100.times { update_system(b, 0.1) }
Halley’s Comet
✦
Edmond Halley (1656--1742) was the
first astronomer to recognize that the
bright “star” appearing roughly every 76
years was the same object
❖
recorded by Chinese, Greek, Persian,
others as early as 240 BC
❖
Halley analyzed data from 1531, 1607, 1682
❖
in 1705 he claimed they were all from the
same comet, and predicted it would
reappear in 1758
A Milestone in Computation
✦
In 1758 Alexis Clairaut (1713--1765)
carried out a sequence of calculations
of when Halley’s comet would reappear
❖
started with basic formula for an ellipse
❖
added corrections for gravitational effects
from Jupiter and Saturn
❖
Clairaut and two friends worked for 5 months
doing computation by hand
prediction: orbit of 1682-1758 would be
618 days longer than previous orbit
comet would be closest to Sun in April 1759 (± 30 days)
actual date: Mar 13, 1759
Milestone (cont’d)
✦
Clairaut provided one of the first major
tests of Newton’s laws of motion
❖
planets, moons, etc are not on “tracks”
that follow predefined paths
❖
all movement is determined by the
gravitational attraction between
objects
What Clairaut Did
✦
✦
Clairaut started with a formula for an
elliptical orbit
❖
when the comet was close to Jupiter
or Saturn he computed the effects of
gravity and adjusted the position
❖
he then made future calculations based
on the updated position
The accuracy of this method depends
on how often the adjustment is made
❖
a single adjustment would not help much
❖
too many adjustments would be impossible
to compute
Time Steps
✦
An example of how frequency
of updates affects accuracy
is shown at right
❖
goal: figure out when an
object will hit the ground
❖
use Newton’s law of
gravitational attraction to
compute position
❖
both simulations start with
object at 100m elevation
❖
the simulation that updates
the position more often gives
a more accurate prediction
(4.5 sec to impact vs 4.05 sec)
∆t = 1.35 sec
∆t = 0.10 sec
Equations of Motion
✦
For many situations there are simple
equations that will predict what will happen
✦
A body moving with a constant velocity:
❖
example: Eugene to Portland on the freeway
❖
180 km ÷ 115 km/hr = 1.56 hr
✦
A body falling with constant acceleration:
✦
Using g = 9.8 m/sec2, the time to fall from a
height of 100 meters should be
∆t = 0.10 sec:
t = 4.50
The 3-Body Problem
✦
As soon as Newton’s theory was published mathematicians tried to create
equations that would predict the locations of three bodies
✦
A falling object is part of a two-body system
✦
❖
the Earth orbiting the Sun is another
two-body system
❖
equations can describe the relative locations
of the two bodies
Earth, Moon, and Sun form a three-body system
❖
there are no equations of motion for three-body systems
❖
it is impossible to solve an equation to exactly predict
the future location of any of the three bodies
Aside: Music of the Spheres
✦
Kepler was an astrologer (not unusual at the time) and also a mystic
✦
He was convinced there was a regular
geometric progression in the orbits
of the planets
✦
❖
the farther a planet is from the Sun, the
slower it moves
❖
ratios of maximum and minimum speeds
appeared to be the same as some well
known harmonies
❖
he published his theories in a book called
Harmonices Mundi (Harmony of the Worlds)*
But according to Newtonian physics
planets and other bodies do not move
with any sort of periodic or predictable
motion
* also called “The music of the spheres”
The N-Body Problem
✦
In general, the problem of describing the motion of a series of 3 or more
bodies is known as the N-body problem
the only way to solve the N-body problem is through computation
✦
In mathematical terms, the process is
known as “numeric integration”
❖the
approach taken by Clairaut in 1758
❖today:
computer simulation
✦
Start with known locations of all
the bodies
✦
Iterate:
❖compute
❖calculate
the forces acting on each body
where the body will be
a (small) point in time later
Lab Project
✦
✦
Our goal for this chapter:
❖
a general introduction to computer simulation
❖
discussion of some of the issues of setting up a “model” on a computer and using it to
learn more about a real system
Project:
❖
solar system simulation
❖
use Ruby methods that implement
time-stepped solution of the
motion of planets
❖
experiments with size of time step,
other variables
Walking in Circles
✦
The first phase of the project: learn more about time steps
✦
Experiment:
✦
❖
imagine a robot explorer on a distant planet
❖
we want to send it out on an expedition, have it return to its starting point
We can transmit simple instructions, e.g.
advance(t)
turn(a)
✦
We want to send instructions to have the
robot travel in a circular path
❖
correcting the path takes time, so we want
to use as few instructions as possible
Robot Viewer
✦
Call a method named view_robot to open a graphics window in your IRB
session:
>> include SphereLab
=> Object
>> view_robot
=> true
✦
The robot is the green arrowhead
❖when
we advance the robot it will
move in the direction pointed to
by the arrow
Planting the Flag
✦
✦
Our goal: tell the robot to move in a circular path
❖
plant a “flag” in the middle of the map
❖
move back to the starting point
❖
start moving north, keep adjusting
the heading to move in a circular path
❖
we want the path to be as close as
possible to a perfect circle
The map is a 400 x 400 meter square
❖
the origin (0,0) is at the lower left
❖
the flag should go at (200,200)
Planting the Flag (cont’d)
✦
The robot starts out near the left
edge of the map:
>> robot.location
=> [40, 200]
✦
Some other information about
the initial state:
>> robot.heading
=> 360.0
due north
>> robot.speed
=> 10.0
10 m/hr
the speed it will move when
we tell it to advance
Planting the Flag (cont’d)
remember: d = r × t
✦
distance to center: 200 - 40 = 160 m
To plant the flag:
at 10 m/hr, advance 16 hours
>> robot.turn(90)
=> 90
>> robot.advance(16)
=> 16
>> robot.plant_flag
=> 0
>> robot.turn(180)
=> 180
>> robot.advance(16)
=> 16
>> robot.turn(90)
=> 90
Live Demo
the flag is indicated by
the dot
Tracking the Robot
✦
To keep a visual record of the robot’s path:
>> robot.track(:on)
=> :on
✦
Now when you call advance you will
see a line marking the robot’s
progress:
>> robot.advance(5)
=> 5
The idea of drawing a path using a
robot comes from “turtle graphics”
Moving the Robot
✦
Start making the circle by advancing and periodically adjusting the path
✦
Call turn to change the heading
by a specified number of degrees:
>> robot.turn(45)
=> 45
>> robot.advance(5)
=> 5
>> robot.turn(45)
=> 45
>> robot.advance(15)
=> 15
Live Demo
Challenge: can you move it in a circle?
Making a Circle
✦
A method named orient
tells the robot to turn until the
center point is at a right angle
>> robot.orient
=> 33.08
✦
To make a circle, just tell the robot
to repeatedly move forward a bit
and reorient after each move
>> robot.orient
view_robot( :flag => [200,200], :track => :on )
n.times { robot.advance(dt); robot.orient }
Live Demo
Circle Examples
✦
Two circles made by the advance-and-reorient strategy are shown below
n = 20
dt = 5
n = 50
dt = 2
n.times { robot.advance(dt); robot.orient }
A smaller time step -- i.e. correcting the trajectory
more frequently -- leads to a smoother path
Time Steps and the Solar System
✦
The same effect will be apparent in the solar system simulations
❖
we will be computing a direction for each planet
❖
a planet will move in a straight line in the computed direction
❖
large time steps mean large differences between computed paths and real paths
Orbit of Mercury with
∆t = 3.5 days
Inner planets with
∆t = 35 days
Acceleration
✦
Our robot moves with constant velocity
✦
Constant velocity is not very common in the “real world”
❖
cars need to start from 0 kph before reaching 115 kph
❖
falling objects start slowly and increase speed the longer they fall
✦
A change in velocity is called acceleration
✦
The equation we saw earlier for calculating the distance an object would
fall is an example of constant acceleration
❖
the acceleration is the result of the Earth’s gravity
❖
gravity causes the velocity to change, but acceleration is constant
Newton’s Law
✦
✦
✦
The equation on the previous slide is a special case
❖
it is for a “two-body” system when one body is the Earth
❖
the other is assumed to be small and close to the Earth’s surface
The general case is described by this formula:
❖
G is the universal gravitational constant
❖
m1 and m2 are the masses of the two bodies
Note that the two bodies pull toward each other
❖
the force pulling one body is matched
by an equal force pulling the other body
in the opposite direction
we’ll come back to this
equation later....
Gravity Experiment
✦
To see how gravity affects the motion of a set
of bodies we will do some experiments with
a two-body system
>> b = make_system(:melon)
=> [melon: 3000.0g <...> <...>,
earth: 5.9736e+24g <...> <...>]
An array with 2 objects
Each object is from a class
named Body
Note: 5.9e+24 is Ruby’s notation
for 5.9 × 1024
Urey Hall at UC San Diego
Gravity Experiment Viewer
✦
To set up the drawing canvas to view the
results of the two-body experiments:
>> view_melon(b)
=> true
✦
The melon is initially on the ground
>> b[0].height
=> 0
✦
Call position_melon to raise it off
the ground:
>> position_melon(b, 35)
=> 35
Any height between 0
and 100 meters
Dropping the Melon
✦
A call to update_melon will calculate the
new positions of the bodies in the system:
>> update_melon(b, 0.5)
=> 0.5
✦
The melon falls about 2.5 meters
in the first 1/2 second:
>> b[0].height
=> 32.54
Was 35 after the call to
position_melon
Dropping the Melon (cont’d)
✦
Call update_melon two more times, and
check the height after each call:
>> update_melon(b, 0.5)
>> b[0].height
=> 27.63
>> update_melon(b, 0.5)
>> b[0].height
=> 20.27
✦
Note that the canvas is updated
to show the melon’s position
after each call to
update_melon
The length of each line segment shows
the distance traveled in each time step
height
distance
35.00
0
32.54
2.46
27.63
4.91
20.27
7.36
acceleration...
Dropping the Melon (cont’d)
✦
Raise the melon back up to 35 meters:
>> position_melon(b, 35)
=> 35
✦
A method named drop_melon will
repeatedly call update_melon until
the height reaches 0:
>> drop_melon(b, 0.5)
=> 2.5
The return value is the total time
Live Demo
In this example, the melon fell for 5
time steps of 1/2 second each
Accuracy
✦
Our computer simulation predicts a drop
time of 2.5 seconds:
>> drop_melon(b, 0.5)
=> 2.5
✦
But compare this to the time predicted
by the equation:
Smaller Time Step
✦
Using a shorter time step gives a much
better result:
>> position_melon(b, 35)
=> 35
>> drop_melon(b, 0.01)
=> 2.67
>> Math.sqrt( 35 / 4.9 )
=> 2.672612
Performance
✦
The drop_melon method iterates
until the melon “hits the ground”
❖
we don’t know how many iterations
it will make
❖
we could attach a “counting probe”
to have it count iterations
❖
but we don’t need a probe -- we can
just figure it out from the return
value
>> drop_melon(b, 0.01)
=> 2.67
time = #steps × step size
#steps = time / step size
this simulation ran for
267 time steps
if we could “zoom in”
we would see 267
dashes in this path...
Horizontal Motion
✦
Another way to see that downward
velocity is continually increasing is to
set a horizontal velocity
>> view_melon(b)
erase old dashes
=> true
when they’re spread out
horizontally it’s easier to
see the dashes get longer
on each time step
>> position_melon(b, 100)
=> 100
>> b[0].velocity.x = 5
=> 5
melon moves right at 5 m/sec
>> drop_melon(b, 0.1)
=> 4.5
Live Demo
4.5 sec × 5 m/sec = 22.5 m
Melon Moves Earth
✦
According to Newton’s equations, the melon exerts a gravitational force on
the Earth
✦
Does the Earth move during our simulation?
>> b = make_system(:melon)
=> [melon: 3000.0g ... earth: ... ]
>> b[1].position
=> <0.0, 0.0, 0.0>
(x, y, z) coordinates
of body b[1]
>> position_melon(b, 100)
>> drop_melon(b, 0.1)
>> b[1].position
=> <0.0, 5.10e-20, 0.0>
yes, by 5 × 10-20 m ...
... about 1/30,000th the diameter of a hydrogen atom
Vectors
✦
A convenient way to describe motion in a system is to use vectors
✦
In drawings, vectors are often shown as arrows
✦
Example: velocity of the falling melon
❖
an arrow pointing to the right indicates
horizontal motion (from the “shove” we
give it when we let it go)
❖
a downward pointing arrow indicates
the motion from the force of gravity
❖
the length of an arrow indicates speed
❖
horizontal motion has a constant speed,
but downward motion increases at each
time step (due to acceleration by gravity)
Vector Addition
✦
The overall motion can be found by adding the two vectors
❖
in the diagram, align the tail of one vector
with the head of the other
velocity vectors are used to
update positions at each
time step
Vector Objects
✦
Vectors are implemented in SphereLab by a class named Vector
✦
Each body has a position and a velocity, both represented by Vector
objects:
>> b = make_system(:melon)
=> [melon: 3000.0g <0.0, 6371000.0, 0.0> <0.0, 0.0, 0.0>,
earth: 5.9736e+24g <0.0, 0.0, 0.0> <0.0, 0.0, 0.0>]
>> b[0].position
=> <0.0, 6371000.0, 0.0>
>> b[1].position
=> <0.0, 0.0, 0.0>
>> b[0].position.x
=> 0.0
>> b[0].position.y
=> 6371000.0
The melon starts on the ground,
which is 6,371,000 meters from
the center of the Earth
Vector Objects (cont’d)
✦
We can change a vector by using an assignment statement
>> b[0].velocity
=> <0, 0.0, 0.0>
>> b[0].velocity.x = 3
=> 3
>> 3.times { update_melon(b, 0.1); p b[0].velocity }
<3.0,
-0.98151, 0.0>
<2.9999, -1.96302, 0.0>
<2.9999, -2.94453, 0.0>
velocity.y increases in magnitude
velocity.x stays constant
N-Body Simulation
✦
✦
✦
We’re now ready to see how to simulate the motions of 3 or more bodies
❖
instead of the equation for a falling body
(which uses the constant g) we need the
more general equation, shown at right
❖
G is the universal gravitational constant
The force that pulls two bodies toward
each other depends on
❖
the mass of each body
❖
the distance between the two bodies
Given this information we can compute two
force vectors
❖
the force then determines the acceleration
of each body
Forces Are Additive
✦
The gravitational forces acting on a body can also be described by vectors
✦
An essential concept in n-body simulation:
❖
the net force acting on a body is the sum of all the pairwise forces
longer arrow =
stronger force
Additive Force Demo
✦
To see how other bodies affect one of the bodies in a system:
>> b = make_system(:fdemo)
=> [f1: ..., s0: ..., s4: ...]
:fdemo stands for “force demo”
>> view_system(b, :pendown => :track)
=> true
✦
There are 6 Body objects in this system
❖the
body named f1 will “fall” toward the other 5
❖bodies
✦
s1 to s5 are static -- in this simulation they won’t move
As you step through the simulation, watch how the direction of f1 changes
according to the forces exerted on it by the other bodies
Additive Force Demo
✦
Before starting the simulation, make copies of the first body:
>> f1 = b[0]
=> f1: 1e+13 kg (81.472,145.85,0) (0,0,0)
>> f2 = f1.clone
=> f1: 1e+13 kg (81.472,145.85,0) (0,0,0)
>> f3 = f1.clone
=> f1: 1e+13 kg (81.472,145.85,0) (0,0,0)
✦
A method named update_one will simulate the motion of one body,
leaving the others stationary:
10.times { update_one(b[0], b[1..5], 1.0); sleep(0.1) }
=> 10
Live Demo
Additive Force Demo
✦
Notice how the “falling body” (f1)
moves toward the others:
❖
at each time step the simulator
computes 5 force vectors, one
for each other body
❖
the 5 vectors are added together
to give the cumulative force acting
on f1
❖
f1 is moved ahead in the direction
determined by the cumulative force
❖
the cycle repeats from f1’s new
position
Live Demo
No Equation Predicts this Motion
✦
Because there are more than two
bodies in this system there is no
equation to predict where the
selected body will be at any point
in the future
❖
we need to “numerically integrate”
by computing forces at small time
intervals
❖
update position and recompute
forces at each time step
Chaos
✦
Here’s a short experiment that demonstrates how tiny changes in starting
conditions can have large effects on the final outcome
>> f2.position.x += 1
=> 82.4717
>> f3.position.x += 2
=> 83.4717
>> f2.graphic.fill = "green"
=> "green"
>> f3.graphic.fill = "yellow"
=> "yellow"
>> 10.times {
update_one(f3, b[1..5], 1.0);
update_one(f2, b[1..5], 1.0);
update_one(f1, b[1..5], 1.0);
sleep(0.1)
}
N-Body Algorithm
✦
The complete simulation needs to compute the sum of forces acting on all
n bodies
for i in 0...nb
for j in (i+1)...nb
Body.interaction( b
The force pulling bodies[i]
toward bodies[j] is the same as
the force pulling bodies[j] toward
bodies[i]
N-Body Algorithm
✦
After calculating the forces, for each body:
❖
compute acceleration due to sum of forces acting on it
❖
update velocity and position vectors
bodies.each do |b|
b.move(time)
# apply the accumulated
Type Source.listing(:step_system) to
see the complete Ruby code
Solar System
✦
To check the accuracy of the watermelon experiment we had an equation
❖
run the simulation, e.g.
>> drop_melon(b, 0.5)
=> 2.5
❖
✦
compare result with
For the N-body methods, we can test the program by having it simulate the
motion of the planets
❖
our evaluation will be an informal visual inspection -- just make sure they move in elliptical
orbits
❖
a better test would be to compare simulated motion against actual measurements
Solar System
✦
Call make_system to set up the experiment:
>> b = make_system(:solarsystem)
=> [sun: 1.9891e+30g ... , pluto: 1.314e+22g ...]
✦
The call to make_system will create an array of 10 bodies
❖one
for the Sun and each of the 9 planets (including Pluto)
❖initial
positions and velocities are from an ephemeris service at JPL
Solar System Dynamics: http://ssd.jpl.nasa.gov/
❖data
✦
from Jan 1, 1970
Call view_system to draw the bodies:
>> view_system(b, :dash => 1)
=> true
(result on the next slide)
Solar System
The method that displays
the bodies automatically
sets the scale so all of
them are visible
The circles are not to scale...
Solar System
It’s easier to see orbits if
you view just the Sun and
inner planets (bodies 0
through 4)
view_system(b[0..4], :dash => 1)
Solar System
Here is a screen shot from
a simulation using a time
step size of 1 Earth day
(86400 seconds)
Earth makes one full orbit
Mercury and Venus have
shorter years -- at this
scale the orbits overlap
Live Demo
365.times { update_system(b, 86459) }
Solar System Experiment
✦
✦
✦
One of the claims made early in the chapter:
❖
planets, asteroids, and comets have elliptical orbits because the Sun is so much larger
than any other body
❖
Kepler was mostly right -- orbits are (generally) ellipses
❖
but we need Newton’s equations to get the finer details or to plot trajectories for
spacecraft
An experiment to test this claim:
❖
increase the mass of one of the bodies, e.g. Venus or Mars
❖
see what effect this has on the motions of all the other bodies
Before we try it out -- what is your prediction?
Live Demo
Solar System
Chaos
Can you explain what’s
going on here?
view_system(b[0..4], :dash => 1)
Scalability
✦
✦
Using notation we saw earlier in the term,
we would say this algorithm for solving
the n-body problem is
❖
it has the same basic structure as insertion sort
❖
one iteration nested inside another
More advanced algorithms are
❖
✦
for i in 0...nb
for j in (i+1)...nb
....
end end
use a “binary tree” instead of an array to organize bodies
Other methods use ideas from fluid dynamics
❖
next slide: QuickTime movie of collision of two galaxies (our own Milky Way and
Andromeda)
❖
based on a simulation with 100,000,000 bodies, ran for 4 days on a supercomputer with
1100 CPUs
QuickTime™ and a
YUV420 codec decompressor
are needed to see this picture.
http://www.galaxydynamics.org/
Modeling and Simulation
✦
Computer modeling is a very powerful method for helping to understand
complex situations
✦
Widely used in
❖
science (e.g. “computational science”)
❖
engineering
❖
medicine
❖
pharmaceutics (“rational drug design”)
❖
business (“computational finance”)
Computer Generated Imagery (CGI)
✦
Animation and special effects
generated by computers are
other examples of modeling
and simulation
❖
the images of dinosaurs in the
movie Jurassic Park were all
created by simulation
❖
for each scene, the dinosaurs
were represented by objects
in a computer model
❖
“ray-tracing” algorithms computed
how light would reflect from these
objects and real objects in the scene
Note the reflection of the dinosaur on
the floor and on the side of the table...
Issues in Modeling and Simulation
✦
✦
How reliable are computer models?
❖
are the values produced by the computer simulation accurate?
❖
what are potential sources of errors (difference between reality and simulation)?
Some places where errors may be introduced:
❖
software -- simple programming bugs
❖
discretization -- size of time step, size of cell in a wire frame, etc
❖
round-off errors -- floating point numbers have a finite resolution
❖
modeling errors -- leaving out an essential component of the real system
❖
‣
‣
example: assuming bodies have a constant mass
comets lose a little of their mass on each approach to the sun
using the wrong scientific or mathematical model
‣
‣
‣
example: Newton’s model does not explain changes in Mercury’s orbit
full explanation requires ideas from general relativity
see “precession (astronomy)” at Wikipedia...
Summary
✦
✦
✦
The ideal situation for understanding a system is through a set of
mathematical equations
❖
use analytical techniques to solve the equations
❖
examples:
For many systems it is not possible to solve the equations analytically
❖
“numeric integration” provides an estimate
❖
example: motions of 3 or more bodies due to gravity
Computer simulation uses calculations based on the mathematical models
❖
for n-body systems, use Newton’s law to compute pairwise gravitational effects
❖
for each body, add the effects induced by all the other bodies