Sage Notebook for the Notes on "Theoretical Mechanics"

Jürgen Vollmer (Uni Leipzig, summer 2020)

License: Creative Commons Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)

exploring the solutinons of a 1D EOM

Setup

base-name for output files

i.e. the common part of all output files

In [1]:
baseName = "./EOM_1D_explorer_"

plot layout

In [2]:
xMin = -2
xMax =  2
vMin = -2
vMax =  2

plot_params = dict(xmin = xMin,  xmax = xMax,   # interval of x-axis
                   ymin = vMin,  ymax = vMax,   #             y-axis
                   aspect_ratio=1,              # equal scaling of x and y-axis (circles look like circles!)
                                     
                   ticks=(1,1),                 # ticks with separation of $(\Delta x, \Delta y)$
                   axes_labels=(r"$x$",       # axis labels for  x-axis
                                r"$\dot x$"),  # axis label  for  y-axis

                   fontsize=16,                 # size of text font

                   figsize=(4,4)                # figure size
                  )

update axis labels

In [3]:
plot_params.update( axes_labels = (r"$h$",        # axis labels for  x-axis
                                   r"$\dot h$")   # axis label  for  y-axis
                  )

EOM

EOM in function notation

In [4]:
# declare variables 
x,v = var("x,v")

# define function
def EOM (x,v) :
     return (v, -1-v)

EOM in vector notation X = (position, velocity)

In [5]:
def dGamma_dt(X, t=0): 
    return [X[1], -1-X[1]]

## solution will be evaluated for times
t = srange(0, 40, 0.01)

## numerical solution of the ODE
from scipy import integrate
def trajectory( IC, **plt_para ) :
  TrajectoryPts = integrate.odeint( dGamma_dt, IC, t )
  return line( TrajectoryPts, **plt_para )

Numerical solution of trajectories

In [6]:
##   plot a single solution
## make + show plot

IC = (1,1)
trajectory( IC, color="blue", thickness=2, **plot_params )
Out[6]:

What about

  • other initial positions (i.e. height)?
  • other initial velocities?

Where is time here?

In [7]:
IC1 = (1,0)
IC2 = (0.55, 1)
plt  = trajectory( IC1, color="blue", thickness=2, **plot_params )
plt += trajectory( IC2, color="green", thickness=2, **plot_params )

show(plt)
In [8]:
##   plot many solutions
## make + show plot

pt = plot( 0 ) 

for xIn in srange(xMin-2, xMax, 1) :
    pt += trajectory( (xIn,  vMax), color="green", thickness=2, **plot_params )

for xIn in srange(xMin, xMax+3, 1) :
    pt += trajectory( (xIn,  vMin), color="blue", thickness=2, **plot_params )

pt
Out[8]:

Vector field showing the flow direction and speed

In [9]:
## add plot parameter for arrow shape and positions
plot_params4ps = dict(
    # shaftwidth=3,
    headwidth=10,
    plot_points=10,
    **plot_params )

## make + show plot
pv = plot_vector_field(EOM(x,v), (x,xMin,xMax), (v,vMin,vMax), color="lightgray", **plot_params4ps )
pv
Out[9]:

Merge plots

In [12]:
p = pt+pv+plot(-1, (-2,2), thickness=3, color="red")
p
Out[12]:

save plot

In [14]:
# p.save_image( baseName+"phaseSpace-h.svg" )
# p.save_image( baseName+"phaseSpace-h.pdf" )
p.save_image( baseName+"phaseSpace-h.png", transparent=True)
In [ ]: