What is happening to the pendulum as the critical point is approached?
By giving it just the right initial "kick" (value of \(\omega_0\)), we can make the pendulum swing around to \(\theta=\pi\) radians (or \(180^\circ\)) where it "balances" in an inverted position relative to the equilibrium position where \(\theta=0\). This is in fact another equilibrium position but of course it is unstable. If the pendulum is perturbed even the slightest amount when it is in that position it will start to swing back towards \(\theta=0\). Nonetheless, it is quite remarkable we can investigate such an exquisitely sensitive phenomenon with such a simple simulation. You may have noticed that as you tuned to \(\omega_0 = \omega_0^\star\) the plots of \(\theta(t)\) developed regions where the curve was flatter and flatter for longer and longer periods of time. This is the pendulum "lingering" very near the unstable equilibrium for longer and longer. If we were able to tune infinitely closely to \(\omega_0^\star\), which we can't do because the arithmetic in MATLAB is finite precision, we could in principle make the pendulum balance in the inverted position for an infinite amount of time.
You were asked to investigate the dependence of the critical value of the initial angular speed of the pendulum, \(\omega_0^\star\), as a function of the discrete time spacing, \(\Delta t\).
Recall that you first did the search with level = 7, which
means that there were 65 time steps in the simulation. I asked
you to repeat the search for level = 8 and level = 9 and I
have done that, as well as level = 10, with the following
results:
level omega0* d(omega0*) 4^(level - 7) d(omega0*)
7 1.98348777361164 0.0124290 0.0123
8 1.99591675541708 0.0030651 0.0123
9 1.99898186533854 7.63766e-04 0.0122
10 1.99974563183654
where d(omega0*) = omega0*(level+1) - omega0*(level)
There are two things that are apparent from this table
\(\omega_0^\star\) is converging to some value, and, at the same rate, namely \(O(\Delta t^2)\), as the fundamental dependent variable in the simulation, \(\theta(t^n)\).
Especially given how sensitive the simulation is to changes in the initial conditions as \(\omega_0 \to \omega_0^\star\) this is truly remarkable.
It seems natural to conjecture that the continuum ("exact") value of \(\omega_0^\star\) is exactly 2.
First, observe that one can show that the various energy quantities associated with the pendulum are
\[ T = \frac{\omega^2}{2} \]
\[ V = 1 - \cos(\theta) \]
\[ E = T + V \]
for the kinetic, potential and total energy, respectively, and where we recall that we are working in units in which \(m = g = L = 1\).
Thus for the critical case, where we give the pendulum the precise initial kick needed for it to swing around by an angle \(\pi\) and stop there (in principle, forever, even though it is an unstable equilibrium), we have
\[
E_{\rm initial} = T_{\rm initial} + V_{\rm initial}
= \frac{(\omega_0^\star)^2}{2} + 0
= \frac{(\omega_0^\star)^2}{2}
\]
\[
E_{\rm final} = T_{\rm final} + V_{\rm final}
= 0 + (1 - \cos(\pi)) = 2
\]
Now, since \(E_{\rm initial} = E_{\rm final}\) we have
\[ \frac{(\omega_0^\star)^2}{2} = 2 \] or \[ \omega_0^\star = 2 \]
This is a nice illustration of something that sometimes happens in computational science: we tackle a problem which, for some reason or another, we aren't able to get much insight into through "traditional" pen and paper techniques. So we turn to simulation and from the results of the direct numerical calculations we find something "simple" that begs understanding in a more fundamental way. In this case, the appearance of the number 2 in the results is a clear signal that there is probably an "easier" way of getting that answer than through the detailed computations that we made.
In the 1920's Van der Pol experimented with oscillations in a vacuum tube triode circuit (an analogue electronic circuit), and found interesting behaviour that you will reproduce which is absent from a linear circuit. To model this phenomena he proposed a nonlinear ODE which is now known as the (unforced) Van der Pol equation:
\[ \frac{d^2x}{dt^2} + a (x^2 - 1) \frac{dx}{dt} + x = 0 \]
where \( a \) is an adjustable parameter. Note that for \( a = 0 \) the above reduces to the equation for simple harmonic motion.
Write a MATLAB script, called vdp.m that solves the above
equation using the ode45 integrator. You will need to rewrite
the equation in canonical first order form, and code a MATLAB function,
say fcn_vdp.m, that evaluates the derivatives (right hand
sides) of the resulting
system of ODEs. Refer to the notes and online slides for the simple
harmonic oscillator if any of this seems unfamiliar.
Your script should solve the equation on the domain \[ 0 \le t \le 50 \]
using absolute and relative tolerances both set to 1.0e-6. You can hardcode
the parameter \( a \) to \( a = 5.0 \) in fcn_vdp.m.
Choose initial data
\[
x = 1.0
\]
\[
\frac{dx}{dt} = 4.0
\]
Your script should also make plots of the displacement, \( x(t) \), and
the phase space evolution, \( dx(t)/dt \) vs \( x(t) \).
If you have any concerns about the validity of your solution, consult the instructor before moving on.
Use vdp.m as the basis for another script vdp3.m
that integrates the Van der Pol equation with \( a = 5.0 \) and for the
3 sets of initial conditions \( [x,dx/dt] \) given by
\[
[2, -6]
\]
\[
[-1, 3]
\]
\[
[-3, -8]
\]
Make a single phase space plot that shows the evolution of the 3 sets of
initial data.
What can you say about the behaviour of the oscillator with respect to the initial conditions?
Create one more script file, vdp_ir.m that integrates
the Van der Pol equation for \( a = 2 \) (you may have to modify
vdp_fcn.m) and with 1025 output times uniformly distributed
from 0 to 50. Then compute and plot three levels of independent
residuals, both scaled and unscaled. Use \( O(\Delta t^2) \) finite
difference approximations to compute the independent residuals.
Do the residuals behave as you expect?
Again, refer to the notes and online slides if necessary, or ask for help.