(Pramod Gupta)
Hi, it's great to be here in Portland.
And thanks for stopping by
on this sunny afternoon
for a talk on computational physics.
So the doors are still open
in case somebody came here by mistake.
There still is some time.
Now, how many people
know who Feynman is?
OK, almost everybody.
There are a few who don't.
So for them I'll just, you know,
introduce them who this guy is.
Everybody knows who Newton is,
but I think it's interesting --
[bongo drums being played
in YouTube video]
-- to see who Feynman is.
So he's a guy who is a legend
among the physics community.
But it's not only just his physics
were extraordinary,
but he was, you know, quite a character.
So if you don't know who he is,
you should probably check out Feynman.
You can search for his videos
on gravitational physics and so on.
But he was quite a bongo player also.
And his great lectures on physics,
which he taught in Cal Tech in 1961,
is a classic of physics.
And a few years ago,
Cal Tech put it all online.
So anybody now can go
and check out those lectures
and see, you know, his teaching.
So let's now go
to the main part of the talk...
...which is about...
...using Python
for calculating planetary orbits.
And...so I'm Pramod Gupta.
I'm a research scientist
at the Department of Astronomy
of the University of Washington.
And the thing we'll focus on today
is Kepler's first law:
The planets move around the sun
in an elliptical orbit
with the sun at one of the two foci
of the ellipse.
So Kepler found this
after many years of extraordinarily hard work
and brilliant thinking
from Tycho Brahe's data.
And you need to keep in mind
that the data was taken
with respect to the earth,
and Kepler had to think about how it is
with respect to the sun,
which was an extraordinary problem.
And just as a reminder,
this is what an ellipse looks like.
It's sort of like a squished circle,
but it has a very precise definition.
The distance -- the sum
of the distance from the two foci
is a constant.
And so Newton applied his laws of motion
and the law of gravity
to the orbit of the planets
and derived Kepler's laws.
And this again
was another extraordinary piece of work.
But the derivation
requires advanced mathematics,
so even hundreds of years later,
it's not an easy problem.
And it's generally taught
to physics students towards,
you know, maybe their third
or fourth year of their study.
And the question is that:
is there some other way?
This is Feynman's blackboard
from when he was teaching his lectures.
And it's interesting
what he points out in the last line.
Now, you know, Feynman was not --
he made a lot of great discoveries
but this particular point --
thing which he's saying,
that is not his discovery.
It was well known before.
But usually you learned about it
if you went on to do advanced physics.
You did not learn about this
in your freshman physics class
where you did mostly much more
simpler problems.
So this is what he's saying.
So for those of you
who have forgotten
cursive handwriting,
he's saying,
"If the force law is known
"as a function of position
and velocities..."
"...then the equations --"
so he's referring to Newton's equations --
"can be solved in close approximation
step by step in time by arithmetic."
OK, so he's saying
we can do this all by arithmetic.
And, you know,
all of us here know arithmetic
from our grade school
and elementary school.
And so,
how should we go about
this calculating of the orbit?
And, you know, Feynman tells us,
"This work can be done rather easily
"by using a table of squares,
cubes, reciprocals.
"And then we need only
multiply x by 1 by r cubed,
"which we do on a slide rule."
So, how many people here
brought slide rules with them?
[laughter]
Oh, there are a few still.
OK, so for those of you
who don't know
what a slide rule looks like,
so this is the thing.
So the secret is: get a slide rule.
But we have computers.
So the computer can do arithmetic.
And if you have
a few hundred million dollars,
this is the computer
you want to buy.
But the amazing thing is that computers
have become so much faster
over the past 50 years since Feynman
gave his lectures
that even the regular
off-the-shelf laptop
can do those computations
pretty quickly.
But we must tell it what to compute.
So...
So the fun and games part
of the talk is over,
and we come to Newton's law of gravity.
And essentially what Newton
is telling us is that
the force between two objects,
in this case a planet and the sun,
is proportional to the product
of their masses.
The big M is the sun.
The little m is the planet.
And r is the distance of the planet
from the sun.
So that tells us the magnitude,
the numerical value of the force.
But it's a vector, which means
that we need to think
about the direction.
So the direction of the force
is from the planet to the Sun,
which this arrow --
which this arrow over here shows.
And we need to take two components.
So from your high school
geometry class, you can see
that this big triangle over here,
this triangle,
and this little triangle over here,
the force triangle
and the position triangle
are similar triangles.
And you can use that
to actually show
that the x component
of the force
is the magnitude of the force,
and then you multiply it by x divided by r,
and similarly
for the y component of the force.
Now Newton not only had
his law of gravity,
he also has a second law,
F = ma, which is --
which, along with e = mc squared,
is probably the most famous
equation in physics.
You know, you see it on t-shirts
and places like that.
And you need to equate these two forces.
And what happens is
that the little m drops out,
which actually is sort of
an amazing coincidence.
And it was -- it took
over 200 years later
for Einstein to actually explain
why that happens.
But in any case,
you get an expression
for the acceleration,
the x component,
and the y component
of the acceleration.
So we have done some algebra here
which is fairly straightforward.
The question is,
how do we use this information
to compute the orbit?
So we'll do a few more algebraic
manipulations just for convenience.
And I'm closely following chapter 9
of Feynman's lectures.
So I highly recommend, you know,
you go and check it out,
because of his incredible,
very readable way of writing.
He doesn't -- he just tells you
just the right amount of details
which are required
to solve this problem.
So we just use units
such that big G and big M are 1.
This is just for convenience in calculations,
but the computer doesn't care.
You could use the actual big G
and big M values.
And we have
the acceleration over here,
the x component,
it becomes slightly simpler looking.
and the y component.
And one more thing we need to
remember is Pythagoras' theorem.
So the distance from the sun
to the planet
is given by the square root
of x squared plus y squared.
And these are just
a few definitions, really,
that the position
is given by x and y.
The velocity of the planet is given
by Vx and Vy.
And the acceleration is given
by Ax and Ay.
So we just saw that we already
know the acceleration
from Newton's two laws.
In this slide, all you need
to keep into account is this.
So if you remember your calculus,
the other stuff will make sense to you.
But all you need
to remember is the --
one, two, three, four --
the fourth line.
So if you know the position
of the particle
(or the planet, in this case)
at time t
and you know the velocity,
you can calculate the position
at a time t + ϵ (epsilon),
where ϵ is the traditional
mathematical notation,
a Greek symbol for a very small value.
And similarly, you can find the velocity
at a time t + ϵ,
which is again
the fourth line over there.
It's not showing up
in highlighting there.
And what this tells you is
that if you know the position
and the velocity
at some given time,
you can use
this arithmetic procedure
to keep calculating it
for the next slightly later time,
and then the slightly
later time and so on.
And that's what we'll do.
Now, Feynman uses leapfrog integration,
which is more accurate.
And it also has some nice properties.
It conserves energy
so that if you go around the whole orbit,
you don't lose energy.
You don't spiral into the sun.
And it also has time reversal symmetry,
which means that if you start
from a given position
and do this calculation,
and then you change the sign of ϵ,
you make it ‒ϵ,
and you start going back, you'll come back
to the same position that you started.
So a lot of other techniques
don't have that property.
Now, Feynman doesn't mention it
because he's --
you know, he's a great teacher,
but I just [indistinct]
putting these two lines here.
And what this means is
that if you do leapfrog integration,
you do the same thing we did earlier.
But we calculate the velocity
at a time offset by ϵ/2.
What that means is that you do
the position calculations
at integer values of ϵ
and you do the velocity calculations
at half integer values.
So ϵ/2, 3ϵ/2 and so on.
And so you can see why
they call it leapfrog.
Essentially you go --
the position values leapfrog
over the velocity values,
and the velocity values leapfrog
over the position values, OK.
And this gives you
much more accuracy.
For those of you who remember
your calculus classes,
it is because you're considering
the velocity
at the middle of the time interval
and not at either end.
And so we know the initial position
and velocity,
and using that, we can calculate
the acceleration over here.
And we can also calculate
the velocity of ϵ/2.
OK. So we need this
offset value of the velocity,
as I just mentioned,
because we are going to do
this leapfrog integration
just like Feynman did it.
And so this is the first iteration.
It might it might look a bit complicated
but it's fairly simple.
You can think of it kind of like
in a recursive sort of way.
So the position at a time ϵ
is equal to the position at a time 0
plus ϵ times the velocity at time ϵ/2.
OK. So it's -- you know, if you know --
it's like a lot of things in programming.
If you know x for i,
then you can know x for i + 1
using some assignment statement.
And this is just given here
for completeness,
which is the general case,
but I'm not going to talk about it
since it's much easier actually
to just look at this formula
and do it manually
using a calculator
or a computer,
as we can easily do now.
And I'll choose the initial conditions
which Feynman used.
Basically what he used
was some kind of conditions
which are very convenient
for the calculation,
which is essentially saying
that you start the planet
at x value of 0.5
and a 0.0 value for y.
So it's right on the x-axis.
And you give it some velocity
in the y direction.
So it's at the x-axis
and, you know, you give it some value
in the y direction for its velocity.
And one more thing we need to do
is choose the step size epsilon.
So, Feynman -- so, note that epsilon
is nowhere in the physics, OK?
Newton didn't know
anything about epsilon.
And Feynman used epsilon = 0.1
to minimize computation.
So -- but the smaller
the epsilon value that you use,
the better accuracy that you get.
And since computers are so
incredibly fast nowadays,
we'll use epsilon = 0.001
for the same initial conditions
that he used.
For other initial conditions,
one may need an even smaller epsilon.
And this is one of the things
about computational physics,
that you have all these epsilons
and mesh size and things like that.
And you need to make sure
that they are small enough
for the problem that you're looking at.
And this is Feynman's solution.
And it's very interesting, you know.
So they did it manually back in the day
with slide rules and things like that.
And, you know, to save --
well, slide rule time, not computer time,
to save the slide rule time,
he just did it for a half orbit.
So he didn't do it all the way. So it's
not necessarily obvious from this
that this is a closed orbit.
But you can see a lot of
interesting features here,
that when you're close to the sun,
it travels longer distances here.
And when it's far away
from the sun,
the planet travels
much shorter distances here.
And that's related
to Kepler's second law,
which I'll not talk about.
But you can also get
Kepler's third law
doing these kind of calculations.
And now what I would like to do
is talk a bit about the code.
Since all of you use Python
on a daily basis,
you might be wondering:
"All this is fine, but how complicated is it
"to put it actually in real Python code,
in a functioning code,
"not just hand-waving arguments?"
And so all we need to do
is set up epsilon
and the number of steps
you want to take in your loop.
So that's just one simple thing
we need to do.
We use the initial condition.
So we set up -- you know,
all of you do a lot of programming,
so you know
that initialization is very critical.
And you need to initialize
things properly.
And we'll use the same conditions
that Feynman used.
But this is just so that we can
match up our calculation
with his calculations.
But you can always, you know,
play around.
That's the fun part of all this,
that you can write your own code
and then play around
with these initial conditions.
We allocate some memory
for the lists we'll use.
We initialize the starting values.
And again, this is just
putting into Python code
the stuff which I wrote earlier
in mathematical form.
And you can see that it's
very straightforward.
You initialize your position,
you initialize your acceleration,
and you initialize your velocity.
Note that the velocity
is offset by ϵ/2.
So the 0 index for velocity
means that you're already
ϵ/2 progressed in time.
And this is the main loop.
So this is really
the only main critical part,
where we are using --
we are calculating x of i,
when we already know x of i ‒ 1
and Vx of i ‒1.
So essentially, as Feynman is saying,
you know, all we need is arithmetic.
All that the computer is doing here
is some straightforward arithmetic.
And it's doing it many times.
We are doing it 10,000 times.
Feynman in his calculation
did it, you know, maybe 20 or so,
so a much smaller number of times
as you saw in that plot.
And at the end, we can plot the orbit
with these matplotlib commands.
And one of the motivations
for this talk is the fact that --
well, one of them was
the lectures are now available online
for anybody to look at,
and the other motivation is that it's --
over the last few years,
you have these nice packages,
nicely built distributions for Python,
like Anaconda and Canopy and so on.
So, you know, you can install them
anywhere and you automatically get
all these things,
like matplotlib and so on.
And you don't need to be a big expert
in installing stuff.
With that, I would like to first,
you know, run some of the code.
So this is going to be a bit --
it's always interesting
doing a live running of code.
I'll just show you the code
for a second so that...
You can see it's essentially
the same code
that I showed you in the slides.
But I just want to show it also
in the actual core,
so you know this is
a core we'll do a live run of.
So there's no extra fancy stuff.
And...
Let's run it for Newton's...
...er, for Feynman's initial conditions.
And you saw it took
less than a second,
something which would have
taken a long time.
And you can see
it's a nice ellipse shape.
And that tells you actually
that Newton's laws, his law of gravity
and his law of motion
actually contain these things.
As Feynman used to like to say
that Schrödinger's equation --
that's the equation
for quantum mechanics --
contains frogs, but we can't solve it
in enough complexity
to actually show 
that it, you know, contains frogs.
But here we can show
that Newton's equations
contain Kepler's laws.
And we have a nice ellipse shape here.
And one thing I like
about matplotlib
is that you can do
these zooming things.
Say, if you click here,
you can zoom into this rectangle.
And you can see
that this (0, 0) position
is where the sun is.
That's one of the foci.
And there's another foci
somewhere here.
And you get a nice elliptical shape.
But this is not all you can do
with all this.
You can actually say
what happens.
So this is much easier this way
with the computational approach
than with the...
...necessarily than with the more
advanced mathematical approach.
What you could do is, you could say,
"What happens if I increase the..."
"...velocity?"
And you can see
the planet escapes orbit.
So you can actually see
that if your planet has too much energy
or if whatever is orbiting
the sun has too much energy,
it's going to escape.
And what happens if you...
...decrease the velocity,
which means
you're decreasing the energy also?
OK, that looks weird.
Now let's see what that means.
If you zoom in...
...you find that it's
not a closed orbit.
So what's happening is that the orbit
seems to be shifting.
And it never comes back
to the same point.
But we have not discovered
any kind of new effect here.
The problem is...
...that our epsilon
is too big for this particular case.
So we need to decrease epsilon.
And now we get
this nice ellipse shape.
And note that it's...
...it's very elliptical,
so to speak.
This is the position
of the sun over here.
So when it's close to the sun,
this planet --
so it's unlikely a planet
would ever orbit like this,
but things like comets and so on
have very elliptical orbits
where they go off way out very far away
and then they come back.
They're very close to the --
relatively close to the sun.
So this is a kind of orbit
which may not be
so common for a planet,
but it could happen
to other kind of objects.
So let's go back to...
...the presentation.
So we saw the code.
And I'd like to conclude
with Feynman's own words here,
where he says that, "Now,
armed with the tremendous power
"of Newton's laws, we can not only
calculate such simple motions."
So he was talking about a particle
on a string and so on.
But he's saying, "We can not only
calculate such simple motions,
"but also, given only a machine
to handle the arithmetic --"
and nowadays,
all of us have such a machine
to handle the arithmetic --
"even the tremendously complex
motions of the planets,
"to as high a degree
of precision as we wish!"
So essentially, in his Chapter 9
he shows how you could do
the same thing
for the whole solar system,
which in his days
was like nine planets.
Now we only regard it
as eight planets.
But basically, you can put in
all the planets, calculate all the forces,
do this computation
step by step,
and make a model
of the solar system,
a realistic model,
on your computer.
And so with that,
I'd like to conclude.
And I hope some of you
at least go and look up
Feynman's Chapter 9.
I'd like to conclude this talk
and open it up for questions.
[applause]
(session chair)
So we are open for questions.
(audience member)
Hi. That was really interesting.
Thank you for the talk.
One question I had is,
you showed the ellipse
that sort of degrades
because your epsilon
wasn't small enough
to simulate it correctly.
Do you have any thoughts or are there
any, like, standard ways of telling --
if you get an unusual result
from a numerical simulation
that's not what you expected,
are there ways of knowing
whether that's a fault
with your simulation code
if it's a fault with your parameters?
Like, how do you tell the difference
between an unexpected effect
that actually really is there
and an unexpected effect
that's caused by how you happen
to be doing a simulation?
(Pramod Gupta)
Right, so that's an important question.
And generally in computational
physics, often you don't --
So here we knew the answer, but
generally you may not know the answer.
And the way to tackle that essentially
is to vary these parameters
like epsilon, which are not
really there in the physics.
So if you see something funny,
you try something finer, say.
So, you know,
you decrease your epsilon.
And unless you notice a pattern
where no matter how the small epsilon is,
you keep seeing the same thing,
then you think
that maybe this is a real effect.
But generally that's the way
you go about solving a problem
which you don't know the answer to.
(audience member)
Thank you.
(audience member)
Hi. So this might be a long answer,
but you mentioned
that the leapfrog method
preserves energy, and suppose...
So, does -- didn't -- if you don't
use the leapfrog method,
do not preserve energy
even if you use, like, a lower accuracy?
And if so, why does
the leapfrog method
work so well
at preserving energy?
(Pramod Gupta)
So the thing about the leapfrog method
is that it has this, you know,
time reversal symmetry.
So, optimal machine precision, OK.
So -- and that is what causes it
to have this extra property
of preserving the energy.
So if you calculated the energy,
which is, you know,
1/2 mv squared plus the
potential energy due to gravity,
it would not trend to a lower
or possibly even a higher way.
It would stay steady.
Maybe a little bit of oscillation,
but it will not, like, you know...
...create a situation where
your planet sort of goes into the sun.
And so essentially...
...compared to a lot
of other techniques,
this particular thing has --
the leapfrog method has this property
because of the
time reversal symmetry.
(audience member)
All right. Thanks.
(session chair)
Any other questions?
So, let's thank Pramod. Thank you.
[applause]
