## Numerical integration of Newton's second law¶

Repeat of kind of thing that is done in Excel in PHYS 211 lab, and in Mathematica in PHYS 221 lab.

The Euler algorithm is

\begin{eqnarray} t_\text{new} &=& t_\text{old} + \Delta t \\ x_\text{new} &=& x_\text{old} + v_\text{old}\, \Delta t \\ v_\text{new} &=& v_\text{old} + a_\text{old}\, \Delta t \end{eqnarray}

December 2016 Marty Ligare

In [1]:
import scipy as sp

import matplotlib as mpl       # As of July 2017 Bucknell computers use v. 2.x
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in the  notebook.
# For figures in separate windows, comment out following line and uncomment
# the next line
# Must come before defaults are changed.
%matplotlib notebook
#%matplotlib

# As of Aug. 2017 reverting to 1.x defaults.
# In 2.x text.ustex requires dvipng, texlive-latex-extra, and texlive-fonts-recommended,
# which don't seem to be universal
# See https://stackoverflow.com/questions/38906356/error-running-matplotlib-in-latex-type1cm?
mpl.style.use('classic')

# M.L. modifications of matplotlib defaults using syntax of v.2.0
# Changes can also be put in matplotlibrc file, or effected using mpl.rcParams[]
plt.rc('figure', figsize = (6, 4.5))            # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True)             # Adjusts supblot parameters for new size


#### Define functions¶

In [2]:
# Linear Restoring Force
def f(x):
return -k*x


#### Set constants and initialize arrays¶

In [3]:
m = 1                         # mass
k = sp.pi**2                  # spring constant
x0 = 1                        # initial position
v0 = 0                        # initial velocity
tf = 10                       # final time
dt = 0.001                    # time step
np = int(tf/dt)               # Number of points (number of intervals = np - 1)

t = sp.linspace(0, tf, np)    # Create array of t-values
x = sp.zeros(len(t))          # Create array for values of x
v = sp.zeros(len(t))          # Create array for values of v
x[0] = x0                     # Fix initial value of velocity
v[0] = v0                     # Fix initial value of velocity


#### Integrate¶

In [4]:
for i in range(1,np-1):
x[i] = x[i-1] + v[i-1]*dt        # x_new = x_old + v_old*dt
v[i] = v[i-1] + f(x[i-1])*dt/m   # v_new = v_old + a_old*dt

In [5]:
plt.figure()
plt.xlabel('$t$')      # Label for horizontal axis
plt.ylabel("$x$")      # Label for vertical axis
plt.title("Linear Restoring Force",fontsize=14)
plt.grid(True)
plt.axhline(0,color='green')  # Makes solid green x-axis
plt.plot(t,x);

In [6]:
plt.figure()
plt.xlabel('$t$')      # Label for horizontal axis
plt.ylabel("$v$")      # Label for vertical axis
plt.title("Linear Restoring Force",fontsize=14)
plt.grid(True)
plt.axhline(0,color='green')  # Makes solid green x-axis
plt.plot(t,v);


### Version Information¶

version_information is from J.R. Johansson (jrjohansson at gmail.com)
See Introduction to scientific computing with Python:
http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-0-Scientific-Computing-with-Python.ipynb

If version_information has been installed system wide (as it has been on linuxremotes), continue with next cell as written. If not, comment out top line in next cell and uncomment the second line.

In [7]:
%load_ext version_information

#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py

Loading extensions from ~/.ipython/extensions is deprecated. We recommend managing extensions like any other Python packages, in site-packages.

In [8]:
version_information scipy, matplotlib

Out[8]:
SoftwareVersion
Python3.6.1 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]
IPython6.1.0
OSLinux 3.10.0 327.36.3.el7.x86_64 x86_64 with redhat 7.2 Maipo
scipy0.19.1
matplotlib2.0.2
Tue Aug 01 11:13:59 2017 EDT
In [ ]: