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}For the harmonic oscillator in this problem there is noticeable growth in amplitude, demonstrating instability in Euler's method, even for $\Delta t = 0.001$.
September 2019
Marty Ligare
In [2]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
In [3]:
# Following is an Ipython magic command that puts figures in the notebook.
%matplotlib notebook
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
mpl.style.use('classic')
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 [4]:
# Linear Restoring Force
def f(x):
return -k*x
Set constants and initialize arrays¶
In [5]:
m = 1 # mass
k = np.pi**2 # spring constant
x0 = 1 # initial position
v0 = 0 # initial velocity
t0 = 0 # initial time
tf = 10 # final time
n = 10001 # number of points
dt = (tf-t0)/(n-1) # time step
t = np.linspace(0, tf, n) # Create array of t-values
x = np.zeros(len(t)) # Create array for values of x
v = np.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 [6]:
for i in range(1,n-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 [7]:
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 [8]:
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 for more information and instructions for package installation.
version_information is installed on the linux network at Bucknell
In [9]:
%load_ext version_information
In [10]:
version_information numpy, matplotlib
Out[10]:
In [ ]: