# Intro to Jupyter and Projectile Motion

I wrote this piece as part of a presentation I gave to the robotics team at Benson High School in Portland, OR. My thanks to my mentor John Delacy for working with me in high school and fostering my love for science and astronomy, and for inviting me to help work with the next generation of scientists and engineers. I loved getting to see their robots and hear about the competitions.

# Intro to Jupyter Notebooks

Currently we are using a jupyter notebook. This format can support either Julia, Python, or R. The format is quite similiar to that used in the propriety software Mathematica.

We have two types of cells:

- Markdown cells
- Code Cells

### Markdown Cells

This is a **Markdown Cell**. Here we write text or add pictures to describe what we are doing. In raw format, everything is easily readable, but the computer can also render it to look even nicer with headings, lists,**bold**,*italic*, tables, links, and other formatting.

I can even write equations.

Equations use LaTeX syntax. LaTeX is a document preparation language used by those in STEM. It’s how we make all those awesome looking papers. It also automatically handles our citations for us :) Check out Overleaf to create documents online.

I can also write code inline `x+=1`

or in block style

```
for ii in 1:n
x+=1
end
```

### Coding Cells

You can evaluate coding cells by pushing **Shift-Enter**.

Next to Coding Cells, `In[ ]`

specifies the code that will be executed, and `Out[]`

specifies the output of the given cell. The number specifies the order in which the code got evaluated.

**Important**: Global variables get edited each time you evaluate a new cell. So order of evaluation matters.

### Julia’s Packages

Julia uses external packages, much like Python, to supplement its core functionality. **PyPlot** imports plotting functionality from Python’s PyPlot with the package PyCall. Julia maintains a complete list of supported external packages. While many are quite specialized or out of date, some are incredibly useful, like different plotting tools, curve fitting, integrators, differential equations, accelerators, working with different file types, and more.

If you need to run a jupyter notebook on your machine instead of on JuliaBox, you would be using the IJulia package.

```
# Update packages
Pkg.update()
# Add a new package to your repository
# Unnecessary if already added
#Pkg.add("PyPlot")
# Load PyPlot into
using PyPlot
```

# Projectile Motion

BUT

NO! How do we solve this?

First let’s write down the infinitesimal, exact version of the equations. Sorry if you haven’t had Calculus yet. Just block out this part and keep going.

In order to put this into an equation, we take the derivative and break it into a courser grained version

Now because $\Delta x$ and $\Delta t$ have actual sizes instead of being infinitesimally small, the computer can deal with them.

Now we take lots of baby steps of $\Delta t$ over our time interval to change the position with the changing velocity.

We can also think of this as finding a small enough interval such that we can treat the y velocity as if it’s constant.

Bonus note: Different types of algorithms,like symplectic, evaluate the velocity at different time points.

### A cell of parameters we need to enter

This first cell is parameters that we need to decided upon. I’ve put in numbers that are fairly resonable and tested.

You can change the value of them and see what happens to the final answer. Hopefully nothing will break ;)

```
θ=π/4 #angle with respect to horizontal
v0=10 # intial velocity
x0=0; # initial position
y0=0;
t0=0 #initial time
tf=2 #final time
dt=1e-3 #time step size
```

### Preliminary calculation of important variables

This is again more parameters, but this cell shouldn’t be changed.

Unless you want to go to the moon or Mars…

```
g=9.8 #m/s^2
vx0=v0*cos(θ) # x component of velocity at initial time step
vy0=v0*sin(θ) # y component of velocity at initial time step
t=collect(t0:dt:tf) # creates an array that holds each time point
nt=length(t) #the number of time steps we willl have
```

### Initialization of Variables

As we march along for a bunch of time steps, we will be computing our position and velocity, but we will need some place to put those numbers. It’s more efficient and better practice to create a place to put those numbers beforehand. So that’s what we are doing here.

We are also setting the first value in the position/ velocity vectors to be their initial values.

```
#initializing empty vectors that will hold position and velocity at each time step
x=zeros(t)
y=zeros(t)
vx=zeros(t)
vy=zeros(t)
x[1]=x0
y[1]=y0
vx[1]=vx0
vy[1]=vy0
```

### Loop of Time Steps

Now we do our actual calculation. We march along, taking tiny little baby steps, calculating our new positions and our new velocities and storing them in our vectors.

```
for ii in 2:nt
x[ii]=x[ii-1]+vx[ii-1]*dt
y[ii]=y[ii-1]+vy[ii-1]*dt
vx[ii]=vx[ii-1]
vy[ii]=vy[ii-1]-g*dt
end
```

## Plotting

We use PyPlot to display our results here. It’s imported into Julia from Python and has great documentation and versitility.

### Tips from Expierence

Always include x and y labels, title, legends, and relevant units **on the graph**.

The graph might seem obvious to you now, but the labeling might not seem obvious to you next week, next month, or next year. And it probably won't seem obvious to someone else looking at your work. So I'll be taking a few extra lines to make sure I include all that.

```
plot(x,y)
xlabel("x")
ylabel("y")
title("Simple Projectile Motion")
```

```
yexact=-g/2.*t.^2.+vy0.*t
plot(t,y,linewidth=6,label="Numerical")
plot(t,yexact,label="Exact")
xlabel("t")
ylabel("y")
legend()
title("Comparing Numerical and Exact Calculations")
```

```
plot(t,yexact-y)
xlabel("t")
ylabel("Y Error")
title("Error over time of Numerical solution at dt= $dt")
```

```
plot(t,vy,label="vy")
plot(t,vx,label="vx")
xlabel("t")
ylabel("velocity")
legend()
title("Velocities over time")
```

## But what if there is a surface?

If we place a surface at $y=0$, or any other location, the ball won’t just keep on falling forever. We can choose two types of actions when it encounters the surface

**Elastic Collision:**The ball bounces with the same amount of kinetic energy, just in the opposite direction**Inelastic Collision:**The ball looses a fraction of its energy in the collision.- At an extreme of this case, the ball looses all its energy.

We use an `if`

statement to determine when it encounters the surface. We’ll just do an elastic collision, so we can just change $v_y$ and not have to worry about total energy, etc.

```
for ii in 2:nt
x[ii]=x[ii-1]+vx[ii-1]*dt
y[ii]=y[ii-1]+vy[ii-1]*dt
if y[ii]<0
println("Hit the surface at t: ",t[ii],"\t x: ",x[ii])
vy[ii-1]=-vy[ii-1]
y[ii]=y[ii-1]+vy[ii-1]*dt
end
vx[ii]=vx[ii-1]
vy[ii]=vy[ii-1]-g*dt
end
```

```
Hit the surface at t: 1.445 x: 10.217692988145673
```

```
plot(x,y)
xlabel("x")
ylabel("y")
title("Projectile Motion with a Surface")
```

```
plot(t,vy,label="vy")
plot(t,vx,label="vx")
xlabel("t")
ylabel("velocity")
legend()
title("Velocity over time with bounce")
```

## Why do we do analytics at all?

So we just found our intercepts and made a bunch of nice graphs without ever having to do any algebra. Why do we force you to slog through it?

Because you need to know what to expect to know when it’s completely going wrong or when something might fail numerically.

This problem is pretty robust, but if you try other problems like the Van der Pol problem, or one’s that display any instability or chaos, we would really need to callibrate to the right time step size. More specialized algorithms, like the ones I talk about in my ODE post, also allow you more accuracy for larger step sizes.

Let’s do the same things but for a range of different time steps to show some of the limitations of numerical analysis.

```
# Lets choose our step sizes
dta=[.001,.01,.1,.2]
```

```
4-element Array{Float64,1}:
0.001
0.01
0.1
0.2
```

```
# The length of each time series
# All arrays will be the same length, but some will be padded
na=floor(Int,(tf-t0)./dta)
# The place holding arrays
ta=zeros(Float64,maximum(na),length(dta))
xa=zeros(ta)
ya=zeros(ta)
vxa=zeros(ta)
vya=zeros(ta)
# Where we start
xa[1,:]=x0
ya[1,:]=y0
vxa[1,:]=vx0
vya[1,:]=vy0
# We perform one loop for each step size
for jj in 1:length(dta)
# This is the step size for the loop
dt=dta[jj]
# We only use the first na[jj] of the arrays
ta[1:na[jj],jj]=linspace(t0,tf,na[jj])
#The same loop we had before
for ii in 2:na[jj]
xa[ii,jj]=xa[ii-1,jj]+vxa[ii-1,jj]*dt
ya[ii,jj]=ya[ii-1,jj]+vya[ii-1,jj]*dt
if ya[ii,jj]<0
vya[ii-1,jj]=-vya[ii-1,jj]
ya[ii,jj]=ya[ii-1,jj]+vya[ii-1,jj]*dt
end
vxa[ii,jj]=vxa[ii-1,jj]
vya[ii,jj]=vya[ii-1,jj]-g*dt
end
end
```

```
for ii in 1:length(dta)
plot(xa[1:na[ii],ii],ya[1:na[ii],ii],label=dta[ii])
end
xlabel("x")
ylabel("y")
legend()
title("Projectile Motions for different time step sizes")
```

```
yexact=-g/2.*ta.^2.+vy0.*ta
yerr=ya-yexact
ylim(-.5,4)
for ii in 1:length(dta)
plot(ta[1:na[ii],ii],yerr[1:na[ii],ii],label=dta[ii])
end
annotate("The Bounce",
xy=[1.5, 1.4 ],
xytext=[.5, 1.5],
xycoords="data",
arrowprops=Dict("facecolor"=>"green"))
xlabel("time")
ylabel("y error")
legend()
title("Error in y position for different time step sizes")
```

## But what about Air Resistence?

Real objects encounter air resistence proportional to velocity. This can’t be solved analytically, but our code can handle it really easily.

```
R=.01
dt=1e-3
```

```
t=collect(t0:dt:tf) # creates an array that holds each time point
nt=length(t) #the number of time steps we willl have
#initializing empty vectors that will hold position and velocity at each time step
x=zeros(t)
y=zeros(t)
vx=zeros(t)
vy=zeros(t)
x[1]=x0
y[1]=y0
vx[1]=vx0
vy[1]=vy0
E=zeros(t)
E[1]=.5*(vx0^2+vy0^2)
```

```
for ii in 2:nt
x[ii]=x[ii-1]+vx[ii-1]*dt
y[ii]=y[ii-1]+vy[ii-1]*dt
if y[ii]<0
#println("Hit the surface at t: ",t[ii],"\t x: ",x[ii])
vy[ii-1]=-vy[ii-1]
y[ii]=y[ii-1]+vy[ii-1]*dt
end
sc=sqrt(vx[ii-1]^4+vy[ii-1]^4)/sqrt(vx[ii-1]^2+vy[ii-1]^2)
vx[ii]=vx[ii-1]-R*vx[ii-1]*sc*dt
vy[ii]=vy[ii-1]-g*dt-R*vy[ii-1]*sc*dt
E[ii]=g*y[ii]+.5*(vx[ii]^2+vy[ii]^2)
end
```

```
plot(x,y)
```

```
yexact=-g/2.*t.^2.+vy0.*t
plot(t,y,label="Air Resistence")
plot(t,yexact,label="Exact No Air")
ylim(-.5,3)
xlabel("t")
ylabel("y")
legend()
title("Comparing With and Without Air")
```

```
plot(t,vy)
```

```
plot(t,vx)
```

```
plot(t,E)
xlabel("time")
ylabel("Energy")
title("Decaying Energy with Air Resistence")
```

# Terminal Velocity

Since we already have the piece of code for air resistance written, let’s just run it for a freely falling object to determine terminal velocity.

```
tf=10
```

```
t=collect(t0:dt:tf) # creates an array that holds each time point
nt=length(t) #the number of time steps we will have
#initializing empty vectors that will hold position and velocity at each time step
y=zeros(t)
vy=zeros(t)
E=zeros(t)
y[1]=0
vy[1]=0
E[1]=.5*vy0^2
```

```
for ii in 2:nt
y[ii]=y[ii-1]+vy[ii-1]*dt
vy[ii]=vy[ii-1]-g*dt-R*sign(vy[ii-1])*vy[ii-1]^2*dt
E[ii]=g*y[ii]+.5*vy[ii]^2
end
```

```
plot(t,vy,label="Velocity")
plot(t,vy[end]*ones(t),label=round(vy[end],2))
xlabel("time")
ylabel("y velocity")
legend()
title("Approaching terminal velocity")
```

```
plot(t,y)
xlabel("time")
ylabel("Position")
title("Particle Falling in Air Resistence")
```

```
plot(t,E)
xlabel("time")
ylabel("Energy")
title("Energy of Falling Particle")
```

## The Termination

Now that we have see terminal velocity, I will terminate this post. See you next time :)