Kernel: SageMath 9.7

C37 Pendulum

Assumptions

  1. F=ma\overrightarrow{F} = m \overrightarrow{a} (Newton's seonc law of motion)

  2. F=Fgravity+Frod+Fdrag+Fexternal\overrightarrow{F} = \overrightarrow{F}_{\text{gravity}} + \overrightarrow{F}_{\text{rod}} + \overrightarrow{F}_{\text{drag}} + \overrightarrow{F}_{\text{external}} (forces are additive and these are the only ones in this situation)

  3. Fgravity=mg(0,1)\overrightarrow{F}_{\text{gravity}} = mg(0, -1) (gravity is a constant downward pull)

  4. Frod \overrightarrow{F}_{\text{rod}} is such that the distance between pivot and bob is constant

  5. Fdrag=bv\overrightarrow{F}_{\text{drag}} = -b\overrightarrow{v} (air and pivot resistance)

  6. Fexternal=Fcos(ωt)(0,1)\overrightarrow{F}_{\text{external}} = F \cos(\omega t) (0,1) (action on the pivot)

Model Derivation

Consider the (x,y)(x, y) position of the bob as a function of θ\theta.

position=p=l(sin(θ),cos(θ))velocity=v=lθ.(cos(θ),sin(θ))acceleration=a=lθ..(cos(θ),sin(θ))+l(θ.)2(sin(θ),cos(θ))\qquad \begin{array}{lcccl} \text{position} &=& \overrightarrow{p} &=& l(\sin(\theta), -\cos(\theta)) \\ \text{velocity} &=& \overrightarrow{v} &=& l\overset{.}{\theta}(\cos(\theta), \sin(\theta)) \\ \text{acceleration} &=& \overrightarrow{a} &=& l\overset{..}{\theta}(\cos(\theta), \sin(\theta)) + l(\overset{.}{\theta})^2(-\sin(\theta), \cos(\theta)) \\ \end{array}

By assumptions 1 and 2:

ma=Fgravity+Frod+Fdrag+Fexternal\qquad m \overrightarrow{a} = \overrightarrow{F}_{\text{gravity}} + \overrightarrow{F}_{\text{rod}} + \overrightarrow{F}_{\text{drag}} + \overrightarrow{F}_{\text{external}}

In the direction of motion

mlθ..=mgsin(θ)+0blθ.+Fcos(ωt)sin(θ)\qquad m l \overset{..}{\theta} = -m g \sin(\theta) + 0 - b l \overset{.}{\theta} + F \cos(\omega t) \sin(\theta)

As a system,

θ.=vv.=glsin(θ)bmv+Fmlcos(ωt)sin(θ)\qquad \begin{array}{lcl} \overset{.}{\theta} &=& v \\ \overset{.}{v} &=& -\dfrac{g}{l} \sin(\theta) - \dfrac{b}{m} v + \dfrac{F}{ml} \cos(\omega t) \sin(\theta) \end{array}

Choose units so that g=l=m=1g = l = m = 1:

θ.=vv.=sin(θ)bv+Fsin(θ)cos(ωt)\qquad \begin{array}{lcl} \overset{.}{\theta} &=& v \\ \overset{.}{v} &=& -\sin(\theta) - b v + F \sin(\theta) \cos(\omega t) \end{array}

Phase Diagrams

var('u v t'); l = 5 b = 0; F = 0; w = 1; u0 = 1; tmax = 30 udot = v vdot = -sin(u) - b*v terr = contour_plot(v^2/2 - cos(u), (u, -l, l), (v, -l, l)) ics = [(0.5,0), (1,0), (2,0), (3.1,0), (0,3), (0,4), (0,-3), (0,-4)] solns = streamline_plot([udot, vdot], (u, -l, l), (v, -l, l), start_points = ics, color = 'green') unull = line([(-l,0), (l,0)], color = 'red') vnull = implicit_plot(vdot == 0, (u, -l, l), (v, -l, l), color = 'blue') show(terr + solns + unull + vnull, axes = False)

Numerical Solutions

To obtain a numerical solution, the parameters, initial conditions, and the time domain must be specified numerically. In the following, I use uu for θ\theta and ww for ω\omega.

Start with an unforced and no drag situation, and consider different initial conditions u(0)=u0u(0) = u_0 and v(0)=0v(0) = 0.

var('u v t') b = 0; F = 0; w = 1; u0 = 1; tmax = 30 udot = v vdot = -sin(u) - b*v + F*cos(w*t)*sin(u) points = desolve_system_rk4([udot, vdot], [u, v], ivar=t, ics=[0,u0,0], step=0.1, end_points=tmax)
uplot = list_plot([[td,ud] for td,ud,vd in points], color='red') vplot = list_plot([[td,vd] for td,ud,vd in points], color='blue') uplot + vplot

What is the period?

udata = [ud for td,ud,vd in points] udata.index(max(udata[1:]))
list_plot([[ud,vd] for td,ud,vd in points])

Now add some forcing at a period near the natural period.

var('u v t') b = 0; F = 0.1; w = 2*pi/6.7; u0 = 1; tmax = 500 udot = v vdot = -sin(u) - b*v + F*cos(w*t)*sin(u) points = desolve_system_rk4([udot, vdot], [u, v], ivar=t, ics=[0,u0,0], step=0.1, end_points=tmax) list_plot([[ud,vd] for td,ud,vd in points])