Scientific Computing in Physics

Why Physics Uses Computing •  Some physics problems have simple equations that can be solved “analytically” (i.e. with pen and paper). These usually involve simple forces and few objects being acted upon: •  e.g.: ball under the influence of gravity: F y Solution: g

m

velocity: position:

= −mg v(t) = v(0) − gt

1 2 y(t) = y(0) + v(0)t − gt 2

This is a very IDEALIZED problem: •  assumes no drag (friction) due to the air •  only 1 object, force acts on entire object equally (body force)

Why Physics Uses Computing •  Forces can be much more complicated: •  e.g. fluid in Earth’s core:

   1   2 F = −∇P +ρ g + µ∇ v + ∇(∇ ⋅ v) − 2 ρΩ × v 3

•  System may include many components, each affected by different forces: •  e.g. creating the moon from collision of Earth & Mars-sized body:

•  These problems can’t be solved with pen & paper. •  Instead we need the rapid calculation power and large memory of the biggest computers available

How Does a Computer Solve Physics Problems? •  Computers don’t know about physics (or any other science) •  They are good at basic math operations (addition, multiplication, exponents…). They are also very fast at doing these operations. •  Recommended plan: 1)  Determine the equations governing the physics we are studying 2)  rewrite them in terms of simple mathematical operations (there may be MANY of these operations but that’s okay because computers are fast, e.g.: •  Fastest computer in the world (Top 500 list) can do ~2x1016 operations per second! •  Scinet here at UofT (#66 on Top 500 list) can do ~3x1014 operations per second! 3)  tell the computer to do the mathematical operations

Lets try it on a simple example…

How Does a Computer Solve Physics Problems? •  Example: Ball under the influence of gravity: F

= −mg

Keep in mind we know the solution already: g

m

velocity: position:

v(t) = v(0) − gt

1 2 y(t) = y(0) + v(0)t − gt 2

•  Lets pretend we don’t know the solution (since we won’t for other problems). How do we get the computer to calculate the solution? 1)  Determine the equations governing the physics we are studying: (these equations true as long as Δt very very small)

Newton’s Law:

F = −mg F a(t) = m

Δv a(t) = Δt Δy v(t) = Δt

Δt : small interval of time Δv : change in velocity over Δt Δy : change in position over Δt

How Does a Computer Solve Physics Problems?

F = −mg F a(t) = m

v

Δv a(t) = Δt Δy v(t) = Δt

vf

Δv vi

ti

Δt

tf

t

How Does a Computer Solve Physics Problems? 2)  rewrite them in terms of simple mathematical operations

F = −mg F a(t) = m We get:

Δv a(t) = Δt Δy v(t) = Δt

v f = vi − gΔt y f = yi + v f Δt

Rewrite our equations for a(t) and v(t) using:

Δv = v f − vi Δy = y f − yi Notice these aren’t the analytic solutions:

v(t) = v(0) − gt

1 2 y(t) = y(0) + v(0)t − gt 2

These are now equations using basic math operations on how to update the velocity and position after each time step.

How Does a Computer Solve Physics Problems? 3)  tell the computer to do the mathematical operations

v f = vi − gΔt y f = yi + v f Δt 1.  2.  3.  4.  5. 

Start with t=0, vi, yi defined at t=0. Define Δt Take 1 step in time (1Δt ) Calculate new vf and yf from formulas above Repeat steps 3 – 5 for however many time steps we need

Lets use the python programming language to do this…

Free Fall Code •  click on the “vidle for vpython icon •  under the “File” menu open “freefall.py” in the folder “scicomp” on the desktop •  Save as: freefall_v1.py (in case you need to start over) •  Run your program using “F5” or from the “Run” menu •  You should get the following:

Free Fall Code •  Lets look at the program: (Any line starting with # is a “comment”, the code does not read it) ball = sphere(pos=(0,10,0), radius=0.5, color=color.red)! ground = box(pos=(0,-0.1,0),size=(24,0.2,5),color=color.green)! •  These lines create the ball and ground ball.vel=vector(0,0,0)! ball.acc=vector(0,-9.8,0)! •  These lines define the ball’s initial velocity and acceleration

Since that is the only part of the code active right now, no time steps are taken so the ball never updates its velocity or position and stays at the same location (that’s why nothing happens)

Free Fall Code •  Okay, so lets make something happen: 1.  Start with t=0 2.  Define Δt •  In the code, remove the # sign in front of the following lines: time=0.0! dt = 0.005! ! 3.  Take 1 step in time (1 Δt ) 4.  Calculate new vf and yf •  In the code, remove the # sign in front of the following lines (make sure to keep the indentation): time=time+dt ! ball.vel.y = ball.vel.y + ball.acc.y*dt! ball.pos.y = ball.pos.y + ball.vel.y*dt! just these formulas:

v f = vi − gΔt y f = yi + v f Δt

Free Fall Code 1.  Repeat steps 3 – 4 for however many time steps we need •  The best way to do this with a computer is using a “loop”. A “while” loop repeats all the commands indented in the lines after the while statement as long as some condition is met. •  Remove the # sign in front of the line: while ball.pos.y >= 0.0:! •  this starts our loop. The commands indented afterwards: time=time+dt ! ball.vel.y = ball.vel.y + ball.acc.y*dt! ball.pos.y = ball.pos.y + ball.vel.y*dt! •  will be calculated each time through the loop •  the loop will keep repeating until the position of the ball < 0 (the ball hits the ground). •  Run your program.

Free Fall Code •  Did it run so fast you couldn’t see it? The code has its own “units” of time which is pretty fast so we can slow it down by using a “rate” command. Remove the # in front of the line: rate(1.0/dt)! •  Now run your program. You should be able to watch the motion now. •  Try running your program with different initial y velocities and different initial y positions by changing the following lines: ball = sphere(pos=(0,10,0), radius=0.5, color=color.red)! change these numbers ball.vel=vector(0,0,0)!

Projectile Motion Code •  Open the program “projectile.py”. Save it as “projectile_v1.py”. •  Notice that currently it is almost identical to the free fall code. If you run it you will see it is the same. •  Projectile motion also allows initial velocities in the “x” (horizontal) direction so instead of dropping the ball we will throw the ball. •  Your mission: Add some new lines to this code to give the ball an initial velocity in the x direction and to update the x components of the position and velocity in time. •  For the initial x velocity you will want to alter the following line (the value is up to you, you can change it later): ball.vel=vector(0,20,0)!

x-component

y-component

Projectile Motion Code •  Is there any acceleration in the x direction?

Projectile Motion Code •  Is there any acceleration in the x direction? No, gravity only works in the y direction. So in the x direction:

v f = vi x f = xi + v f Δt •  Do we ever need to update the velocity in the x direction?

Projectile Motion Code •  Is there any acceleration in the x direction? No, gravity only works in the y direction. So in the x direction:

v f = vi x f = xi + v f Δt •  Do we ever need to update the velocity in the x direction? No, its always the same. It’s a constant = initial x velocity. •  Do we need to update the position in the x direction?

Projectile Motion Code •  Is there any acceleration in the x direction? No, gravity only works in the y direction. So in the x direction:

v f = vi x f = xi + v f Δt •  Do we ever need to update the velocity in the x direction? No, its always the same. It’s a constant = initial x velocity. •  Do we need to update the position in the x direction? Yes! According to the above formula.

Projectile Motion Code •  Here is the line in the code for the y position. ball.pos.y = ball.pos.y + ball.vel.y*dt! •  Create a very similar line for the “x” position (replace all the y’s with x’s) and add it to the code right after the above line. Don’t forget to keep the line indented! •  Run your code. Did you get projectile motion? •  Try varying your initial x and y velocities to see what you get. •  You have now created a projectile motion code! Recap: •  We defined the laws of physics for our object •  We rewrote the laws into simple mathematical operations •  We made the computer do the operations over and over again to update the motion of the object.

Contact Forces •  gravity is an easy force to work with because it acts at all times and on all of the body •  What about modeling another type of force? For example: •  Contact force due to a collision: Very hard! But luckily we have some physics “conservation” laws to help us.

A

vA=constant

What happens when they collide?

B

vB=0

Physics of Collisions It depends! •  What are the balls made of? e.g. Tennis balls collide differently than balls made of peanut butter (try this at home) •  We will consider a special kind of collision: an “elastic collision”. In this case, the energy is conserved in the collision

before collision

after collision

vA,i=constant

A

vB,f

A

vB,i=0

B

B

vB,f

Physics of Collisions •  Use momentum conservation and energy conservation laws to determine final velocities: momentum = mass * velocity

p = m*v Momentum Conservation Law: Consider all objects in system. If no outside forces acting, then total momentum of system stays constant even if the objects collide.

ptot = mA * vA + mB * vB

is constant in time

à Initial momentum (before collision) = Final momentum (after collision)

ptot ,i = ptot , f ⇒ mA vA,i + mB vB,i = mA vA, f + mB vB, f

Physics of Collisions Momentum conservation: For an “elastic” collision:

mA vA,i + mB vB,i = mA vA, f + mB vB, f

1 1 1 1 2 2 2 mA v A,i + mB vB,i = mA v A, f + mB v 2B, f 2 2 2 2

We will only consider collisions where vB,i = 0 (it makes the math easier). Then we can solve the above two equations for the final velocities after the collision:

vB, f

2mA = vA,i (mA + mB )

vA, f

mB = vA,i − vB, f mA

Lets work on a code to model these elastic collisions with the above equation. Open “collision.py”. Save as “collision_v1.py”

Collisions Code •  Code creates 2 balls (ballA & ballB): ballA=sphere(pos=(-2,0,0),radius=0.5,color=color.red)! ballB=sphere(pos=(2,0,0),radius=0.5,color=color.blue)! •  sets the mass of each ball: mA=1.0! mB=4.0! •  defines the balls’ initial velocities: ballA.vel=vector(2,0,0)! ballB.vel=vector(0,0,0)! •  then updates the positions of the balls using a while loop: while 1==1:! rate(1/dt)! ballA.pos=ballA.pos+ballA.vel*dt ballB.pos=ballB.pos+ballB.vel*dt! •  Run the program, what happens?

!!

!

!

Collisions Code •  The balls don’t collide! Well that’s because we haven’t told the computer about the collision. •  To model the collision we need to tell the computer : 1.  how to determine if a collision has occurred 2.  how to calculate the velocities after the collision Determining when a collision occurs in terms of things the computer knows: A B R ballA.pos

ballB.pos

Think about a condition on the positions of the balls that would mean a collision occurs. Write your equation on the whiteboard once your group has agreed on it.

Collisions Code

B

A R

ballA.pos

ballB.pos

ballB.pos.x – ballA.pos.x = ballA.radius+ballB.radius ! •  We only want to change the velocities if a collision occurs. So we will use an “if” statement to determine this •  In the “while loop” remove the # in front of the line: if ?????:!

Collisions Code if ?????:! •  Replace the ????? with: ballB.pos.x – ballA.pos.x == ballA.radius+ballB.radius ! ***Make sure to leave the “ : ” after the question marks*** ***Notice the double “ == “ *** •  Now lets put in the equations for the velocities after the collision:

vB, f

2mA = vA,i (mA + mB )

vA, f

mB = vA,i − vB, f mA

ballB.vel=2.0*mA/(mA+mB)*ballA.vel ! ballA.vel=ballA.vel-mB/mA*ballB.vel ! •  Run your program. What happens?

Collisions Code •  The balls don’t collide!!!!!!

Why not!!!!!!

•  The answer lies in how computers store data and do mathematical operations. Basically, small errors creep (around the 16th decimal place!) so that makes the == condition impossible to match. ballB.pos.x – ballA.pos.x == ballA.radius+ballB.radius ! •  Luckily, there is a way around this. We can say a collision occurs when the condition below is true for the FIRST time: ballB.pos.x – ballA.pos.x