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)

Phase Portrait of the Mathematical Pendulum

Setup

base-name for output files

i.e. the common part of all output files

baseName = "./EOM_mathPendulum__"

plot layout

plot_params = dict(xmin = -4,    xmax = 4,      # interval of x-axis
                   ymin = -3,    ymax = 3,      #             y-axis
                   aspect_ratio=1,              # equal scaling of x and y-axis (circles look like circles!)
                                     
                   ticks=(pi,2),                 # ticks with separation of $(\Delta x, \Delta y)$
                   tick_formatter=pi,
                   axes_labels=(r"$\theta$",         # axis labels for  x-axis
                                r"$\dot \theta/\omega$"),   # axis label  for  y-axis

                   fontsize=16,                 # size of text font

                   figsize=(4,4)                # figure size
                  )

EOM and energy

as function of the variable

theta,v = var("theta,v")

EOM

def EOM (theta,v) :
    return (v, -sin(theta))

energy

def E (theta, v) :
    return -cos(theta) + v**2/2

Plot of energy function with horizontal plane at height zPlane

zPlane = 1.5

## add plot parameters for 3d plot
plot_params3d     = dict( plot_params )
param_update = dict(axes_labels = ("theta","v","E"),
                    opacity=0.87
                )
plot_params3d.update( param_update )

## make + show plot
p3d  = plot3d(E(x,v), (x,-4, 4), (v,-3,3), **plot_params3d )
p3d += plot3d(zPlane, (x,-4, 4), (v,-3,3), color="orange", **plot_params3d )
p3d

save plot

p3d.save_image( baseName+"energy-surface.png" )

interactive version to move plane

@interact
def _(zPlane=(0,2,0.1)):
   ## add plot parameters for 3d plot
   plot_params3d     = dict( plot_params )
   param_update = dict(axes_labels = ("x","v","E"),
                       opacity=0.87
                       )
   plot_params3d.update( param_update )

   ## make + show plot
   p3d  = plot3d(E(x,v), (x,-4, 4), (v,-2.5,2.5), **plot_params3d )
   p3d += plot3d(zPlane, (x,-4, 4), (v,-2.5,2.5), color="orange", **plot_params3d )
   p3d.show()

Widget Javascript not detected.  It may not be installed or enabled properly.

Potential

## add plot parameters for 3d plot
plot_params_Pot     = dict(xmin = -4,    xmax = 4,      # interval of x-axis
                   ymin = -1.5,    ymax = 1.5,      #             y-axis
                   aspect_ratio=1,              # equal scaling of x and y-axis (circles look like circles!)
                                     
                   ticks=(pi,1),                 # ticks with separation of $(\Delta x, \Delta y)$
                   tick_formatter=pi,
                   axes_labels=(r"",         # axis labels for  x-axis
                                r"$U$"),   # axis label  for  y-axis

                   fontsize=16,                 # size of text font

                   figsize=(4,4),                # figure size
                   thickness = 4
                 )

ppot = plot( -cos(x), (x,-4,4), **plot_params_Pot )

show(ppot)

Contour lines for const energy

## add plot parameters for color bar
plot_params4pc = dict(plot_params)
param_update = dict(
    ymin = -3, 
    ymax = 5.5,
    axes_labels=(r"$\theta$",         # axis labels for  x-axis
                 r"$\dot \theta/\omega$       "),   # axis label  for  y-axis
    ticks=([-pi,0,pi],[-2,0,2, 3.5,4.5,5.5]),                 # ticks with separation of $(\Delta x, \Delta y)$
    tick_formatter=[(r"$-\pi$", r"$0$", r"$\pi$"), ("-2","0","2",r"$-MgL$","",r"$MgL$")],
    colorbar=False, 
    contours=[-0.995, -0.5, 0, .5, 1, 2, 4],
    linewidths = 2.5,
    fill=False,
    labels=True,
    axes = True,
    # label_color='black',
    cmap='hsv')
plot_params4pc.update(param_update)

## make + show plot
pc = contour_plot(E(x,v), (x,-4, 4), (v,-3,3.3), **plot_params4pc )

pc += plot(4.5-cos(x), (x,-4,4), thickness=4)
pc += plot((3.3,4.5), (x,-4,4), color="black" )
pc += text( r"$U$", (-0.5,0.9), fontsize=22, rotation=90, zorder=0, axis_coords=True, color="black")

pc

Flow of EOM

ps = streamline_plot(EOM(x,v), (x,-4,4), (v,-3,3), **plot_params )
ps

Vector field showing the flow direction and speed

## add plot parameter for arrow shape and positions
plot_params4ps = dict(
    # shaftwidth=3,
    headwidth=5,
    plot_points=10,
    **plot_params )

## make + show plot
pv = plot_vector_field(EOM(x,v), (x,-4,4), (v,-3,3), **plot_params4ps )
pv.show()

## make + show plot with gray arrows
pvg = plot_vector_field(EOM(x,v), (x,-4,4), (v,-3,3), color="gray", **plot_params4ps )

Merge plots

p = pvg+pc
p

save plot

p.save_image( baseName+"phaseSpace.svg" )
# p.save_image( baseName+"phaseSpace.pdf" )

homoclines

integral(1/cos(theta/2), theta).show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\log\left(\sin\left(\frac{1}{2} \, \theta\right) + 1\right) - \log\left(\sin\left(\frac{1}{2} \, \theta\right) - 1\right)$$

integral(1/sqrt(2+2*cos(theta)), theta).simplify_full().show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \log\left(2 \, \sin\left(\frac{1}{2} \, \theta\right) + 2\right) - \frac{1}{2} \, \log\left(-2 \, \sin\left(\frac{1}{2} \, \theta\right) + 2\right)$$

integral(1/sqrt(1+cos(theta)), theta).simplify_full().show()

$$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, \sqrt{2} \log\left(2 \, \sin\left(\frac{1}{2} \, \theta\right) + 2\right) - \frac{1}{2} \, \sqrt{2} \log\left(-2 \, \sin\left(\frac{1}{2} \, \theta\right) + 2\right)$$

t = var("t")
plot_params = dict(ticks=( [-10,-5,0,5,10], [-pi,0,pi] ),                 # ticks with separation of $(\Delta x, \Delta y)$
                   tick_formatter=[(r"$-10$", r"$-5$", r"$0$", r"$5$", r"$10$"), (r"$-\pi$", r"$0$", r"$\pi$")],
                   axes_labels=(r"$\omega\, (t-t_0)$",         # axis labels for  x-axis
                                r"$\theta(t)$"),   # axis label  for  y-axis

                   fontsize=16,                 # size of text font

                   figsize=(4,3),                # figure size
                   frame=True,
                   axes=True
                  )

Ptheta  = plot( 2*arcsin(tanh(t)), (t,-10,10), thickness=3, **plot_params)
Ptheta += plot( (-pi, pi), (t,-10,10), color="gray", linestyle="dotted",
              ymin=-3.5, ymax=3.5)

show(Ptheta)

save plot

Ptheta.save_image( baseName+"homocline.svg" )
Ptheta.save_image( baseName+"homocline.pdf" )