Solving ordinary differential equations exactly can be very challenging. In this project we were asked to examine the behavior of a pendulum system over time, using the trapezoid method as a means of approximating the solution at each time step. The trapezoid method is an improvement on Euler's method, and details on it can be found here: Trapezoid Notes. For an undamped pendulum, the first order differential equation is given by: \begin{align} y_1' &= y_2\\ y_2' &= -\frac{g}{l} \sin y_1 \end{align}
In Part 1, we were just asked to modify the given pend.m file to allow for damped motion. The equation governing this motion is given by: \begin{align} y_1' &= y_2\\ y_2' &= -\frac{g}{l} \sin y_1 - dy_2 \end{align} Unlike an actual pendulum which would very easily fall from an initial condition of $y_1(0)=\pi$, $y_2(0)=0$, a pendulum started in this position in this program will stay there for all time. This is not a physical solution. The modified file can be found here: pend.m
For Part 2, we needed to build a forced damped oscillator, for which the equation is given by: \begin{align} y_1' &= y_2\\ y_2' &= -\frac{g}{l} \sin y_1 - dy_2 + A\sin t \end{align} Here, our damping was d=1 and the forcing parameter was A=10. Our task was to see what would happen using different step sizes in the trapezoid method. The larger the step size, the larger the error in the approximate solution, and hence the behavior exhibited with a step size of h=0.005 was much different than with a step size of h=0.1. With a step size of h=0.005, all solutions eventually settled into some periodic motion, but for h=0.1, the behavior was much more erratic.
For Part 3 we had to set our forcing parameter A=10 and our step size h=0.005. We were told that there were two periodic attractors, and our job was to find two initial conditions $(y_1,y_2)=(a,0)$ and $(b,0)$, where $|a-b|\leq 0.1$ that are attracted to different periodic trajectories. We found that $a=\pi - 0.05$ and $b=\pi - 0.05$ resulted in trajectories that were essentially mirror images of each other. Each trajectory took too long to settle down, and the resulting video files were too large.
Next, we adapted pend.m to build a damped pendulum with oscillating pivot, the
trajectory of which is governed by the ODE:
\begin{align}
y_1' &= y_2\\
y_2' &= -(\frac{g}{l} + A\cos2\pi t)\sin y_1 - dy_2
\end{align}
Setting d=0.1 and the length l=2.5m, we were to find as accurately as possible the
of the forcing strength A such that the inverted pendulum becomes stable. This was a
more difficult problem, and at first we tried values by hand and found A=[19,25] was
one possible range. To try to get more accurate boundaries, we created the code
eqlbrm.m that runs with a given A value, and then increases or
decreases A (depending on the boundary and whether stable or unstable), to try to
converge on the boundary. Using this method, we narrowed the forcing constant to the
interval A=[18.365, 25.52], but outside of the upper bound were other values that
exhibited stable behavior for the inverted pendulum.
To try to explore the boundaries in a different way, we used the loop
p3_A_loop.txt along with our ydot function:
ydot.m. In the loop we use ode45 to solve ydot, and then check
for the first A value such that the minimum $\theta$ value is zero. Here we got the
interval A=[18.356, 25.526], and it took a little while to realize the intervals did
not match because ode45 uses more advanced methods than the trapezoid method to solve
the ode, and hence gets a more accurate result. We did plot the motion given by ode45
in the file exctish.m.
Going back through these two parts, I was unable to duplicate the non-matching results
that I got from the two methods the first time. The values given by the ode45 method
appear to be valid.
Much like in Part 6, we now had o find the minimum A value such that $\theta=0$
becomes an unstable equilibrium. Again, the two methods gave two different results,
and with eqlbrm.m we found A=14.7, and with the loop/ode45 method we found A=14.615.
We were then told to build the double pendulum, for which the motion is governed by the equations: \begin{align} y_1' &= y_2\\ y_2' &= \frac{-3g\sin y_1-g\sin (y_1-2y_3)-2\sin (y_1-y_3)(y^2_4+y^2_2\cos (y_1-y_3))}{3-\cos (2y_1 -2y_3)} - dy_2\\ y_3' &= y_4\\ y_4' &= \frac{2\sin (y_1-y_3)[2y^2_2+2g\cos y_1+y^2_4\cos (y_1-y_3)]}{3-\cos (2y_1 -2y_3)}\\ \end{align} Our file can by view here: dblpend.m. We were unable to stabilize the double pendulum.