{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Introduction to Computational Dynamics\n",
"\n",
"In this notebook we'll do some actual computer programming, using functions and\n",
"loops to compute the motion of an object subjected to a force.\n",
"\n",
"If you're new to programming (or just new to Python), you should first review\n",
"the notebook titled \"[Functional Programming](howto-functional_programming.html)\".\n",
"\n",
"## A. Set up\n",
"\n",
"You will need to import the following packages.\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"from numpy import *\n",
"import matplotlib.pyplot as plt"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## B. Theory\n",
"\n",
"A complex physical system may be converted to a computational Physics problem\n",
"when a direct solution is difficult or impossible. With a computational model,\n",
"a computer can work out how the system evolves by adding together many small\n",
"time steps. \n",
"\n",
"To simulate the continuous functions $x(t)$ and $v(t)$,\n",
"we will convert integral equations, such as\n",
"\n",
"$$x(t) = x_0 + \\int_{0}^t v(t) \\; dt$$\n",
"\n",
"$$v(t) = v_0 + \\int_{0}^t a(t) \\; dt$$\n",
"\n",
"into approximate [finite difference](https://en.wikipedia.org/wiki/Finite_difference) \n",
"equations, such as\n",
"\n",
"$$x_{i+1} = x_i + v_i \\Delta t$$\n",
"\n",
"$$v_{i+1} = v_i + a_i \\Delta t$$\n",
"\n",
"where $i$ labels the many small time steps taken by the computer.\n",
"The result is a sequence of $x$ positions \n",
"$\\left \\{ x_0,x_1,x_2,... \\right \\}$ \n",
"at times \n",
"$\\left \\{t_0, t_1,t_2,... \\right \\}$.\n",
"\n",
"At each time step $i$, the current value of the force is found, which \n",
"may depend on the current position and velocity.\n",
"From this, the current value of the acceleration is obtained\n",
"using Newton's Second Law: $a_i = F/m$.\n",
"With acceleration $a_i$ we find the next velocity $v_{i+1}$.\n",
"With the velocity $v_i$ we find the next position $x_{i+1}$.\n",
"This process is then repeated in a loop for as many time steps as we need.\n",
"\n",
"For motion in 2 or 3 dimensions, we repeat this independently for\n",
"\n",
"- $x$, $v_x$, $a_x$ and $F_x$ \n",
"- $y$, $v_y$, $a_y$ and $F_y$\n",
"- $z$, $v_z$, $a_z$ and $F_z$ \n",
"\n",
"In the sections below, we find the position of a body as follows:\n",
"\n",
"1. **Write a force function**: Program a function that returns the force on the\n",
" body. It may depend on the body's location, its mass, and/or on other\n",
" quantities (such as velocity or some material properties).\n",
"2. **Define the time domain**: initial time $t_i$, final time $t_f$, and the\n",
" time interval $\\Delta t$. The time interval should be small enough that the\n",
" velocity of the body does not change significantly over any step from $t_i$\n",
" to $t_{i+1}$.\n",
"3. **Define the initial conditions**: provide values for the initial position\n",
" and velocity of the body ($\\vec{x}_0$ and $\\vec{v}_0$) and the body's mass\n",
" ($m$).\n",
"4. **Perform a computational loop** over many small time steps.\n",
" At each new time, find the new force, then the new acceleration,\n",
" then the new velocity and position.\n",
" Repeat the loop every $\\Delta t$ until you arrive at the final time of interest. \n",
"\n",
"If you save the position and velocity every time step, you may then plot the\n",
"trajectory or otherwise study the motion.\n",
"\n",
"## C. Projectile Motion Example\n",
"\n",
"Our first computational model is the motion of an object of mass $m$ subject to\n",
"a uniform gravity force, $\\vec{F} = -mg \\hat{\\jmath}$.\n",
"\n",
"Let's follow the procedure listed above for a 1 kg body launched from the\n",
"origin with initial speed and direction $v_0=10$ m/s and $\\theta_0=60^\\circ$.\n",
"\n",
"**1. Write a force function.**\n",
"\n",
"This force is very simple since it is the same at every point in space, and\n",
"does not change in time.\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"g = 9.8 # N/kg # define constant gravity field\n",
"jhat = array((0,1)) # define constant vertical unit vector\n",
"def gravity_force(m):\n",
" Fg = -m*g*jhat\n",
" return Fg"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"**2. Define time domain.**\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"t0,tf = 0, 2.0 # initial and final time (seconds)\n",
"dt = 0.01 # delta t (seconds)"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"**3. Define intial conditions.**\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"v0 = 10 # initial speed (m/s)\n",
"Q0 = 60 # initial direction (degrees)\n",
"rpd = pi/180 # radians per degree conversion factor\n",
"m = 1.0 # mass of body (kg)\n",
"r = array(( 0, 0 )) # initial position vector (m)\n",
"v = v0*array(( cos(Q0*rpd), sin(Q0*rpd) )) # initial velocity vector (m/s)\n",
"t = t0 # initial time"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"We want to keep track of the position $\\vec{r}$ and the velocity $\\vec{v}$ at\n",
"all times. One way to do this is with a \n",
"[Python list](https://realpython.com/python-lists-tuples/), which uses the square\n",
"brackets. We'll make a list of vectors, adding a new one at each time step in a\n",
"loop. We'll use `r_list`, `v_list` and `t_list` for these list\n",
"quantities.\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"r_list = [r,] # start a list of position vectors\n",
"v_list = [v,] # start a list of velocity vectors\n",
"t_list = [t,] # start a list of times"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"**4. Perform the loop.**\n",
"\n",
"There are several ways to proceed. Here I've used a `while` loop:\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"while t <= tf: # keep repeating until time tf\n",
" F = gravity_force(m) # get current force \n",
" a = F/m # current acceleration\n",
" t = t + dt # increment the time\n",
" r = r + v*dt # increment the position\n",
" v = v + a*dt # increment the velocity\n",
" r_list.append(r) # append the new r to list r_list\n",
" v_list.append(v) # .. v v_list\n",
" t_list.append(t) # .. t t_list"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"All the data are now stored in lists of vectors, `r_list`, `v_list` and `t_list`. \n",
"\n",
"It is more convenient to use a two-dimensional array for these quantities,\n",
"where each row is a different time, and the columns hold the $x$ and $y$\n",
"components. For example,\n",
"\n",
"$$\\texttt{rr} = \\begin{bmatrix}\n",
" x_0 & y_0 \\\\ x_1 & y_1 \\\\ x_2 & y_2 \\\\ ... & \n",
"\\end{bmatrix}$$\n",
"\n",
"To access the $x$ and $y$ position at time step 5, we use `rr[5,:]` (row 5, the\n",
"`:` means all columns). To get the entire time-history of the $x$ component we\n",
"use `rr[:,0]` (the `:` means all rows, the `0` means column `0` which is the\n",
"$x$ data).\n",
"\n",
"To convert our lists of vectors into two-dimensional arrays:\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"rr = array(r_list) # for convenience later, convert the\n",
"vv = array(v_list) # lists into arrays\n",
"tt = array(t_list) #"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"That's it. Now `rr` and `vv` contain the position and velocity of the object at\n",
"all times in `tt`. It is often helpful to graph the results.\n",
"\n",
"**Graph the result.**\n",
"\n",
"Here is a graph of $y$ versus $x$ position, colored accoring to time $t$. The\n",
"`matplotlib` plot method we use here is `scatter`. You can see its options in\n",
"the [`scatter` documentation](http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.scatter).\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"x = rr[:,0] # take the first column for x\n",
"y = rr[:,1] # take the second column for y\n",
"fig,ax = plt.subplots() # start making a plot\n",
"ax.scatter(x, y, c=tt, ) # color the dots by time\n",
"ax.grid(True) # show a grid on the plot\n",
"ax.set_xlabel('x (m)')\n",
"ax.set_ylabel('y (m)')\n",
"ax.set_title('projectile in a vacuum')\n",
"ax.set_aspect('equal') # make x- and y-distance the same on screen\n",
"plt.show()"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"You should see a perfect parabola. The computer doesn't know where the ground\n",
"is so it keeps accelerating downward until the final time is reached (this is\n",
"fixed below).\n",
"\n",
"## D. A Loop in a Loop\n",
"\n",
"Now that we have the basic code, we can execute it over and over again in a\n",
"variety of circumstances. \n",
"\n",
"For example, perhaps we'd like to compare the behavior of projectiles on\n",
"different planets with varying gravitational constants. On the moon $g=1.6\n",
"\\;\\rm{N}/\\rm{kg}$; on Jupiter $g=23.5 \\;\\rm{N}/\\rm{kg}$. The code above has\n",
"been written a bit more compactly and put in a loop where `g` is given a few\n",
"different values. In this case, we prepare the plot first, then perform the\n",
"work above for each `g`.\n"
],
"metadata": {}
},
{
"cell_type": "code",
"source": [
"def gravity_force(m):\n",
" Fg = -m*g*array((0,1))\n",
" return Fg\n",
"fig,ax = plt.subplots() # start the plot\n",
"ax.grid(True) \n",
"ax.set_xlabel('x (m)') \n",
"ax.set_ylabel('y (m)') \n",
"ax.set_title('projectiles, changing g') \n",
"ax.set_aspect('equal') \n",
"m = 1.0 # mass of body (kg)\n",
"t0,tf = 0, 2.0 # initial and final time (seconds)\n",
"dt = 0.01 # delta t (seconds)\n",
"r0 = array(( 0, 0 )) # initial position vector (m)\n",
"Q0 = 60 # initial direction (degrees)\n",
"v0 = 10*array(( cos(Q0*pi/180), sin(Q0*pi/180) )) # initial velocity (m/s)\n",
"for g in [1.6,9.8,23.5]:\n",
" r,v,t = r0,v0,t0 # set initial position, velocity and time\n",
" r_list,v_list,t_list = [r],[v],[t] # start lists\n",
" while t <= tf: # keep repeating until time tf\n",
" F = gravity_force(m) # get current force \n",
" a = F/m # current acceleration\n",
" t = t + dt # increment the time\n",
" r = r + v*dt # increment the position\n",
" v = v + a*dt # increment the velocity\n",
" r_list.append(r) # append r,v,t to lists\n",
" v_list.append(v) #\n",
" t_list.append(t) #\n",
" if r[1] < 0: break # stop the while-loop if y<0 (hit the ground)\n",
" rr = array(r_list) # for convenience later, convert the\n",
" vv = array(v_list) # lists into arrays\n",
" tt = array(t_list) #\n",
" # print and plot results:\n",
" print(f'projectile with g = {g:.1f} N/kg:')\n",
" print(f' range {rr[-1,0]:0.2f} m, flight time {t:0.2f} s')\n",
" ax.plot( rr[:,0], rr[:,1], label='{:.1f}'.format(g) )\n",
"ax.legend(loc='upper left') # show a legend on the plot\n",
"plt.show()"
],
"metadata": {},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"Note the range is incorrect for the moon projectile ($g=1.6 \\;\\rm{N}/\\rm{kg}$)\n",
"since it never got a chance to land.\n",
"\n",
"## E. Other Examples\n",
"\n",
"With the methods above you can do quite a bit. \n",
"\n",
"In [Dynamics Examples](dynamics-examples.html) there are a couple more demonstrations. \n",
"\n",
"For example, you can make \"Lissajous figures\" with a 2D spring:\n",
"\n",
"![2D spring](images/2d-spring.png)\n",
"\n",
"Or observe elliptical gravitational orbits:\n",
"\n",
"![gravitational orbit](images/orbit.png)\n"
],
"metadata": {}
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 5
}