More Cubic Fun

As sort of a quick follow-up to the last post, I decided that it was high time to start exploring (of course) the configuration space of cubics. More explorations on this forthcoming, but I'm testing out a new graphing calculator interface called Desmos. Let's see how well it can be shared.

(click the lower right to go to an interactive graph. What may be fun is playing around with the sliders. If I can figure out how to save it to a movie, I'll update it here)

Basically what's going on here is I'm trying to visualize the phase space of depressed cubics $x^3 + px + q$ (In our half-tablespoon example, $p = -3$ and $q = 1$) by seeing what happens when you map the unit disk in the complex plane. (If $p=0$ it's basically the cubing function which triples angles and makes things wind around the origin 3 times as often). It stands to reason that if $p$ and $q$ pass through certain interesting values, a qualitative change may occur (to relate this to some math that may be familiar, think the duscrumminanmanaent, I mean, discriminant (that word is just so annoyingly hard to say), of the quadratic: $b^2 - 4ac$ in $ax^2 + bx+c$: when it passes through 0, the nature of the solution changes). One day I may make a post detailing the connection between this kind visualization and topological understanding of the concept of solvability by radicals. It might be far in the future, though!

Here's a 3D rendering of the image of the unit disk as $p$ varies (with a vertical dimension rendered to distinguish distinct domain points that map to the same point (because the mapping is not 1-1. The mapping is not perfect, as the surface still self-intersects in places, but the direction the surface is facing as one traverses an intersection is distinct enough so you stay on the right path). (This is generated using SageMath, not Desmos, which doesn't yet have 3D stuff, oh well)

$r \leq 1, -3 \leq p \leq 3$

$r \leq 2, -6 \leq p \leq 6$, 

A Calculation in Statistical Mechanics

We've been harping about phase and state space for a long while, now, and have mentioned at more than one point that statistical mechanics is a subject dealing with, in some sense, the ultimate state space, of many gazillions of dimensions. It is amazing that we can even calculate anything about such spaces. So we take a break from tori, and return to the (music of the) spheres, today, to give an example. We consider $N$-dimensional spheres, for $N = 0, 1, 2, 3, \dots$ and their "volumes". Note, of course, by $N$-sphere, we always mean the curved surface of points in $(N+1)$-dimensional space with unit distance from the orgin, and not the whole ball. [The $3$-sphere, as we know, is foliated by nested tori. We aren't going TOO far afield].

To talk about this right, we need a little physics. Most of you probably have heard of the principle of conservation of energy... If you haven't, or even if you have, make sure to read about how Hamiltonian Mechanics is Awesome. Anyway, so if you imagine little billiard balls of uniform mass $m$ freely zipping about, banging into stuff, each particle having position $\mathbf{q}_i$ and momentum $\mathbf{p}_i$, then since we know Hamiltonian Mechanics is Awesome, this says that the two quantities evolve in such a way that the energy of the $i$th particle \[
E_i = \frac{|\mathbf{p}_i|^2}{2m}
\]
stays constant. But if you imagine having a very, very large number of such billiard balls (several billiards of such (1 billiard = $10^{15}$ in most dialects of English, or such equivalent cognates) HAR HAR!), and compute the total energy, summing them all up,\[
H = \sum_i  E_i = \sum_i \frac{|\mathbf{p}_i|^2}{2m}.
\] To say that the total energy is constant, then, is to say \[
\sum_i \frac{|\mathbf{p}_i|^2}{2m} = E
\] for some fixed number $E$. But what does this say? If we write out each \[
|\mathbf{p}_i|^2 = p_{ix}^2 + p_{iy}^2 + p_{iz}^2,
\] then \[
2mE=\sum_i p_{ix}^2 + p_{iy}^2 + p_{iz}^2,
\] meaning that if we sum over all the particles, say $N$ of them, where $N \approx 10^{24}$ (a milliard billiard of billiard balls), this means we're summing the squares of a $3N$-dimensional vector and equating it to a constant. So if we ask ourselves:

What is the collection of all possible momenta of $N$ particles that have one constant energy level $E$,

this is the same as saying,

What are all the points on the $(3N-1)$-dimensional sphere of radius $\sqrt{2mE}$?

(the -1 from the fact that you're expressing the surface as the solution set to one equation: the boundary of the $k$-ball is the $(k-1)$-sphere) Interesting, you might say, but is there any occasion which actually requires you to truly do calculations treating it precisely as a geometrical $(3N-1)$-dimensional object?

Well, what's one of the first things you think about when talking geometry? Areas and volumes, of course. There are situations in which you need to integrate over such energy "surface", or within the $N$-ball contained within, for $N \sim 10^{24}$. Basically, to find quantities like entropy, one needs to "count some number of cells" inside the $3N$-ball . This yields the total "number" of energy states with level less than or equal to $E$. "Number" in quotations because classically, of course, the actual number of states is the number of points on a sphere, which is of course uncountably infinite. We instead mean, per unit $3N$-dimensional volume, which actually leads to the density of states. The distinction is a little blurry, because, with a little help of quantum mechanics, we can think of each state as a cell of volume $(2\pi \hbar)^{3N}$, and each distinct state of momentum, which must be quantized, ultimately arises from boundary conditions in solving Schrödinger's equation, or eigenstuff.

Suffice to say, however, let's talk as if calculating the volume of the $3N$-ball does indeed count the total number of momentum states with energy less than or equal to $E$, i.e., let's calculate \[
\frac{1}{(2\pi \hbar)^{3N}}\int_{|\mathbf{p}|^2 \leq 2mE} d\mathbf{q}d\mathbf{p} =\frac{V^N}{(2\pi \hbar)^{3N}}\int_{|\mathbf{p}|^2 \leq 2mE} d\mathbf{p} .
\]
where $V$ is the ordinary 3D volume of the spatial region where all our particles live. For the momentum integral, we need the formula (whose derivation is fodder for another post!) for the volume of an $K$-ball \[
V_K(R) = \frac{\pi^{K/2}}{\Gamma\left(1+\frac{K}{2}\right)}R^K
\]
where $\Gamma(1+ K/2)$ is $(K/2)!$, the factorial of $K/2$, if $K$ is even, and $\pi^{1/2} K!!/2^{K/2}$ where $K!!$ is the double factorial $1\cdot 3 \cdot 5 \cdots k$, for $k$ odd.

Plugging this in with $R = \sqrt{2mE}$ and $K=3N$ we have the total volume of a $3N$-ball in momentum space: \[
V_{3N}(\sqrt{2mE}) = \frac{\pi^{3N/2}}{\Gamma\left(1+\frac{3N}{2}\right)} (2mE)^{3N/2}
\]
To make headway here, we can use Stirling's Approximation for the Gamma function:\[
\Gamma\left(1+\frac{3N}{2}\right) \approx \sqrt{3N\pi} \left(\frac{3N}{2e}\right)^{3N/2}
\] This gives \[
V_{3N}(\sqrt{2mE}) \approx \frac{1}{\sqrt{3\pi N}}\left(\frac{4\pi meE}{3N} \right)^{3N/2}
\] or, the total number of states is \[
\#\text{states} = \frac{V^N V_{3N}(\sqrt{2mE})}{(2\pi \hbar)^{3N}} \approx \frac{1}{\sqrt{3\pi N}}\left(\frac{4\pi meEV^{2/3}}{3N (2\pi \hbar)^2} \right)^{3N/2}
\]
Considering the smallness of $\hbar$ ($10^{-34}$ in everyday units) and the exponent $3N/2$ being on the order of $10^{24}$, this could be a staggeringly huge number (but that's what you expect: for lots of particles, they really can occupy even more states. For example, two particles behaving independently is the square the total number of states of for one particle). Thus to measure entropy, we take the natural log of all of this (and add a factor of the Boltzmann constant, $k$ which is on the order of $10^{-23}$ in ordinary units), we get \[
k\ln(\#\text{states}) = \frac{3Nk}{2} \left( \ln \left(\frac{4\pi mEV^{2/3}}{3N(2\pi \hbar)^2}\right) +1\right) -\frac{1}{2}\ln(3\pi N)k.
\]
Since those last terms are small potatoes (not being multiplied by $N$), we usually just drop it to arrive at the entropy of our system of particles as
\[
S = \frac{3Nk}{2}  \left( \ln \left(\frac{4\pi mEV^{2/3}}{3N(2\pi \hbar)^2}\right)\right)+\frac{3Nk}{2}.
\] And for you chemists, $Nk = nR$, where $n$ is the number of moles and $R$ is the gas constant:
\[
S = \frac{3nR}{2}\left( \ln \left(\frac{4\pi mEV^{2/3}}{3N(2\pi \hbar)^2}\right)\right)+\frac{3nR}{2}.
\]
That's really nasty. How about we try to plug numbers into this and see what it is? We'll do that more concretely next time. But suffice to say, let's get back to what's interesting: we calculated the volume of some high-dimensional stuff!

Hamiltonian Mechanics is Awesome (Double Pendulum, Part 4)

In this post, we work toward discretizing the equations for (and simulating) the double pendulum, considered in Parts 1, 2, and 3 of this series. It turns out that one very important qualitative behavior that we want is for the total energy to be conserved (one simulation I ran of this, directly discretizing the equations for $\ddot{\boldsymbol\theta}$, eventually ran amok over time, getting faster and faster, clearly something that would not occur in a real system). Thus, one place to look is symplectic methods, which we've remarked on before. However, at least with the symplectic methods I've learned, they seem more conceptually geared toward solving differential equations involving angular position and momentum, rather than angular position and velocity. This is not strictly true, as I have seen symplectic methods at work with position and velocity. However, this might be an artifact of the relation between velocity and momentum often being very trivial (e.g., just multiplying by some mass). In any case, the position-and-momentum viewpoint is more than important enough to be worth diving into now. This is referred to as Hamiltonian mechanics. We look at $H = T+V$ instead of $L = T-V$ and consider our variables as generalized position and momentum instead of generalized position and velocity (more generally, Hamiltonian mechanics can be derived from Lagrangian mechanics via the Legendre transform).

In fact, most basic physics courses tend to focus on the sum of kinetic and potential energies, saying energy is "converted" between these two forms. They just never used the word "Hamiltonian." It's not clear, at first, what the difference between the two approaches is, especially with simple examples where it's pretty trivial to go from one to the other. In mathematical or veteran Nested Tori reader terms, the difference is that of being equations on the tangent bundle vs. the cotangent bundle (called phase space). Momentum takes a life of its own in more complicated systems (and in fact, in quantum mechanics, it's not clear what something can even be the velocity of, so something called "momentum" crops up as its own fundamental quantity without reference to the velocity of anything).

So, what is this mystical momentum for our system? We've already calculated something which we called "angular momentum" for the double pendulum: $\mathbf{p} = \partial L/\partial \dot{\boldsymbol \theta}$, the partial "gradient" vector with respect to the angular velocities. This is what it is generally. (As a check, we note for a single particle of mass $m$ and kinetic energy $\frac 1 2 m v^2$, $\partial L/\partial \dot{\mathbf{x}}=\partial L/\partial {\mathbf{v}}$ is indeed the usual $m\mathbf{v}$.) For our example, we derived
\[
\mathbf{p} =  \frac{\partial L}{\partial \dot{\boldsymbol{ \theta}}} = \begin{pmatrix} (m_1 +m_2) \ell_1^2 & m_2\ell_1\ell_2\cos(\theta_2-\theta_1) \\ m_2 \ell_1 \ell_2 \cos(\theta_2 -\theta_1) & m_2 \ell_2^2\end{pmatrix} \begin{pmatrix}\dot \theta_1 \\ \dot\theta_2 \end{pmatrix} = I\dot{\boldsymbol{\theta}}
\]
That's almost something like the case $\mathbf{p} = m\mathbf{v}$ for a single particle, with the rotational analogues of mass and velocity, though the "mass" here is much more complicated (depending on the positions). In particular, as a vector, it doesn't have to point in the same direction as the velocity, unlike in the particle case (that misalignment is also the cause of washing machines thumping with an unbalanced load). With this in mind, we rewrite the kinetic energy in terms of $\mathbf{p}$ as follows: assuming the matrix $I$ is invertible, we write $\dot{\boldsymbol\theta} = I^{-1} \mathbf{p}$, substitute in $T = \frac 1 2 \dot{\boldsymbol{\theta}}^t I \dot{\boldsymbol{\theta}}$, and write $\mathbf{q} = \boldsymbol{\theta}$ (the usual notation), to get
\[
H(\mathbf{q},\mathbf{p}) = T+V= \frac{1}{2} (\mathbf{p}^t I^{-1}) I (I^{-1} \mathbf{p}) + V(\mathbf q) =\frac{1}{2} \mathbf{p}^t I^{-1}  \mathbf{p} + V(\mathbf q).
\]

Hamilton's Equations

Now that we have derived $H$, the Hamiltonian, in order to relate it to our actual equations of motion, we have Hamilton's equations:
 \[
 \dot{\mathbf{q}} = \frac{\partial H}{\partial \mathbf{p}}
 \] \[
\dot{\mathbf{p}} = -\frac{\partial H}{\partial \mathbf{q}}
\]
which are very nice: $\mathbf{p}$ and $\mathbf{q}$ are almost "dual" to each other in some sense. There's just that extra minus sign; this can be remembered by recalling that $H$ includes potential energy and thus differentiating with respect to $\mathbf{q}$ gives the gradient. But force is conventionally given as the negative of the gradient of the potential in physics. We thus find, plugging in the above:
\[
\dot{\mathbf{q}} = I^{-1} \mathbf{p}
\]
and
\[
\dot{\mathbf{p}} = -\frac 1 2 \mathbf{p}^t \frac{\partial I^{-1}}{\partial \mathbf{q}} \mathbf{p} -\frac{\partial V}{\partial \mathbf{q}}.
\]
The part that really complicates things and makes this different from the simpler examples is that the moment of inertia $I$ is dependent on $\mathbf{q}$, so we can't ignore that term when differentiating with respect to $\mathbf{q}$ (it is often just the $-\partial V/\partial\mathbf{q}$ term by itself). How do we calculate that derivative of $I^{-1}$? We step back down to components:
\[
\dot{p}_i = -\frac 1 2 \mathbf{p}^t \frac{\partial I^{-1}}{\partial \theta_i} \mathbf{p} -\frac{\partial V}{\partial \theta_i}= +\frac 1 2 \mathbf{p}^t I^{-1}\frac{\partial I}{\partial \theta_i}I^{-1} \mathbf{p} -\frac{\partial V}{\partial \theta_i}
\] \[
= \frac{1}{2} \dot{\mathbf q}^t \frac{\partial I}{\partial \theta_i}\dot{\mathbf q}-\frac{\partial V}{\partial \theta_i}
\]
where we have employed the identity that the derivative of the inverse of a matrix is always multiplication, by the inverse on both sides, of the derivative of the matrix (and a sign reversal). Let's write out what all of this is. We note that $\partial I/\partial \theta_1 = - \partial I/\partial \theta_2$, because the only way the angles in appear in $I$ is the difference $\theta_2 - \theta_1$. Thus, we have
\[
\dot {\mathbf q} = I^{-1} \mathbf{p},
\]
and using that as an intermediate variable (because to write it all out explicitly is a real mess, and we will need to use it when discretizing it, anyway):
\[
\dot{p}_1 = \frac{1}{2}\dot{\mathbf{q}}^t \begin{pmatrix} 0 & m_2\ell_1\ell_2 \sin(\theta_2-\theta_1)  \\ m_2\ell_1\ell_2 \sin(\theta_2-\theta_1) &0\end{pmatrix} \dot{\mathbf{q}} + (m_1 + m_2) g\ell_1 \sin \theta_1
\] \[
\dot{p}_2 = \frac{1}{2}\dot{\mathbf{q}}^t \begin{pmatrix} 0 & -m_2\ell_1\ell_2 \sin(\theta_2-\theta_1)  \\ -m_2\ell_1\ell_2 \sin(\theta_2-\theta_1) &0\end{pmatrix} \dot{\mathbf{q}} +  m_2 g\ell_2 \sin \theta_2.
\]
In the final part, we discretize these equations. That'll have to wait for next time.

The Euler-Lagrange equations for the Double Pendulum (Config Spaces, Part 3)

In this post, continuing the explorations of the double pendulum (see Part 1 and Part 2) we concentrate on deriving its equation of motion (the Euler-Lagrange equation). These differential equations are the heart of Lagrangian mechanics, and indeed really what one tries to get to when applying the methods (it's essentially a way of getting Newton's 2nd Law for complicated systems). And why it is awesome. This is, like the previous part, going to be a more technical than visual, but we'll get back to that in Part 4.

Recall from Part 2 that for the double pendulum,
\[
L(\theta,\dot \theta) = T - V = \frac 1 2 (m_1+m_2) \ell_1^2 \dot{\theta}_1^2 +\frac 1 2 m_2\left( 2 \dot{\theta}_1 \dot{\theta}_2\ell_1 \ell_2 \cos(\theta_2 - \theta_1) + \ell_2^2 \dot{\theta}_2^2\right)\] \[+ (m_1+m_2) g \ell_1 \cos \theta_1 + m_2 g \ell_2 \cos \theta_2,
\]
where, $\theta_1$, $\theta_2$ are the angles from vertical of the two pendulums, $\dot \theta_1$, $\dot \theta_2$ are their angular speeds, $\ell_1$, $\ell_2$ are the length of the rods, and $m_1$, $m_2$ are their masses. This looks very imposing, so let's remark a little on the structure of that equation. In simpler situations, one often sees different squared speed terms, but here we have the interaction of the two parts of the system, in terms with products like $\dot \theta_1 \dot \theta_2$, as well as the $\cos(\theta_2 - \theta_1)$. This distinguishes the system from two completely isolated (uncoupled) pendulums; the Lagrangian already captures the aspect that energy gets transferred back and forth between the two pendulums.

The equations of motion are derived from considering an optimization problem, namely, what paths $\gamma$ through state space connecting two states optimizes the action integral
\[
\int_a^b L(\gamma(t),\dot\gamma(t)) dt.
\]
In a process quite similar to what's done in first-year calculus, we "take a derivative and set it to zero", which derives the equations
\[
\frac{\partial L}{\partial \theta_i} - \frac d {dt} \frac{\partial L}{\partial \dot \theta_i} = 0
\]
where $\partial L/\partial \dot \theta_i$ means the partial derivative with respect to $\dot \theta_i$ as if it were an independent variable. In fact, one way is to start off with a formula involving two independent variables $\theta$ and $\omega$; only after solving for the extremal path do we actually find $\omega = \dot \theta$ (this is a generalization called the Hamilton-Pontryagin principle).

$\partial L/\partial \dot {\boldsymbol \theta}$, angular momentum, and the moment of inertia

Let's now do the computation $\frac{\partial L}{\partial \dot \theta_1}$ (this quantity actually is the angular momentum of the first pendulum, $p_1$)—treating everything except $\dot \theta_1$ as a constant:
\[
p_1 = \frac{\partial L}{\partial \dot \theta_1} = (m_1 + m_2) \ell_1^2 \dot \theta_1 + m_2  \ell_1\ell_2
 \cos(\theta_2-\theta_1)\dot \theta_2 .
\]
Note that unlike the more basic examples of angular momentum, this isn't simply dependent only on the (angular) velocity $\dot \theta_1$, but also on $\dot \theta_2$ and nonlinearly on $\theta_1$ and $\theta_2$. Similarly,
\[
p_2 = \frac{\partial L}{\partial \dot \theta_2} = m_2 \ell_2^2 \dot \theta_2 + m_2\ell_1\ell_2 \cos(\theta_2 - \theta_1) \dot \theta_1.
\]
This is more compactly expressed as
\[
\frac{\partial L}{\partial \dot{\boldsymbol{\theta}}}=\begin{pmatrix}p_1\\p_2\end{pmatrix} = \begin{pmatrix} (m_1 +m_2) \ell_1^2 & m_2\ell_1\ell_2\cos(\theta_2-\theta_1) \\ m_2 \ell_1 \ell_2 \cos(\theta_2 -\theta_1) & m_2 \ell_2^2\end{pmatrix} \begin{pmatrix}\dot \theta_1 \\ \dot\theta_2 \end{pmatrix}.
\]
This helps keep the unruliness under one umbrella. The square matrix in the above will be important, so let's call it $I$ (it is something like a moment of inertia for our system—but it's more complicated, principally because there is an interaction between its two component parts... It's not even a rigid body!). It should also be noted that the kinetic energy is $T = \frac 1 2 \dot{\boldsymbol{\theta}}^t I \dot{\boldsymbol \theta}$, so it still is quadratic in $\dot {\boldsymbol \theta}$, just like good ol' $\frac{1}{2} mv^2$.

The forces, $\partial L/\partial \boldsymbol{\theta}$

Due to the nonlinear dependence of the Lagrangian on the angular variables, the derivatives are trickier. But we forge on, because sticking to just the simplest examples won't help us gain as much insight into how things are used in practice. We have
\[
\frac{\partial L}{\partial \theta_1} = m_2 \ell_1 \ell_2 \dot \theta_1 \dot \theta_2 \sin(\theta_2-\theta_1) - (m_1 + m_2) g \ell_1 \sin \theta_1
\]

and
\[
\frac{\partial L}{\partial \theta_2} = -m_2 \ell_1 \ell_2 \dot \theta_1 \dot \theta_2 \sin(\theta_2-\theta_1) - m_2 g \ell_2 \sin \theta_2.
\]

Putting it together

Now we take the time derivatives of the momenta we derived above and subtract and set to zero to derive the Euler-Lagrange equations. That's straightfoward calculus; but different from simpler examples in that, of course, the chain rule must be applied to the $\cos(\theta_2 - \theta_1)$, so there will be an extra $\dot \theta_1$ or $\dot \theta_2$ term:

\[
\frac{d}{dt} \frac{\partial L}{\partial \dot{\boldsymbol{\theta}}} =  \begin{pmatrix} 0 & -m_2\ell_1\ell_2\sin(\theta_2-\theta_1)(\dot \theta_2 - \dot \theta_1) \\ -m_2 \ell_1 \ell_2 \sin(\theta_2 -\theta_1)(\dot \theta_2 - \dot \theta_1) & 0\end{pmatrix} \begin{pmatrix}\dot \theta_1 \\ \dot\theta_2 \end{pmatrix}\]
\[+\begin{pmatrix} (m_1 +m_2) \ell_1^2 & m_2\ell_1\ell_2\cos(\theta_2-\theta_1) \\ m_2 \ell_1 \ell_2 \cos(\theta_2 -\theta_1) & m_2 \ell_2^2\end{pmatrix} \begin{pmatrix}\ddot \theta_1 \\ \ddot\theta_2 \end{pmatrix}.
\]
The nice thing about it is that when it is subtracted from $\partial L/\partial \boldsymbol \theta$, the cross terms with the $\dot \theta_1 \dot\theta_2$ cancel. Also, the matrix multiplying the $\ddot{\theta}_i$'s is none other than our moment-of-inertia matrix $I$. This finally gives us

\[
0=\frac{\partial L}{\partial \boldsymbol\theta} - \frac{d}{dt} \frac{\partial L}{\partial \dot{\boldsymbol\theta}} =\begin{pmatrix}-(m_1+m_2) g \ell_1 \sin \theta_1 - m_2 \ell_1\ell_2 \sin(\theta_2-\theta_1)\dot{\theta}_2^2 \\ -m_2 g \ell_2 \sin \theta_2 +m_2 \ell_1\ell_2 \sin(\theta_2-\theta_1)\dot{\theta}_1^2  \end{pmatrix} -I \begin{pmatrix} \ddot \theta_1 \\ \ddot \theta_2 \end{pmatrix}
\]
as the Euler-Lagrange equations. As the final equation of motion, we bring the $\ddot {\boldsymbol \theta}$ term to the other side:
\[
I \begin{pmatrix} \ddot \theta_1 \\ \ddot \theta_2 \end{pmatrix}  = \begin{pmatrix}-(m_1+m_2) g \ell_1 \sin \theta_1 - m_2 \ell_1\ell_2 \sin(\theta_2-\theta_1)\dot{\theta}_2^2 \\ -m_2 g \ell_2 \sin \theta_2 +m_2 \ell_1\ell_2 \sin(\theta_2-\theta_1)\dot{\theta}_1^2  \end{pmatrix}
\]
(keep in mind that $I$ varies with $\theta_1$ and $\theta_2$, however). We discretize and solve this next time. For now, simply take note of the interesting nonlinear interactions: in the positions as the sine of the difference, and quadratic nonlinearity in the derivatives (note the reversal: the $\ddot\theta_1$ gets a $\dot \theta_2^2$ and vice versa).