Computational Dynamics -- Examples¶
This notebook provides a few examples of the techniques discussed in "Introduction to Computational Dynamics".
from numpy import *
import matplotlib.pyplot as plt
A. Two-dimensional Spring¶
Suppose a particle in a plane is subjected to a linear restoring force in $x$ and in $y$ -- this would be a sort of two-dimensional spring. Here $k_x$ and $k_y$ are the spring constants (the stiffness or strength of the spring) in each direction. The resulting path of motion is known as a Lissajous figure (see one being drawn by a pendulum here).
def spring_force_2D(k,r):
"""Given the spring constants in the x and y directions (as a 2-element vector)
and the current position, this function returns the (2D) force on a
particle by a two-dimensional spring.
"""
force = -k * r # multiplies arrays component-wise
return force
t0,tf = 0, 32.0 # initial and final time (seconds)
dt = 0.001 # delta t (seconds)
r0 = array(( 1, 0 )) # initial position vector (m)
v0 = array(( 1, 1 )) # initial velocity (m/s)
kx,ky = 1,0.64 # spring constants, (N/m)
k = array(( kx, ky ))
m = 1.0 # mass of body (kg)
r,v,t = r0,v0,t0 # set initial position, velocity and time
rr,vv,tt = [r],[v],[t] # start lists
i = 0 # index starts at 0
while t <= tf: # keep repeating until time tf
F = spring_force_2D(k,r) # get current force
a = F/m # current acceleration
t = t + dt # increment the time
r = r + v*dt # increment the position
v = v + a*dt # increment the velocity
rr.append(r) # append r,v,t to lists
vv.append(v) #
tt.append(t) #
rr = array(rr) # for convenience later, convert the
vv = array(vv) # lists into arrays
tt = array(tt) #
# plot results:
fig,ax = plt.subplots()
ax.grid(True)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('2D Spring')
ax.set_aspect('equal')
ax.scatter( rr[:,0], rr[:,1], c=tt, ) # color the dots by time
#plt.savefig('2d-spring.png')
plt.show()
You can make some pretty figures by playing around with the spring constants in the code above (especially when $k_x = c^2 k_y$, where $c$ is a ratio of integers).
B. Gravitational orbits¶
Newton's Law of Universal Gravitation is
$$ \vec{F}_G = G \frac{m_1 m_2}{r^2} \hat{r} $$
The force of gravity produces a circular orbit if the velocity is just right.
def gravity_force(k,m,r):
"""Given a constant k (which could be G*M for Newtonian gravity)
and the current position, this function returns an inverse-square
force vector directed radially inward.
"""
r_mag = sqrt(sum(r**2)) # magnitude of r vector
rhat = r/r_mag # radial unit vector
force = -k*m / r_mag**2 * rhat
return force
k = 1.0 # constant gives overall gravity strength
m = 1.0 # mass of orbiter (doesn't actually matter)
t0,tf = 0, 10 # initial and final time (seconds)
dt = 0.001 # delta t (seconds)
r0 = array(( 1, 0 )) # initial position vector (m)
v0 = array(( 0, 1.0 )) # initial velocity (m/s)
r,v,t = r0,v0,t0 # set initial position, velocity and time
rr,vv,tt = [r],[v],[t] # start lists
i = 0 # index starts at 0
while t <= tf: # keep repeating until time tf
F = gravity_force(k,m,r) # current force vector
a = F/m # current acceleration
t = t + dt # increment the time
r = r + v*dt # increment the position
v = v + a*dt # increment the velocity
rr.append(r) # append r,v,t to lists
vv.append(v) #
tt.append(t) #
rr = array(rr) # for convenience later, convert the
vv = array(vv) # lists into arrays
tt = array(tt) #
# plot results:
fig,ax = plt.subplots()
ax.grid(True)
ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_title('Gravitational Orbit')
ax.set_aspect('equal')
ax.scatter( rr[:,0], rr[:,1], c=tt, ) # color the dots by time
ax.plot( [0,],[0,],'o',color='orange') # show a dot at central mass
#plt.savefig('orbit.png')
plt.show()
Try changing the initial velocity to see elliptical orbits. You can even find escape velocity -- where the object does not return (if $k=1$ it will be a little over 1.4). For greater speed you may need a longer time domain.