# start and final times
t0 = 0.0
tn = 1.0
# Period
t = np.linspace(t0, tn, n, endpoint=False)
SqrSignal = np.zeros(n)
SqrSignal = 1.0+signal.square(2*np.pi*5*t)
plt.plot(t, SqrSignal)
plt.ylim(-0.5, 2.5)
plt.show()
For the sinusoidal example studied in the previous subsections the period is \( \tau=2\pi/\omega \). However, higher harmonics can also satisfy the periodicity requirement. In general, any force that satisfies the periodicity requirement can be expressed as a sum over harmonics, $$ \begin{equation} F(t)=\frac{f_0}{2}+\sum_{n>0} f_n\cos(2n\pi t/\tau)+g_n\sin(2n\pi t/\tau). \label{_auto62} \end{equation} $$
From the previous subsection, one can write down the answer for \( x_{pn}(t) \), by substituting \( f_n/m \) or \( g_n/m \) for \( F_0/m \) into Eq.s \eqref{eq:fastdriven1} or \eqref{eq:fastdriven2} respectively. By writing each factor \( 2n\pi t/\tau \) as \( n\omega t \), with \( \omega\equiv 2\pi/\tau \), $$ \begin{equation} \label{eq:fourierdef1} F(t)=\frac{f_0}{2}+\sum_{n>0}f_n\cos(n\omega t)+g_n\sin(n\omega t). \end{equation} $$
The solutions for \( x(t) \) then come from replacing \( \omega \) with \( n\omega \) for each term in the particular solution in Equations \eqref{eq:partform} and \eqref{eq:Ddrive}, $$ \begin{eqnarray} x_p(t)&=&\frac{f_0}{2k}+\sum_{n>0} \alpha_n\cos(n\omega t-\delta_n)+\beta_n\sin(n\omega t-\delta_n),\\ \nonumber \alpha_n&=&\frac{f_n/m}{\sqrt{((n\omega)^2-\omega_0^2)+4\beta^2n^2\omega^2}},\\ \nonumber \beta_n&=&\frac{g_n/m}{\sqrt{((n\omega)^2-\omega_0^2)+4\beta^2n^2\omega^2}},\\ \nonumber \delta_n&=&\tan^{-1}\left(\frac{2\beta n\omega}{\omega_0^2-n^2\omega^2}\right). \end{eqnarray} $$
Because the forces have been applied for a long time, any non-zero damping eliminates the homogenous parts of the solution, so one need only consider the particular solution for each \( n \).
The problem will considered solved if one can find expressions for the coefficients \( f_n \) and \( g_n \), even though the solutions are expressed as an infinite sum. The coefficients can be extracted from the function \( F(t) \) by $$ \begin{eqnarray} \label{eq:fourierdef2} f_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~F(t)\cos(2n\pi t/\tau),\\ \nonumber g_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~F(t)\sin(2n\pi t/\tau). \end{eqnarray} $$
To check the consistency of these expressions and to verify Eq. \eqref{eq:fourierdef2}, one can insert the expansion of \( F(t) \) in Eq. \eqref{eq:fourierdef1} into the expression for the coefficients in Eq. \eqref{eq:fourierdef2} and see whether $$ \begin{eqnarray} f_n&=?&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~\left\{ \frac{f_0}{2}+\sum_{m>0}f_m\cos(m\omega t)+g_m\sin(m\omega t) \right\}\cos(n\omega t). \end{eqnarray} $$
Immediately, one can throw away all the terms with \( g_m \) because they convolute an even and an odd function. The term with \( f_0/2 \) disappears because \( \cos(n\omega t) \) is equally positive and negative over the interval and will integrate to zero. For all the terms \( f_m\cos(m\omega t) \) appearing in the sum, one can use angle addition formulas to see that \( \cos(m\omega t)\cos(n\omega t)=(1/2)(\cos[(m+n)\omega t]+\cos[(m-n)\omega t] \). This will integrate to zero unless \( m=n \). In that case the \( m=n \) term gives $$ \begin{equation} \int_{-\tau/2}^{\tau/2}dt~\cos^2(m\omega t)=\frac{\tau}{2}, \label{_auto63} \end{equation} $$
and $$ \begin{eqnarray} f_n&=?&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2} dt~f_n/2\\ \nonumber &=&f_n~\checkmark. \end{eqnarray} $$
The same method can be used to check for the consistency of \( g_n \).
Consider the driving force: $$ \begin{equation} F(t)=At/\tau,~~-\tau/2 < t < \tau/2,~~~F(t+\tau)=F(t). \label{_auto64} \end{equation} $$
Find the Fourier coefficients \( f_n \) and \( g_n \) for all \( n \) using Eq. \eqref{eq:fourierdef2}.
Only the odd coefficients enter by symmetry, i.e. \( f_n=0 \). One can find \( g_n \) integrating by parts, $$ \begin{eqnarray} \label{eq:fouriersolution} g_n&=&\frac{2}{\tau}\int_{-\tau/2}^{\tau/2}dt~\sin(n\omega t) \frac{At}{\tau}\\ \nonumber u&=&t,~dv=\sin(n\omega t)dt,~v=-\cos(n\omega t)/(n\omega),\\ \nonumber g_n&=&\frac{-2A}{n\omega \tau^2}\int_{-\tau/2}^{\tau/2}dt~\cos(n\omega t) +\left.2A\frac{-t\cos(n\omega t)}{n\omega\tau^2}\right|_{-\tau/2}^{\tau/2}. \end{eqnarray} $$
The first term is zero because \( \cos(n\omega t) \) will be equally positive and negative over the interval. Using the fact that \( \omega\tau=2\pi \), $$ \begin{eqnarray} g_n&=&-\frac{2A}{2n\pi}\cos(n\omega\tau/2)\\ \nonumber &=&-\frac{A}{n\pi}\cos(n\pi)\\ \nonumber &=&\frac{A}{n\pi}(-1)^{n+1}. \end{eqnarray} $$
More text will come here, chpater 5.7-5.8 of Taylor are discussed during the lectures. The code here uses the Fourier series discussed in chapter 5.7 for a square wave signal. The equations for the coefficients are are discussed in Taylor section 5.7, see Example 5.4. The code here visualizes the various approximations given by Fourier series compared with a square wave with period \( T=0.2 \), witth \( 0.1 \) and max value \( F=2 \). We see that when we increase the number of components in the Fourier series, the Fourier series approximation gets closes and closes to the square wave signal.
import numpy as np
import math
from scipy import signal
import matplotlib.pyplot as plt
# number of points
n = 500
# start and final times
t0 = 0.0
tn = 1.0
# Period
T =0.2
# Max value of square signal
Fmax= 2.0
# Width of signal
Width = 0.1
t = np.linspace(t0, tn, n, endpoint=False)
SqrSignal = np.zeros(n)
FourierSeriesSignal = np.zeros(n)
SqrSignal = 1.0+signal.square(2*np.pi*5*t+np.pi*Width/T)
a0 = Fmax*Width/T
FourierSeriesSignal = a0
Factor = 2.0*Fmax/np.pi
for i in range(1,500):
FourierSeriesSignal += Factor/(i)*np.sin(np.pi*i*Width/T)*np.cos(i*t*2*np.pi/T)
plt.plot(t, SqrSignal)
plt.plot(t, FourierSeriesSignal)
plt.ylim(-0.5, 2.5)
plt.show()
The material here was discussed during the lecture of February 19 and 21. It is also covered by Taylor in section 5.8.
Consider a particle at rest in the bottom of an underdamped harmonic oscillator, that then feels a sudden impulse, or change in momentum, \( I=F\Delta t \) at \( t=0 \). This increases the velocity immediately by an amount \( v_0=I/m \) while not changing the position. One can then solve the trajectory by solving Eq. \eqref{eq:homogsolution} with initial conditions \( v_0=I/m \) and \( x_0=0 \). This gives $$ \begin{equation} x(t)=\frac{I}{m\omega'}e^{-\beta t}\sin\omega't, ~~t>0. \label{_auto65} \end{equation} $$
Here, \( \omega'=\sqrt{\omega_0^2-\beta^2} \). For an impulse \( I_i \) that occurs at time \( t_i \) the trajectory would be $$ \begin{equation} x(t)=\frac{I_i}{m\omega'}e^{-\beta (t-t_i)}\sin[\omega'(t-t_i)] \Theta(t-t_i), \label{_auto66} \end{equation} $$
where \( \Theta(t-t_i) \) is a step function, i.e. \( \Theta(x) \) is zero for \( x < 0 \) and unity for \( x>0 \). If there were several impulses linear superposition tells us that we can sum over each contribution, $$ \begin{equation} x(t)=\sum_i\frac{I_i}{m\omega'}e^{-\beta(t-t_i)}\sin[\omega'(t-t_i)]\Theta(t-t_i) \label{_auto67} \end{equation} $$
Now one can consider a series of impulses at times separated by \( \Delta t \), where each impulse is given by \( F_i\Delta t \). The sum above now becomes an integral, $$ \begin{eqnarray}label{eq:Greeny} x(t)&=&\int_{-\infty}^\infty dt'~F(t')\frac{e^{-\beta(t-t')}\sin[\omega'(t-t')]}{m\omega'}\Theta(t-t')\\ \nonumber &=&\int_{-\infty}^\infty dt'~F(t')G(t-t'),\\ \nonumber G(\Delta t)&=&\frac{e^{-\beta\Delta t}\sin[\omega' \Delta t]}{m\omega'}\Theta(\Delta t) \end{eqnarray} $$
The quantity \( e^{-\beta(t-t')}\sin[\omega'(t-t')]/m\omega'\Theta(t-t') \) is called a Green's function, \( G(t-t') \). It describes the response at \( t \) due to a force applied at a time \( t' \), and is a function of \( t-t' \). The step function ensures that the response does not occur before the force is applied. One should remember that the form for \( G \) would change if the oscillator were either critically- or over-damped.
When performing the integral in Eq. \eqref{eq:Greeny} one can use angle addition formulas to factor out the part with the \( t' \) dependence in the integrand, $$ \begin{eqnarray} \label{eq:Greeny2} x(t)&=&\frac{1}{m\omega'}e^{-\beta t}\left[I_c(t)\sin(\omega't)-I_s(t)\cos(\omega't)\right],\\ \nonumber I_c(t)&\equiv&\int_{-\infty}^t dt'~F(t')e^{\beta t'}\cos(\omega't'),\\ \nonumber I_s(t)&\equiv&\int_{-\infty}^t dt'~F(t')e^{\beta t'}\sin(\omega't'). \end{eqnarray} $$
If the time \( t \) is beyond any time at which the force acts, \( F(t'>t)=0 \), the coefficients \( I_c \) and \( I_s \) become independent of \( t \).
Consider an undamped oscillator (\( \beta\rightarrow 0 \)), with characteristic frequency \( \omega_0 \) and mass \( m \), that is at rest until it feels a force described by a Gaussian form, $$ \begin{eqnarray*} F(t)&=&F_0 \exp\left\{\frac{-t^2}{2\tau^2}\right\}. \end{eqnarray*} $$
For large times (\( t>>\tau \)), where the force has died off, find \( x(t) \).\\ Solve for the coefficients \( I_c \) and \( I_s \) in Eq. \eqref{eq:Greeny2}. Because the Gaussian is an even function, \( I_s=0 \), and one need only solve for \( I_c \), $$ \begin{eqnarray*} I_c&=&F_0\int_{-\infty}^\infty dt'~e^{-t^{\prime 2}/(2\tau^2)}\cos(\omega_0 t')\\ &=&\Re F_0 \int_{-\infty}^\infty dt'~e^{-t^{\prime 2}/(2\tau^2)}e^{i\omega_0 t'}\\ &=&\Re F_0 \int_{-\infty}^\infty dt'~e^{-(t'-i\omega_0\tau^2)^2/(2\tau^2)}e^{-\omega_0^2\tau^2/2}\\ &=&F_0\tau \sqrt{2\pi} e^{-\omega_0^2\tau^2/2}. \end{eqnarray*} $$
The third step involved completing the square, and the final step used the fact that the integral $$ \begin{eqnarray*} \int_{-\infty}^\infty dx~e^{-x^2/2}&=&\sqrt{2\pi}. \end{eqnarray*} $$
To see that this integral is true, consider the square of the integral, which you can change to polar coordinates, $$ \begin{eqnarray*} I&=&\int_{-\infty}^\infty dx~e^{-x^2/2}\\ I^2&=&\int_{-\infty}^\infty dxdy~e^{-(x^2+y^2)/2}\\ &=&2\pi\int_0^\infty rdr~e^{-r^2/2}\\ &=&2\pi. \end{eqnarray*} $$
Finally, the expression for \( x \) from Eq. \eqref{eq:Greeny2} is $$ \begin{eqnarray*} x(t>>\tau)&=&\frac{F_0\tau}{m\omega_0} \sqrt{2\pi} e^{-\omega_0^2\tau^2/2}\sin(\omega_0t). \end{eqnarray*} $$
Let us end our discussion of oscillations with another classical case, the pendulum.
The angular equation of motion of the pendulum is given by Newton's equation and with no external force it reads $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \label{_auto68} \end{equation} $$ with an angular velocity and acceleration given by $$ \begin{equation} v=l\frac{d\theta}{dt}, \label{_auto69} \end{equation} $$ and $$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \label{_auto70} \end{equation} $$
We do however expect that the motion will gradually come to an end due a viscous drag torque acting on the pendulum. In the presence of the drag, the above equation becomes $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \label{_auto71} \end{equation} $$ where \( \nu \) is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then $$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \label{_auto72} \end{equation} $$ with \( A \) and \( \omega \) two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.
We define $$ \omega_0=\sqrt{g/l}, $$ the so-called natural frequency and the new dimensionless quantities $$ \hat{t}=\omega_0t, $$ with the dimensionless driving frequency $$ \hat{\omega}=\frac{\omega}{\omega_0}, $$ and introducing the quantity \( Q \), called the quality factor, $$ Q=\frac{mg}{\omega_0\nu}, $$ and the dimensionless amplitude $$ \hat{A}=\frac{A}{mg} $$
We have $$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$ This equation can in turn be recast in terms of two coupled first-order differential equations as follows $$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$ and $$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$
These are the equations to be solved. The factor \( Q \) represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude \( \hat{A} \) is measured in units of the maximum possible gravitational torque while \( \hat{\omega} \) is the angular frequency of the external torque measured in units of the pendulum's natural frequency.
The gravitational potential energy and forces involving two masses \( a \) and \( b \) are $$ \begin{eqnarray} U_{ab}&=&-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|},\\ \nonumber F_{ba}&=&-\frac{Gm_am_b}{|\boldsymbol{r}_a-\boldsymbol{r}_b|^2}\hat{r}_{ab},\\ \nonumber \hat{r}_{ab}&=&\frac{\boldsymbol{r}_b-\boldsymbol{r}_a}{|\boldsymbol{r}_a-\boldsymbol{r}_b|}. \end{eqnarray} $$
Here \( G=6.67\times 10^{-11} \) Nm$^2$/kg$^2$, and \( F_{ba} \) is the force on \( b \) due to \( a \). By inspection, one can see that the force on \( b \) due to \( a \) and the force on \( a \) due to \( b \) are equal and opposite. The net potential energy for a large number of masses would be $$ \begin{equation} U=\sum_{a < b}U_{ab}=\frac{1}{2}\sum_{a\ne b}U_{ab}. \label{_auto73} \end{equation} $$
Thus far, we have considered the trajectory as if the force is centered around a fixed point. For two bodies interacting only with one another, both masses circulate around the center of mass. One might think that solutions would become more complex when both particles move, but we will see here that the problem can be reduced to one with a single body moving according to a fixed force by expressing the trajectories for \( \boldsymbol{r}_1 \) and \( \boldsymbol{r}_2 \) into the center-of-mass coordinate \( \boldsymbol{R}_{\rm cm} \) and the relative coordinate \( \boldsymbol{r} \), $$ \begin{eqnarray} \boldsymbol{R}_{\rm cm}&\equiv&\frac{m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2}{m_1+m_2},\\ \nonumber \boldsymbol{r}&\equiv&\boldsymbol{r}_1-\boldsymbol{r_2}. \end{eqnarray} $$
Here, we assume the two particles interact only with one another, so \( \boldsymbol{F}_{12}=-\boldsymbol{F}_{21} \) (where \( \boldsymbol{F}_{ij} \) is the force on \( i \) due to \( j \). The equations of motion then become $$ \begin{eqnarray} \ddot{\boldsymbol{R}}_{\rm cm}&=&\frac{1}{m_1+m_2}\left\{m_1\ddot{\boldsymbol{r}}_1+m_2\ddot{\boldsymbol{r}}_2\right\}\\ \nonumber &=&\frac{1}{m_1+m_2}\left\{\boldsymbol{F}_{12}+\boldsymbol{F}_{21}\right\}=0.\\ \ddot{\boldsymbol{r}}&=&\ddot{\boldsymbol{r}}_1-\ddot{\boldsymbol{r}}_2=\left(\frac{\boldsymbol{F}_{12}}{m_1}-\frac{\boldsymbol{F}_{21}}{m_2}\right)\\ \nonumber &=&\left(\frac{1}{m_1}+\frac{1}{m_2}\right)\boldsymbol{F}_{12}. \end{eqnarray} $$
The first expression simply states that the center of mass coordinate \( \boldsymbol{R}_{\rm cm} \) moves at a fixed velocity. The second expression can be rewritten in terms of the reduced mass \( \mu \). $$ \begin{eqnarray} \mu \ddot{\boldsymbol{r}}&=&\boldsymbol{F}_{12},\\ \frac{1}{\mu}&=&\frac{1}{m_1}+\frac{1}{m_2},~~~~\mu=\frac{m_1m_2}{m_1+m_2}. \end{eqnarray} $$
Thus, one can treat the trajectory as a one-body problem where the reduced mass is \( \mu \), and a second trivial problem for the center of mass. The reduced mass is especially convenient when one is considering gravitational problems because then $$ \begin{eqnarray} \mu \ddot{r}&=&-\frac{Gm_1m_2}{r^2}\hat{r}\\ \nonumber &=&-\frac{GM\mu}{r^2}\hat{r},~~~M\equiv m_1+m_2. \end{eqnarray} $$
For the gravitational problem, the reduced mass then falls out and the trajectory depends only on the total mass \( M \).
The kinetic energy and momenta also have analogues in center-of-mass coordinates. The total and relative momenta are $$ \begin{eqnarray} \boldsymbol{P}&\equiv&\boldsymbol{p}_1+\boldsymbol{p}_2=M\dot{\boldsymbol{R}}_{\rm cm},\\ \nonumber \boldsymbol{q}&\equiv&\mu\dot{\boldsymbol{r}}. \end{eqnarray} $$
With these definitions, a little algebra shows that the kinetic energy becomes $$ \begin{eqnarray} T&=&\frac{1}{2}m_1|\boldsymbol{v}_1|^2+\frac{1}{2}m_2|\boldsymbol{v}_2|^2\\ \nonumber &=&\frac{1}{2}M|\dot{\boldsymbol{R}}_{\rm cm}|^2 +\frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2\\ \nonumber &=&\frac{P^2}{2M}+\frac{q^2}{2\mu}. \end{eqnarray} $$
The standard strategy is to transform into the center of mass frame, then treat the problem as one of a single particle of mass \( \mu \) undergoing a force \( \boldsymbol{F}_{12} \). Scattering angles can also be expressed in this frame, then transformed into the lab frame. In practice, one sees examples in the literature where \( d\sigma/d\Omega \) expressed in both the "center-of-mass" and in the "laboratory" frame.
Kepler's laws state that a gravitational orbit should be an ellipse with the source of the gravitational field at one focus. Deriving this is surprisingly messy. To do this, we first use angular momentum conservation to transform the equations of motion so that it is in terms of \( r \) and \( \theta \) instead of \( r \) and \( t \). The overall strategy is to
Recognizing that the numerator of the third term is the velocity squared, and that it can be written in polar coordinates, $$ \begin{equation} v^2=\dot{x}^2+\dot{y}^2=\dot{r}^2+r^2\dot{\theta}^2, \label{_auto74} \end{equation} $$
one can write \( \ddot{r} \) as $$ \begin{eqnarray} \label{eq:radialeqofmotion2} \ddot{r}&=&\frac{F_x\cos\theta+F_y\sin\theta}{m}+\frac{\dot{r}^2+r^2\dot{\theta}^2}{r}-\frac{\dot{r}^2}{r}\\ \nonumber &=&\frac{F}{m}+\frac{r^2\dot{\theta}^2}{r}\\ \nonumber m\ddot{r}&=&F+\frac{L^2}{mr^3}. \end{eqnarray} $$
This derivation used the fact that the force was radial, \( F=F_r=F_x\cos\theta+F_y\sin\theta \), and that angular momentum is \( L=mrv_{\theta}=mr^2\dot{\theta} \). The term \( L^2/mr^3=mv^2/r \) behaves like an additional force. Sometimes this is referred to as a centrifugal force, but it is not a force. Instead, it is the consequence of considering the motion in a rotating (and therefore accelerating) frame.
Now, we switch to the particular case of an attractive inverse square force, \( F=-\alpha/r^2 \), and show that the trajectory, \( r(\theta) \), is an ellipse. To do this we transform derivatives w.r.t. time to derivatives w.r.t. \( \theta \) using the chain rule combined with angular momentum conservation, \( \dot{\theta}=L/mr^2 \). $$ \begin{eqnarray} \label{eq:rtotheta} \dot{r}&=&\frac{dr}{d\theta}\dot{\theta}=\frac{dr}{d\theta}\frac{L}{mr^2},\\ \nonumber \ddot{r}&=&\frac{d^2r}{d\theta^2}\dot{\theta}^2 +\frac{dr}{d\theta}\left(\frac{d}{dr}\frac{L}{mr^2}\right)\dot{r}\\ \nonumber &=&\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2 -2\frac{dr}{d\theta}\frac{L}{mr^3}\dot{r}\\ \nonumber &=&\frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2 -\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2 \end{eqnarray} $$
Equating the two expressions for \( \ddot{r} \) in Eq.s \eqref{eq:radialeqofmotion2} and \eqref{eq:rtotheta} eliminates all the derivatives w.r.t. time, and provides a differential equation with only derivatives w.r.t. \( \theta \), $$ \begin{equation} \label{eq:rdotdot} \frac{d^2r}{d\theta^2}\left(\frac{L}{mr^2}\right)^2 -\frac{2}{r}\left(\frac{dr}{d\theta}\right)^2\left(\frac{L}{mr^2}\right)^2 =\frac{F}{m}+\frac{L^2}{m^2r^3}, \end{equation} $$
that when solved yields the trajectory, i.e. \( r(\theta) \). Up to this point the expressions work for any radial force, not just forces that fall as \( 1/r^2 \).
The trick to simplifying this differential equation for the inverse square problems is to make a substitution, \( u\equiv 1/r \), and rewrite the differential equation for \( u(\theta) \). $$ \begin{eqnarray} r&=&1/u,\\ \nonumber \frac{dr}{d\theta}&=&-\frac{1}{u^2}\frac{du}{d\theta},\\ \nonumber \frac{d^2r}{d\theta^2}&=&\frac{2}{u^3}\left(\frac{du}{d\theta}\right)^2-\frac{1}{u^2}\frac{d^2u}{d\theta^2}. \end{eqnarray} $$
Plugging these expressions into Eq. \eqref{eq:rdotdot} gives an expression in terms of \( u \), \( du/d\theta \), and \( d^2u/d\theta^2 \). After some tedious algebra, $$ \begin{equation} \frac{d^2u}{d\theta^2}=-u-\frac{F m}{L^2u^2}. \label{_auto75} \end{equation} $$
For the attractive inverse square law force, \( F=-\alpha u^2 \), $$ \begin{equation} \frac{d^2u}{d\theta^2}=-u+\frac{m\alpha}{L^2}. \label{_auto76} \end{equation} $$
The solution has two arbitrary constants, \( A \) and \( \theta_0 \), $$ \begin{eqnarray} \label{eq:Ctrajectory} u&=&\frac{m\alpha}{L^2}+A\cos(\theta-\theta_0),\\ \nonumber r&=&\frac{1}{(m\alpha/L^2)+A\cos(\theta-\theta_0)}. \end{eqnarray} $$
The radius will be at a minimum when \( \theta=\theta_0 \) and at a maximum when \( \theta=\theta_0+\pi \). The constant \( A \) is related to the eccentricity of the orbit. When \( A=0 \) the radius is a constant \( r=L^2/(m\alpha) \), and the motion is circular. If one solved the expression \( mv^2/r=-\alpha/r^2 \) for a circular orbit, using the substitution \( v=L/(mr) \), one would reproduce the expression \( r=L^2/(m\alpha) \).
The form describing the elliptical trajectory in Eq. \eqref{eq:Ctrajectory} can be identified as an ellipse with one focus being the center of the ellipse by considering the definition of an ellipse as being the points such that the sum of the two distances between the two foci are a constant. Making that distance \( 2D \), the distance between the two foci as \( 2a \), and putting one focus at the origin, $$ \begin{eqnarray} 2D&=&r+\sqrt{(r\cos\theta-2a)^2+r^2\sin^2\theta},\\ \nonumber 4D^2+r^2-4Dr&=&r^2+4a^2-4ar\cos\theta,\\ \nonumber r&=&\frac{D^2-a^2}{D+a\cos\theta}=\frac{1}{D/(D^2-a^2)-a\cos\theta/(D^2-a^2)}. \end{eqnarray} $$
By inspection, this is the same form as Eq. \eqref{eq:Ctrajectory} with \( D/(D^2-a^2)=m\alpha/L^2 \) and \( a/(D^2-a^2)=A \).
Let us remind ourselves about what an ellipse is before we proceed.
import numpy as np
from matplotlib import pyplot as plt
from math import pi
u=1. #x-position of the center
v=0.5 #y-position of the center
a=2. #radius on the x-axis
b=1.5 #radius on the y-axis
t = np.linspace(0, 2*pi, 100)
plt.plot( u+a*np.cos(t) , v+b*np.sin(t) )
plt.grid(color='lightgray',linestyle='--')
plt.show()
The total energy of a particle is $$ \begin{eqnarray} E&=&U(r)+\frac{1}{2}mv_\theta^2+\frac{1}{2}m\dot{r}^2\\ \nonumber &=&U(r)+\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m\dot{r}^2\\ \nonumber &=&U(r)+\frac{L^2}{2mr^2}+\frac{1}{2}m\dot{r}^2. \end{eqnarray} $$
The second term then contributes to the energy like an additional repulsive potential. The term is sometimes referred to as the "centrifugal" potential, even though it is actually the kinetic energy of the angular motion. Combined with \( U(r) \), it is sometimes referred to as the "effective" potential, $$ \begin{eqnarray} U_{\rm eff}(r)&=&U(r)+\frac{L^2}{2mr^2}. \end{eqnarray} $$
Note that if one treats the effective potential like a real potential, one would expect to be able to generate an effective force, $$ \begin{eqnarray} F_{\rm eff}&=&-\frac{d}{dr}U(r) -\frac{d}{dr}\frac{L^2}{2mr^2}\\ \nonumber &=&F(r)+\frac{L^2}{mr^3}=F(r)+m\frac{v_\perp^2}{r}, \end{eqnarray} $$
which is indeed matches the form for \( m\ddot{r} \) in Eq. \eqref{eq:radialeqofmotion2}, which included the centrifugal force.
The following code plots this effective potential for a simple choice of parameters, with a standard gravitational potential \( -\alpha/r \). Here we have chosen \( L=m=\alpha=1 \).
# Common imports
import numpy as np
from math import *
import matplotlib.pyplot as plt
Deltax = 0.01
#set up arrays
xinitial = 0.3
xfinal = 5.0
alpha = 1.0 # spring constant
m = 1.0 # mass, you can change these
AngMom = 1.0 # The angular momentum
n = ceil((xfinal-xinitial)/Deltax)
x = np.zeros(n)
for i in range(n):
x[i] = xinitial+i*Deltax
V = np.zeros(n)
V = -alpha/x+0.5*AngMom*AngMom/(m*x*x)
# Plot potential
fig, ax = plt.subplots()
ax.set_xlabel('r[m]')
ax.set_ylabel('V[J]')
ax.plot(x, V)
fig.tight_layout()
plt.show()
Using the above parameters, we can now study the evolution of the system using for example the velocity Verlet method. This is done in the code here for an initial radius equal to the minimum of the potential well. We seen then that the radius is always the same and corresponds to a circle (the radius is always constant).
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
# Simple Gravitational Force -alpha/r
DeltaT = 0.01
#set up arrays
tfinal = 100.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
# Constants of the model, setting all variables to one for simplicity
alpha = 1.0
AngMom = 1.0 # The angular momentum
m = 1.0 # scale mass to one
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/m/alpha)
# Initial conditions
r0 = rmin
v0 = 0.0
r[0] = r0
v[0] = v0
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up acceleration
a = -alpha/(r[i]**2)+c1/(r[i]**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
anew = -alpha/(r[i+1]**2)+c1/(r[i+1]**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots(2,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Velocity')
ax[1].plot(t,v)
save_fig("RadialGVV")
plt.show()
Changing the value of the initial position to a value where the energy is positive, leads to an increasing radius with time, a so-called unbound orbit. Choosing on the other hand an initial radius that corresponds to a negative energy and different from the minimum value leads to a radius that oscillates back and forth between two values.
Consider a particle of mass \( m \) in a 2-dimensional harmonic oscillator with potential $$ U=\frac{1}{2}kr^2=\frac{1}{2}k(x^2+y^2). $$
If the orbit has angular momentum \( L \), we can find the radius and angular velocity of the circular orbit as well as the b) the angular frequency of small radial perturbations.
We consider the effective potential. The radius of a circular orbit is at the minimum of the potential (where the effective force is zero). The potential is plotted here with the parameters \( k=m=0.1 \) and \( L=1.0 \).
# Common imports
import numpy as np
from math import *
import matplotlib.pyplot as plt
Deltax = 0.01
#set up arrays
xinitial = 1.0
xfinal = 5.0
k = 0.1 # spring constant
m = 0.1 # mass, you can change these
AngMom = 1.0 # The angular momentum
n = ceil((xfinal-xinitial)/Deltax)
x = np.zeros(n)
for i in range(n):
x[i] = xinitial+i*Deltax
V = np.zeros(n)
V = 0.5*k*x*x+0.5*AngMom*AngMom/(m*x*x)
# Plot potential
fig, ax = plt.subplots()
ax.set_xlabel('r[m]')
ax.set_ylabel('V[J]')
ax.plot(x, V)
fig.tight_layout()
plt.show()
The effective potential looks like that of a harmonic oscillator for large \( r \), but for small \( r \), the centrifugal potential repels the particle from the origin. The combination of the two potentials has a minimum for at some radius \( r_{\rm min} \). $$ \begin{eqnarray*} 0&=&kr_{\rm min}-\frac{L^2}{mr_{\rm min}^3},\\ r_{\rm min}&=&\left(\frac{L^2}{mk}\right)^{1/4},\\ \dot{\theta}&=&\frac{L}{mr_{\rm min}^2}=\sqrt{k/m}. \end{eqnarray*} $$
For particles at \( r_{\rm min} \) with \( \dot{r}=0 \), the particle does not accelerate and \( r \) stays constant, i.e. a circular orbit. The radius of the circular orbit can be adjusted by changing the angular momentum \( L \).
For the above parameters this minimum is at \( r_{\rm min}=1 \).
Now consider small vibrations about \( r_{\rm min} \). The effective spring constant is the curvature of the effective potential. $$ \begin{eqnarray*} k_{\rm eff}&=&\left.\frac{d^2}{dr^2}U_{\rm eff}(r)\right|_{r=r_{\rm min}}=k+\frac{3L^2}{mr_{\rm min}^4}\\ &=&4k,\\ \omega&=&\sqrt{k_{\rm eff}/m}=2\sqrt{k/m}=2\dot{\theta}. \end{eqnarray*} $$
Here, the second step used the result of the last step from part (a). Because the radius oscillates with twice the angular frequency, the orbit has two places where \( r \) reaches a minimum in one cycle. This differs from the inverse-square force where there is one minimum in an orbit. One can show that the orbit for the harmonic oscillator is also elliptical, but in this case the center of the potential is at the center of the ellipse, not at one of the foci.
The solution is also simple to write down exactly in Cartesian coordinates. The \( x \) and \( y \) equations of motion separate, $$ \begin{eqnarray*} \ddot{x}&=&-kx,\\ \ddot{y}&=&-ky. \end{eqnarray*} $$
So the general solution can be expressed as $$ \begin{eqnarray*} x&=&A\cos\omega_0 t+B\sin\omega_0 t,\\ y&=&C\cos\omega_0 t+D\sin\omega_0 t. \end{eqnarray*} $$
The code here finds the solution for \( x \) and \( y \) using the code we developed in homework 4.
DeltaT = 0.01
#set up arrays
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
radius = np.zeros(n)
# Constants of the model
k = 0.1 # spring constant
m = 0.1 # mass, you can change these
omega02 = sqrt(k/m) # Frequency
AngMom = 1.0 # The angular momentum
rmin = (AngMom*AngMom/k/m)**0.25
# Initial conditions as compact 2-dimensional arrays
#x0 =rmin*0.5; y0 = sqrt(rmin*rmin-x0*x0)
x0 = 1.0; y0= 1.0
r0 = np.array([x0,y0])
v0 = np.array([0.0,0.0])
r[0] = r0
v[0] = v0
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up the acceleration
a = -r[i]*omega02
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
anew = -r[i+1]*omega02
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
# Plot position as function of time
radius = np.sqrt(r[:,0]**2+r[:,1]**2)
fig, ax = plt.subplots(3,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius squared')
ax[0].plot(t,r[:,0]**2+r[:,1]**2)
ax[1].set_xlabel('time')
ax[1].set_ylabel('x position')
ax[1].plot(t,r[:,0])
ax[2].set_xlabel('time')
ax[2].set_ylabel('y position')
ax[2].plot(t,r[:,1])
fig.tight_layout()
save_fig("2DimHOVV")
plt.show()
With some work using double angle formulas, one can calculate $$ \begin{eqnarray*} r^2&=&x^2+y^2\\ \nonumber &=&(A^2+C^2)\cos^2(\omega_0t)+(B^2+D^2)\sin^2\omega_0t+(AB+CD)\cos(\omega_0t)\sin(\omega_0t)\\ \nonumber &=&\alpha+\beta\cos 2\omega_0 t+\gamma\sin 2\omega_0 t,\\ \alpha&=&\frac{A^2+B^2+C^2+D^2}{2},~~\beta=\frac{A^2-B^2+C^2-D^2}{2},~~\gamma=AB+CD,\\ r^2&=&\alpha+(\beta^2+\gamma^2)^{1/2}\cos(2\omega_0 t-\delta),~~~\delta=\arctan(\gamma/\beta), \end{eqnarray*} $$
and see that radius oscillates with frequency \( 2\omega_0 \). The factor of two comes because the oscillation \( x=A\cos\omega_0t \) has two maxima for \( x^2 \), one at \( t=0 \) and one a half period later.
The following code shows first how we can solve this problem using the radial degrees of freedom only.
DeltaT = 0.01
#set up arrays
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
E = np.zeros(n)
# Constants of the model
AngMom = 1.0 # The angular momentum
m = 0.1
k = 0.1
omega02 = k/m
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/k/m)**0.25
# Initial conditions
r0 = rmin
v0 = 0.0
r[0] = r0
v[0] = v0
E[0] = 0.5*m*v0*v0+0.5*k*r0*r0+0.5*c2/(r0*r0)
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up acceleration
a = -r[i]*omega02+c1/(r[i]**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
anew = -r[i+1]*omega02+c1/(r[i+1]**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
E[i+1] = 0.5*m*v[i+1]*v[i+1]+0.5*k*r[i+1]*r[i+1]+0.5*c2/(r[i+1]*r[i+1])
# Plot position as function of time
fig, ax = plt.subplots(2,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Energy')
ax[1].plot(t,E)
save_fig("RadialHOVV")
plt.show()
The effective force can be extracted from the effective potential, \( U_{\rm eff} \). Beginning from the equations of motion, Eq. \eqref{eq:radialeqofmotion}, for \( r \), $$ \begin{eqnarray} m\ddot{r}&=&F+\frac{L^2}{mr^3}\\ \nonumber &=&F_{\rm eff}\\ \nonumber &=&-\partial_rU_{\rm eff},\\ \nonumber F_{\rm eff}&=&-\partial_r\left[U(r)+(L^2/2mr^2)\right]. \end{eqnarray} $$
For a circular orbit, the radius must be fixed as a function of time, so one must be at a maximum or a minimum of the effective potential. However, if one is at a maximum of the effective potential the radius will be unstable. For the attractive Coulomb force the effective potential will be dominated by the \( -\alpha/r \) term for large \( r \) because the centrifugal part falls off more quickly, \( \sim 1/r^2 \). At low \( r \) the centrifugal piece wins and the effective potential is repulsive. Thus, the potential must have a minimum somewhere with negative potential. The circular orbits are then stable to perturbation.
The effective potential is sketched for two cases, a \( 1/r \) attractive potential and a \( 1/r^3 \) attractive potential. The \( 1/r \) case has a stable minimum, whereas the circular orbit in the \( 1/r^3 \) case is unstable.
If one considers a potential that falls as \( 1/r^3 \), the situation is reversed and the point where \( \partial_rU \) disappears will be a local maximum rather than a local minimum.  The circular orbit is then at a maximum of the effective potential and the orbits are unstable. It is the clear that for potentials that fall as \( r^n \), that one must have \( n>-2 \) for the orbits to be stable.
Consider a potential \( U(r)=\beta r \). For a particle of mass \( m \) with angular momentum \( L \), find the angular frequency of a circular orbit. Then find the angular frequency for small radial perturbations.
For the circular orbit you search for the position \( r_{\rm min} \) where the effective potential is minimized, $$ \begin{eqnarray*} \partial_r\left\{\beta r+\frac{L^2}{2mr^2}\right\}&=&0,\\ \beta&=&\frac{L^2}{mr_{\rm min}^3},\\ r_{\rm min}&=&\left(\frac{L^2}{\beta m}\right)^{1/3},\\ \dot{\theta}&=&\frac{L}{mr_{\rm min}^2}=\frac{\beta^{2/3}}{(mL)^{1/3}} \end{eqnarray*} $$
Now, we can find the angular frequency of small perturbations about the circular orbit. To do this we find the effective spring constant for the effective potential, $$ \begin{eqnarray*} k_{\rm eff}&=&\partial_r^2 \left.U_{\rm eff}\right|_{r_{\rm min}}\\ &=&\frac{3L^2}{mr_{\rm min}^4},\\ \omega&=&\sqrt{\frac{k_{\rm eff}}{m}}\\ &=&\frac{\beta^{2/3}}{(mL)^{1/3}}\sqrt{3}. \end{eqnarray*} $$
If the two frequencies, \( \dot{\theta} \) and \( \omega \), differ by an integer factor, the orbit's trajectory will repeat itself each time around. This is the case for the inverse-square force, \( \omega=\dot{\theta} \), and for the harmonic oscillator, \( \omega=2\dot{\theta} \). In this case, \( \omega=\sqrt{3}\dot{\theta} \), and the angles at which the maxima and minima occur change with each orbit.
The code example here is meant to illustrate how we can make a plot of the final orbit. We solve the equations in polar coordinates (the example here uses the minimum of the potential as initial value) and then we transform back to cartesian coordinates and plot \( x \) versus \( y \). We see that we get a perfect circle when we place ourselves at the minimum of the potential energy, as expected.
# Simple Gravitational Force -alpha/r
DeltaT = 0.01
#set up arrays
tfinal = 8.0
n = ceil(tfinal/DeltaT)
# set up arrays for t, v and r
t = np.zeros(n)
v = np.zeros(n)
r = np.zeros(n)
phi = np.zeros(n)
x = np.zeros(n)
y = np.zeros(n)
# Constants of the model, setting all variables to one for simplicity
alpha = 1.0
AngMom = 1.0 # The angular momentum
m = 1.0 # scale mass to one
c1 = AngMom*AngMom/(m*m)
c2 = AngMom*AngMom/m
rmin = (AngMom*AngMom/m/alpha)
# Initial conditions, place yourself at the potential min
r0 = rmin
v0 = 0.0 # starts at rest
r[0] = r0
v[0] = v0
phi[0] = 0.0
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up acceleration
a = -alpha/(r[i]**2)+c1/(r[i]**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
anew = -alpha/(r[i+1]**2)+c1/(r[i+1]**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
phi[i+1] = t[i+1]*c2/(r0**2)
# Find cartesian coordinates for easy plot
x = r*np.cos(phi)
y = r*np.sin(phi)
fig, ax = plt.subplots(3,1)
ax[0].set_xlabel('time')
ax[0].set_ylabel('radius')
ax[0].plot(t,r)
ax[1].set_xlabel('time')
ax[1].set_ylabel('Angle $\cos{\phi}$')
ax[1].plot(t,np.cos(phi))
ax[2].set_ylabel('y')
ax[2].set_xlabel('x')
ax[2].plot(x,y)
save_fig("Phasespace")
plt.show()
Try to change the initial value for \( r \) and see what kind of orbits you get. In order to test different energies, it can be useful to look at the plot of the effective potential discussed above.
However, for orbits different from a circle the above code would need modifications in order to allow us to display say an ellipse. For the latter, it is much easier to run our code in cartesian coordinates, as done here. In this code we test also energy conservation and see that it is conserved to numerical precision. The code here is a simple extension of the code we developed for homework 4.
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
DeltaT = 0.01
#set up arrays
tfinal = 10.0
n = ceil(tfinal/DeltaT)
# set up arrays
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
E = np.zeros(n)
# Constants of the model
m = 1.0 # mass, you can change these
alpha = 1.0
# Initial conditions as compact 2-dimensional arrays
x0 = 0.5; y0= 0.
r0 = np.array([x0,y0])
v0 = np.array([0.0,1.0])
r[0] = r0
v[0] = v0
rabs = sqrt(sum(r[0]*r[0]))
E[0] = 0.5*m*(v[0,0]**2+v[0,1]**2)-alpha/rabs
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up the acceleration
rabs = sqrt(sum(r[i]*r[i]))
a = -alpha*r[i]/(rabs**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
rabs = sqrt(sum(r[i+1]*r[i+1]))
anew = -alpha*r[i+1]/(rabs**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
E[i+1] = 0.5*m*(v[i+1,0]**2+v[i+1,1]**2)-alpha/rabs
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots(3,1)
ax[0].set_ylabel('y')
ax[0].set_xlabel('x')
ax[0].plot(r[:,0],r[:,1])
ax[1].set_xlabel('time')
ax[1].set_ylabel('y position')
ax[1].plot(t,r[:,0])
ax[2].set_xlabel('time')
ax[2].set_ylabel('y position')
ax[2].plot(t,r[:,1])
fig.tight_layout()
save_fig("2DimGravity")
plt.show()
print(E)
Scattering experiments don't measure entire trajectories. For elastic collisions, they measure the distribution of final scattering angles at best. Most experiments use targets thin enough so that the number of scatterings is typically zero or one. The cross section, \( \sigma \), describes the cross-sectional area for particles to scatter with an individual target atom or nucleus. Cross section measurements form the basis for MANY fields of physics. BThe cross section, and the differential cross section, encapsulates everything measurable for a collision where all that is measured is the final state, e.g. the outgoing particle had momentum \( \boldsymbol{p}_f \). y studying cross sections, one can infer information about the potential interaction between the two particles. Inferring, or constraining, the potential from the cross section is a classic {\it inverse} problem. Collisions are either elastic or inelastic. Elastic collisions are those for which the two bodies are in the same internal state before and after the collision. If the collision excites one of the participants into a higher state, or transforms the particles into different species, or creates additional particles, the collision is inelastic. Here, we consider only elastic collisions.
For Coulomb forces, the cross section is infinite because the range of the Coulomb force is infinite, but for interactions such as the strong interaction in nuclear or particle physics, there is no long-range force and cross-sections are finite. Even for Coulomb forces, the part of the cross section that corresponds to a specific scattering angle, \( d\sigma/d\Omega \), which is a function of the scattering angle \( \theta_s \) is still finite.
If a particle travels through a thin target, the chance the particle scatters is \( P_{\rm scatt}=\sigma dN/dA \), where \( dN/dA \) is the number of scattering centers per area the particle encounters. If the density of the target is \( \rho \) particles per volume, and if the thickness of the target is \( t \), the areal density (number of target scatterers per area) is \( dN/dA=\rho t \). Because one wishes to quantify the collisions independently of the target, experimentalists measure scattering probabilities, then divide by the areal density to obtain cross-sections, $$ \begin{eqnarray} \sigma=\frac{P_{\rm scatt}}{dN/dA}. \end{eqnarray} $$
Instead of merely stating that a particle collided, one can measure the probability the particle scattered by a given angle. The scattering angle \( \theta_s \) is defined so that at zero the particle is unscattered and at \( \theta_s=\pi \) the particle is scattered directly backward. Scattering angles are often described in the center-of-mass frame, but that is a detail we will neglect for this first discussion, where we will consider the scattering of particles moving classically under the influence of fixed potentials \( U(\boldsymbol{r}) \). Because the distribution of scattering angles can be measured, one expresses the differential cross section, $$ \begin{equation} \frac{d^2\sigma}{d\cos\theta_s~d\phi}. \label{_auto77} \end{equation} $$
Usually, the literature expresses differential cross sections as $$ \begin{equation} d\sigma/d\Omega=\frac{d\sigma}{d\cos\theta d\phi}=\frac{1}{2\pi}\frac{d\sigma}{d\cos\theta}, \label{_auto78} \end{equation} $$
where the last equivalency is true when the scattering does not depend on the azimuthal angle \( \phi \), as is the case for spherically symmetric potentials.
The differential solid angle \( d\Omega \) can be thought of as the area subtended by a measurement, \( dA_d \), divided by \( r^2 \), where \( r \) is the distance to the detector, $$ \begin{eqnarray} dA_d=r^2 d\Omega. \end{eqnarray} $$
With this definition \( d\sigma/d\Omega \) is independent of the distance from which one places the detector, or the size of the detector (as long as it is small).
Differential scattering cross sections are calculated by assuming a random distribution of impact parameters \( b \). These represent the distance in the \( xy \) plane for particles moving in the \( z \) direction relative to the scattering center. An impact parameter \( b=0 \) refers to being aimed directly at the target's center. The impact parameter describes the transverse distance from the \( z=0 \) axis for the trajectory when it is still far away from the scattering center and has not yet passed it. The differential cross section can be expressed in terms of the impact parameter, $$ \begin{equation} d\sigma=2\pi bdb, \label{_auto79} \end{equation} $$
which is the area of a thin ring of radius \( b \) and thickness \( db \). In classical physics, one can calculate the trajectory given the incoming kinetic energy \( E \) and the impact parameter if one knows the mass and potential. From the trajectory, one then finds the scattering angle \( \theta_s(b) \). The differential cross section is then $$ \begin{equation} \frac{d\sigma}{d\Omega}=\frac{1}{2\pi}\frac{d\sigma}{d\cos\theta_s}=b\frac{db}{d\cos\theta_s}=\frac{b}{(d/db)\cos\theta_s(b)}. \label{_auto80} \end{equation} $$
Typically, one would calculate \( \cos\theta_s \) and \( (d/db)\cos\theta_s \) as functions of \( b \). This is sufficient to plot the differential cross section as a function of \( \theta_s \).
The total cross section is $$ \begin{equation} \sigma_{\rm tot}=\int d\Omega\frac{d\sigma}{d\Omega}=2\pi\int d\cos\theta_s~\frac{d\sigma}{d\Omega}. \label{_auto81} \end{equation} $$
Even if the total cross section is infinite, e.g. Coulomb forces, one can still have a finite differential cross section as we will see later on.
An asteroid of mass \( m \) and kinetic energy \( E \) approaches a planet of radius \( R \) and mass \( M \). What is the cross section for the asteroid to impact the planet?
allows one to solve for \( b_{\rm max} \), $$ \begin{eqnarray*} b_{\rm max}&=&R\frac{p_f}{p_0}\\ &=&R\frac{\sqrt{2m(E+GMm/R)}}{\sqrt{2mE}}\\ \sigma_{\rm impact}&=&\pi R^2\frac{E+GMm/R}{E}. \end{eqnarray*} $$
This refers to the calculation of \( d\sigma/d\Omega \) due to an inverse square force, \( F_{12}=\pm\alpha/r^2 \) for repulsive/attractive interaction. Rutherford compared the scattering of \( \alpha \) particles ($^4$He nuclei) off of a nucleus and found the scattering angle at which the formula began to fail. This corresponded to the impact parameter for which the trajectories would strike the nucleus. This provided the first measure of the size of the atomic nucleus. At the time, the distribution of the positive charge (the protons) was considered to be just as spread out amongst the atomic volume as the electrons. After Rutherford's experiment, it was clear that the radius of the nucleus tended to be roughly 4 orders of magnitude smaller than that of the atom, which is less than the size of a football relative to Spartan Stadium.
The incoming and outgoing angles of the trajectory are at \( \pm\theta' \). They are related to the scattering angle by \( 2\theta'=\pi+\theta_s \).
In order to calculate differential cross section, we must find how the impact parameter is related to the scattering angle. This requires analysis of the trajectory. We consider our previous expression for the trajectory where we derived the elliptic form for the trajectory, Eq. \eqref{eq:Ctrajectory}. For that case we considered an attractive force with the particle's energy being negative, i.e. it was bound. However, the same form will work for positive energy, and repulsive forces can be considered by simple flipping the sign of \( \alpha \). For positive energies, the trajectories will be hyperbolas, rather than ellipses, with the asymptotes of the trajectories representing the directions of the incoming and outgoing tracks. Rewriting Eq. \eqref{eq:Ctrajectory}, $$ \begin{equation}label{eq:ruthtraj} r=\frac{1}{\frac{m\alpha}{L^2}+A\cos\theta}. \end{equation} $$
Once \( A \) is large enough, which will happen when the energy is positive, the denominator will become negative for a range of \( \theta \). This is because the scattered particle will never reach certain angles. The asymptotic angles \( \theta' \) are those for which the denominator goes to zero, $$ \begin{equation} \cos\theta'=-\frac{m\alpha}{AL^2}. \label{_auto82} \end{equation} $$
The trajectory's point of closest approach is at \( \theta=0 \) and the two angles \( \theta' \), which have this value of \( \cos\theta' \), are the angles of the incoming and outgoing particles. From Fig (to come), one can see that the scattering angle \( \theta_s \) is given by, $$ \begin{eqnarray} \label{eq:sthetover2} 2\theta'-\pi&=&\theta_s,~~~\theta'=\frac{\pi}{2}+\frac{\theta_s}{2},\\ \nonumber \sin(\theta_s/2)&=&-\cos\theta'\\ \nonumber &=&\frac{m\alpha}{AL^2}. \end{eqnarray} $$
Now that we have \( \theta_s \) in terms of \( m,\alpha,L \) and \( A \), we wish to re-express \( L \) and \( A \) in terms of the impact parameter \( b \) and the energy \( E \). This will set us up to calculate the differential cross section, which requires knowing \( db/d\theta_s \). It is easy to write the angular momentum as $$ \begin{equation} L^2=p_0^2b^2=2mEb^2. \label{_auto83} \end{equation} $$
Finding \( A \) is more complicated. To accomplish this we realize that the point of closest approach occurs at \( \theta=0 \), so from Eq. \eqref{eq:ruthtraj} $$ \begin{eqnarray} \label{eq:rminofA} \frac{1}{r_{\rm min}}&=&\frac{m\alpha}{L^2}+A,\\ \nonumber A&=&\frac{1}{r_{\rm min}}-\frac{m\alpha}{L^2}. \end{eqnarray} $$
Next, \( r_{\rm min} \) can be found in terms of the energy because at the point of closest approach the kinetic energy is due purely to the motion perpendicular to \( \hat{r} \) and $$ \begin{equation} E=-\frac{\alpha}{r_{\rm min}}+\frac{L^2}{2mr_{\rm min}^2}. \label{_auto84} \end{equation} $$
One can solve the quadratic equation for \( 1/r_{\rm min} \), $$ \begin{equation} \frac{1}{r_{\rm min}}=\frac{m\alpha}{L^2}+\sqrt{(m\alpha/L^2)^2+2mE/L^2}. \label{_auto85} \end{equation} $$
We can plug the expression for \( r_{\rm min} \) into the expression for \( A \), Eq. \eqref{eq:rminofA}, $$ \begin{equation} A=\sqrt{(m\alpha/L^2)^2+2mE/L^2}=\sqrt{(\alpha^2/(4E^2b^4)+1/b^2} \label{_auto86} \end{equation} $$
Finally, we insert the expression for \( A \) into that for the scattering angle, Eq. \eqref{eq:sthetover2}, $$ \begin{eqnarray} \label{eq:scattangle} \sin(\theta_s/2)&=&\frac{m\alpha}{AL^2}\\ \nonumber &=&\frac{a}{\sqrt{a^2+b^2}}, ~~a\equiv \frac{\alpha}{2E} \end{eqnarray} $$
The differential cross section can now be found by differentiating the expression for \( \theta_s \) with \( b \), $$ \begin{eqnarray} \label{eq:rutherford} \frac{1}{2}\cos(\theta_s/2)d\theta_s&=&\frac{ab~db}{(a^2+b^2)^{3/2}}=\frac{bdb}{a^2}\sin^3(\theta_s/2),\\ \nonumber d\sigma&=&2\pi bdb=\frac{\pi a^2}{\sin^3(\theta_s/2)}\cos(\theta_s/2)d\theta_s\\ \nonumber &=&\frac{\pi a^2}{2\sin^4(\theta_s/2)}\sin\theta_s d\theta_s\\ \nonumber \frac{d\sigma}{d\cos\theta_s}&=&\frac{\pi a^2}{2\sin^4(\theta_s/2)},\\ \nonumber \frac{d\sigma}{d\Omega}&=&\frac{a^2}{4\sin^4(\theta_s/2)}. \end{eqnarray} $$
where \( a= \alpha/2E \). This the Rutherford formula for the differential cross section. It diverges as \( \theta_s\rightarrow 0 \) because scatterings with arbitrarily large impact parameters still scatter to arbitrarily small scattering angles. The expression for \( d\sigma/d\Omega \) is the same whether the interaction is positive or negative.
Consider a particle of mass \( m \) and charge \( z \) with kinetic energy \( E \) (Let it be the center-of-mass energy) incident on a heavy nucleus of mass \( M \) and charge \( Z \) and radius \( R \). Find the angle at which the Rutherford scattering formula breaks down.
Let \( \alpha=Zze^2/(4\pi\epsilon_0) \). The scattering angle in Eq. \eqref{eq:scattangle} is $$ \sin(\theta_s/2)=\frac{a}{\sqrt{a^2+b^2}}, ~~a\equiv \frac{\alpha}{2E}. $$
The impact parameter \( b \) for which the point of closest approach equals \( R \) can be found by using angular momentum conservation, $$ \begin{eqnarray*} p_0b&=&b\sqrt{2mE}=Rp_f=R\sqrt{2m(E-\alpha/R)},\\ b&=&R\frac{\sqrt{2m(E-\alpha/R)}}{\sqrt{2mE}}\\ &=&R\sqrt{1-\frac{\alpha}{ER}}. \end{eqnarray*} $$
Putting these together $$ \theta_s=2\sin^{-1}\left\{ \frac{a}{\sqrt{a^2+R^2(1-\alpha/(RE))}} \right\},~~~a=\frac{\alpha}{2E}. $$
It was from this departure of the experimentally measured \( d\sigma/d\Omega \) from the Rutherford formula that allowed Rutherford to infer the radius of the gold nucleus, \( R \).
Just like electrodynamics, one can define "fields", which for a small additional mass \( m \) are the force per mass and the additional potential energy per mass. The {\it gravitational field} related to the force has dimensions of force per mass, or acceleration, and can be labeled \( \boldsymbol{g}(\boldsymbol{r}) \). The potential energy per mass has dimensions of energy per mass. This is analogous to the electromagnetic potential, which is the potential energy per charge, and the electric field which is the force per charge.
Because the field \( \boldsymbol{g} \) obeys the same inverse square law for a point mass as the electric field does for a point charge, the gravitational field also satisfies a version of Gauss's law, $$ \begin{equation} \label{eq:GravGauss} \oint d\boldsymbol{A}\cdot\boldsymbol{g}=-4\pi GM_{\rm inside}. \end{equation} $$
Here, \( M_{\rm inside} \) is the net mass inside a closed area.
Gauss's law can be understood by considering a nozzle that sprays paint in all directions uniformly from a point source. Let \( B \) be the number of gallons per minute of paint leaving the nozzle. If the nozzle is at the center of a sphere of radius \( r \), the paint per square meter per minute that is deposited on some part of the sphere is $$ \begin{eqnarray} F(r)&=&\frac{B}{4\pi r^2}. \end{eqnarray} $$
Now, let \( F \) also be assigned a direction, so that it becomes a vector pointing along the direction of the flying paint. For any surface that surrounds the nozzle, not necessarily a sphere, one can state that $$ \begin{eqnarray} \label{eq:paint} \oint \boldsymbol{dA}\cdot\boldsymbol{F}&=&B, \end{eqnarray} $$
regardless of the shape of the surface. This follows because the rate at which paint is deposited on the surface should equal the rate at which it leaves the nozzle. The dot product ensures that only the component of \( \boldsymbol{F} \) into the surface contributes to the deposition of paint. Similarly, if \( \boldsymbol{F} \) is any radial inverse-square forces, that falls as \( B/(4\pi r^2) \), then one can apply Eq. \eqref{eq:paint}. For gravitational fields, \( B/(4\pi) \) is replaced by \( GM \), and one quickly "derives" Gauss's law for gravity, Eq. \eqref{eq:GravGauss}.
Consider Earth to have its mass \( M \) uniformly distributed in a sphere of radius \( R \). Find the magnitude of the gravitational acceleration as a function of the radius \( r \) in terms of the acceleration of gravity at the surface \( g(R) \). Assume \( r < R \), i.e. you are inside the surface.
{\bf Solution}: Take the ratio of Eq. \eqref{eq:GravGauss} for two radii, \( R \) and \( r < R \), $$ \begin{eqnarray*} \frac{4\pi r^2 g(r)}{4\pi R^2 g(R)}&=&\frac{4\pi GM_{\rm inside~r}}{4\pi GM_{\rm inside~R}}\\ \nonumber &=&\frac{r^3}{R^3}\\ \nonumber g(r)&=&g(R)\frac{r}{R}~. \end{eqnarray*} $$ The potential energy per mass is similar conceptually to the voltage, or electric potential energy per charge, that was studied in electromagnetism, if \( V\equiv U/m \), \( \boldsymbol{g}=-\nabla V \).
Consider a spherical planet of radius \( r \) a distance \( D \) from another body of mass \( M \). The magnitude of the force due to \( M \) on an small object of mass \( \delta m \) on surface of the planet can be calculated by performing a Taylor expansion about the center of the spherical planet. $$ \begin{equation} F=-\frac{GM\delta m}{D^2}+2\frac{GM\delta m}{D^3}\Delta D+\cdots \label{_auto87} \end{equation} $$
If the \( z \) direction points toward the large object, \( \Delta D \) can be referred to as \( z \). In the accelerating frame of an observer at the center of the planet, $$ \begin{equation} \delta m\frac{d^2 z}{dt^2}=F-\delta ma'+{\rm other~forces~acting~on~} \delta m, \label{_auto88} \end{equation} $$
where \( a' \) is the acceleration of the observer. Because \( \delta ma' \) equals the gravitational force on \( \delta m \) if it were located at the planet's center, one can write $$ \begin{equation} m\frac{d^2z}{dt^2}=2\frac{GM\delta m}{D^3}z+{\rm other~forces~acting~on~}\delta m. \label{_auto89} \end{equation} $$
Here the other forces could represent the forces acting on \( \delta m \) from the spherical planet such as the gravitational force or the contact force with the surface. If \( \theta \) is the angle w.r.t. the \( z \) axis, the effective force acting on \( \delta m \) is $$ \begin{equation} F_{\rm eff}\approx 2\frac{GM\delta m}{D^3}r\cos\theta\hat{z}+{\rm other~forces~acting~on~}\delta m. \label{_auto90} \end{equation} $$
This first force is the "tidal" force. It pulls objects outward from the center of the object. If the object were covered with water, it would distort the objects shape so that the shape would be elliptical, stretched out along the axis pointing toward the large mass \( M \). The force is always along (either parallel or antiparallel to) the \( \hat{z} \) direction.
Consider the Earth to be a sphere of radius \( R \) covered with water,
with the gravitational acceleration at the surface noted by \( g \). Now
assume that a distant body provides an additional constant
gravitational acceleration \( \boldsymbol{a} \) pointed along the \( z \) axis. Find
the distortion of the radius as a function of \( \theta \). Ignore
planetary rotation and assume \( a <
Now consider that the field is no longer constant, but that instead \( a=-kz \) with \( |kR| <
Here, the potential due to the external field is \( (1/2)kz^2 \) so that \( -\nabla U=-kz \). One now needs to solve for \( h(\theta) \). Absorbing all the constant terms from both sides of the equation into one constant \( C \), and because both \( h \) and \( kR \) are small, we can through away terms of order \( h^2 \) or \( kRh \). This gives
$$
\begin{eqnarray*}
gh(\theta)-\frac{1}{2}kR^2\cos^2\theta&=&C,\\
h(\theta)&=&\frac{C}{g}+\frac{1}{2g}kR^2\cos^2\theta,\\
h(\theta)&=&\frac{1}{2g}kR^2(\cos^2\theta-1/3).
\end{eqnarray*}
$$
The term with the factor of \( 1/3 \) replaced the constant and was chosen so that the average height of the water would be zero.
The Sun's mass is \( 27\times 10^6 \) the Moon's mass, but the Sun is 390 times further away from Earth as the Sun. What is ratio of the tidal force of the Sun to that of the Moon.
{\bf Solution}: The gravitational force due to an object \( M \) a distance \( D \) away goes as \( M/D^2 \), but the tidal force is only the difference of that force over a distance \( R \),
$$
F_{\rm tidal}\propto \frac{M}{D^3}R.
$$
Therefore the ratio of force is
$$
\begin{eqnarray*}
\frac{F_{\rm Sun's~tidal~force}}{F_{\rm Moon's~tidal~force}}
&=&\frac{M_{\rm sun}/D_{\rm sun}^3}{M_{\rm moon}/D_{\rm moon}^3}\\
&=&\frac{27\times 10^6}{390^3}=0.46.
\end{eqnarray*}
$$
The Moon more strongly affects tides than the Sun.
The calculus of variations involves
problems where the quantity to be minimized or maximized is an integral.
The usual minimization problem one faces involves taking a function
\( {\cal L}(x) \), then finding the single value \( x \) for which \( {\cal L} \)
is either a maximum or minimum. In multivariate calculus one also
learns to solve problems where you minimize for multiple variables,
\( {\cal L}(x_1,x_2,\cdots x_n) \), and finding the points \( (x_1\cdots
y_n) \) in an \( n \)-dimensional space that maximize or minimize the
function. Here, we consider what seems to be a much more ambitious
problem. Imagine you have a function \( {\cal L}(x(t),\dot{x}(t),t) \),
and you wish to find the extrema for an infinite number of values of
\( x \), i.e. \( x \) at each point \( t \). The function \( {\cal L} \) will not only
depend on \( x \) at each point \( t \), but also on the slope at each point,
plus an additional dependence on \( t \). Note we are NOT finding an
optimum value of \( t \), we are finding the set of optimum values of \( x \)
at each point \( t \), or equivalently, finding the function \( x(t) \).
One treats the function \( x(t) \) as being unknown while minimizing the action
$$
S=\int_{t_1}^{t_2}dt~{\cal L}(x(t),\dot{x}(t),t).
$$
Thus, we are minimizing \( S \) with respect to an infinite number of
values of \( x(t_i) \) at points \( t_i \). As an additional criteria, we will
assume that \( x(t_1) \) and \( x(t_2) \) are fixed, and that that we will
only consider variations of \( x \) between the boundaries. The dependence
on the derivative, \( \dot{x}=dx/dt \), is crucial because otherwise the
solution would involve simply finding the one value of \( x \) that
minimized \( {\cal L} \), and \( x(t) \) would equal a constant if there were no
explicit \( t \) dependence. Furthermore, \( x \) wouldn't need to be
continuous at the boundary.
In the general case we have an integral of the type
$$
S[q]= \int_{t_1}^{t_2} {\cal L}(q(t),\dot{q}(t),t)dt,
$$
where \( S \) is the quantity which is sought minimized or maximized. The
problem is that although \( {\cal L} \) is a function of the general variables
\( q(t),\dot{q}(t),t \) (note our change of variables), the exact dependence of \( q \) on \( t \) is not known.
This means again that even though the integral has fixed limits \( t_1 \)
and \( t_2 \), the path of integration is not known. In our case the unknown
quantities are the positions and general velocities of a given number
of objects and we wish to choose an integration path which makes the
functional \( S[q] \) stationary. This means that we want to find minima,
or maxima or saddle points. In physics we search normally for minima.
Our task is therefore to find the minimum of \( S[q] \) so that its
variation \( \delta S \) is zero subject to specific constraints. The
constraints can be treated via the technique of Lagrangian multipliers
as we will see below.
We assume the existence of an optimum path, that is a path for which
\( S[q] \) is stationary. There are infinitely many such paths. The
difference between two paths \( \delta q \) is called the variation of
\( q \).
We call the variation \( \eta(t) \) and it is scaled by a factor \( \alpha \).
The function \( \eta(t) \) is arbitrary except for
$$
\eta(t_1)=\eta(t_2)=0,
$$
and we assume that we can model the change in \( q \) as
$$
q(t,\alpha) = q(t)+\alpha\eta(t),
$$
and
$$
\delta q = q(t,\alpha) -q(t,0)=\alpha\eta(t).
$$
We choose \( q(t,\alpha=0) \) as the unkonwn path that will minimize \( S \). The value
\( q(t,\alpha\ne 0) \) describes a neighbouring path.
We have
$$
S[q(\alpha)]= \int_{t_1}^{t_2} {\cal L}(q(t,\alpha),\dot{q}(t,\alpha),t)dt.
$$
The condition for an extreme of
$$
S[q(\alpha)]= \int_{t_1}^{t_2} {\cal L}(q(t,\alpha),\dot{q}(t,\alpha),t)dt,
$$
is
$$
\left[\frac{\partial S[q(\alpha)]}{\partial t}\right]_{\alpha=0} =0.
$$
The \( \alpha \) dependence is contained in \( q(t,\alpha) \) and \( \dot{q}(t,\alpha) \) meaning that
$$
\left[\frac{\partial E[q(\alpha)]}{\partial \alpha}\right]=\int_{t_1}^{t_2} \left( \frac{\partial {\cal l}}{\partial q}\frac{\partial q}{\partial \alpha}+\frac{\partial {\cal L}}{\partial \dot{q}}\frac{\partial \dot{q}}{\partial \alpha}\right)dt.
$$
We have defined
$$
\frac{\partial q(x,\alpha)}{\partial \alpha}=\eta(x)
$$
and thereby
$$
\frac{\partial \dot{q}(t,\alpha)}{\partial \alpha}=\frac{d(\eta(t))}{dt}.
$$
Using
$$
\frac{\partial q(t,\alpha)}{\partial \alpha}=\eta(t),
$$
and
$$
\frac{\partial \dot{q}(t,\alpha)}{\partial \alpha}=\frac{d(\eta(t))}{dt},
$$
in the integral gives
$$
\left[\frac{\partial S[q(\alpha)]}{\partial \alpha}\right]=\int_{t_1}^{t_2} \left( \frac{\partial {\cal L}}{\partial q}\eta(t)+\frac{\partial {\cal L}}{\partial \dot{q}}\frac{d(\eta(t))}{dt}\right)dt.
$$
Integrating the second term by parts
$$
\int_{t_1}^{t_2} \frac{\partial {\cal L}}{\partial \dot{q}}\frac{d(\eta(t))}{dt}dt =\eta(t)\frac{\partial {\cal L}}{\partial \dot{q}}|_{t_1}^{t_2}-
\int_a^b \eta(t)\frac{d}{dx}\frac{\partial {\cal L}}{\partial \dot{q}}dt,
$$
and since the first term dissappears due to \( \eta(a)=\eta(b)=0 \), we obtain
$$
\left[\frac{\partial S[q(\alpha)]}{\partial \alpha}\right]=\int_{t_1}^{t_2} \left( \frac{\partial {\cal L}}{\partial q}-\frac{d}{dx}\frac{\partial {\cal L}}{\partial \dot{q}}
\right)\eta(t)dt=0,
$$
which can also be written as
$$
\left[\frac{\partial S[q(\alpha)]}{\partial \alpha}\right]_{\alpha=0}=\int_{t_1}^{t_2} \left( \frac{\partial {\cal L}}{\partial q}-\frac{d}{\
dx}\frac{\partial {\cal L}}{\partial \dot{q}}\right)\delta q(t)dt=\delta S = 0.
$$
The condition for a stationary value is thus a partial differential equation
$$
\frac{\partial {\cal L}}{\partial q}-\frac{d}{dx}\frac{\partial {\cal L}}{\partial \dot{q}}=0,
$$
known as the Euler-Lagrange equation.
Sometimes an auxiliary constraint is added to the problem (beyond
fixing the end poits \( y_1 \) and \( y_2 \)). Just ahead, we will work on the
example of a hanging chain. The shape of the curve minimizes the
potential energy, under the constraint of a fixed length of
chain. Before presenting such an example we first review the method of
Lagrange multipliers as a method for finding minima or maxima under
constraints.
Imagine a function \( f(x_1,x_2\cdots x_n) \) for which you wish to find
the minima. Additionally, you are given a constraint
$$
\begin{eqnarray}
C(x_1\cdots x_n)=0
\end{eqnarray}
$$
The usual condition for a a minimum is
$$
\begin{eqnarray}
\frac{\partial f}{\partial x_i}=0{\rm ,~~or~}\nabla f=0.
\end{eqnarray}
$$
which would be \( n \) equations for the \( n \) variables. The gradient of a
scalar is a vector, so you should think of \( \nabla \) as
\( \boldsymbol{\nabla} \). However, the solution will likely not satisfy the
constraint, i.e. the point at which \( f(x_1\cdots x_n) \) has an extrema,
may not be a point where \( C(x_1\cdots x_n)=0 \).
A necessary condition for the solution is that
$$
\begin{equation}
\nabla f\cdot\boldsymbol{\epsilon}=0,
\label{_auto91}
\end{equation}
$$
for any infinitesimal vector \( \boldsymbol{\epsilon} \) if \( \boldsymbol{\epsilon} \)
satisfies the condition
$$
\begin{equation}
\delta C=\nabla C\cdot\boldsymbol{\epsilon}=0.
\label{_auto92}
\end{equation}
$$
That is to say if I take a small step in a direction that doesn't
change the constraint, then \( f \) must not change if it is an
extrema. Not changing the constraint implies the step is orthogonal to
\( \nabla C \). As there are \( n \) dimensions of \( x \), the vector \( \nabla C \)
defines one direction, and \( \boldsymbol{\epsilon} \) can be in any of the \( n-1 \)
directions orthogonal to \( \nabla C \). If \( \nabla f\cdot\boldsymbol{\epsilon}=0 \)
for ANY of the \( n-1 \) directions of \( \boldsymbol{\epsilon} \) orthogonal to
\( \nabla C \), then
$$
\begin{equation}
\nabla f ~||~ \nabla C.
\label{_auto93}
\end{equation}
$$
Because the two vectors are parallel you can say there must exist some
constant \( \lambda \) such that
$$
\begin{equation}
\nabla(f-\lambda C)=0.
\label{_auto94}
\end{equation}
$$
Here, \( \lambda \) is known as a Lagrange multiplier. Satisfying
this equation is a necessary, but not a sufficient
condition. One could add a constant to the constraint and the gradient
would not change. One must find the correct value of \( \lambda \) that
satisfies the constraint \( C=0 \), rather than \( C= \) some other
constant. The strategy is then to solve
thw above equation then adjust \( \lambda \) until one
finds the \( x_1\cdots x_n \) that gives \( C(x_1\cdots x_n)=0 \).
The method of Lagrange multipliers is counter-intuitive to one's
intuition to use the constraint to reduce the dimensionality of the
problem. Normally, minimizing a function of \( n \) variables, leads to
\( n \) equations and \( n \) unknowns. A constraint could be used, by
substitution, to replace the \( n \) variables with \( n-1 \)
variables. Instead, we add an unknown parameter, \( \lambda \), and change
the equation to \( n+1 \) equations with \( n+1 \) unknowns, with the extra
unknown being the Lagrange multiplier \( \lambda \). Often, it is rather
easy to solve for \( x_1\cdots x_n \). Then one is left with the usually
difficult problem of finding \( \lambda \), often requiring the solution
of a transcendental equation.
Let us try to formalize this. We consider a function of three independent variables \( f(x,y,z) \) . For
the function \( f \) to be an extreme we have
$$
df=0.
$$
A necessary and sufficient condition is
$$
\frac{\partial f}{\partial x} =\frac{\partial f}{\partial y}=\frac{\partial f}{\partial z}=0,
$$
due to
$$
df = \frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz.
$$
In physical problems the variables \( x,y,z \) are often subject to constraints (in our case \( q \) and the orthogonality constraint)
so that they are no longer all independent. It is possible at least in principle to use each constraint to eliminate one variable
and to proceed with a new and smaller set of independent varables.
The use of so-called Lagrangian multipliers is an alternative technique when the elimination of
of variables is incovenient or undesirable. Assume that we have an equation of constraint on the variables \( x,y,z \)
$$
\phi(x,y,z) = 0,
$$
resulting in
$$
d\phi = \frac{\partial \phi}{\partial x}dx+\frac{\partial \phi}{\partial y}dy+\frac{\partial \phi}{\partial z}dz =0.
$$
Now we cannot set anymore
$$
\frac{\partial f}{\partial x} =\frac{\partial f}{\partial y}=\frac{\partial f}{\partial z}=0,
$$
if \( df=0 \) is wanted
because there are now only two independent variables! Assume \( x \) and \( y \) are the independent variables.
Then \( dz \) is no longer arbitrary.
However, we can add to
$$
df = \frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial y}dy+\frac{\partial f}{\partial z}dz,
$$
a multiplum of \( d\phi \), viz. \( \lambda d\phi \), resulting in
$$
df+\lambda d\phi = (\frac{\partial f}{\partial z}+\lambda\frac{\partial \phi}{\partial x})dx+(\frac{\partial f}{\partial y}+\lambda\frac{\partial \phi}{\partial y})dy+(\frac{\partial f}{\partial z}+\lambda\frac{\partial \phi}{\partial z})dz =0.
$$
Our multiplier is chosen so that
$$
\frac{\partial f}{\partial z}+\lambda\frac{\partial \phi}{\partial z} =0.
$$
However, we took \( dx \) and \( dy \) as to be arbitrary and thus we must have
$$
\frac{\partial f}{\partial x}+\lambda\frac{\partial \phi}{\partial x} =0,
$$
and
$$
\frac{\partial f}{\partial y}+\lambda\frac{\partial \phi}{\partial y} =0.
$$
When all these equations are satisfied, \( df=0 \). We have four
unknowns, \( x,y,z \) and \( \lambda \). Actually we want only \( x,y,z \),
\( \lambda \) need not to be determined, it is therefore often called
Lagrange's undetermined multiplier. If we have a set of constraints
\( \phi_k \) we have the equations
$$
\frac{\partial f}{\partial x_i}+\sum_k\lambda_k\frac{\partial \phi_k}{\partial x_i} =0.
$$
Consider a particle constrained to move along a path (like a bead
moving without friction on a wire) and you need to design a path from
\( x=y=0 \) to some final point \( x_f,y_f \). Assume there is a constant
force in the \( x \) direction, \( F_x=mg \). Design the path so that the time
the bead travels is a minimum.
The net time is
$$
T=\int \frac{d\ell}{v}=\int_0^{x_f} dx~\frac{\sqrt{1+y'^2}}{\sqrt{2gx}}={\rm minimum}.
$$
Here we made use of the fact that \( d\ell=\sqrt{dx^2+dy^2} \) and that
the velocity is determined by \( KE=mv^2/2=mgx \). The Euler equations can
be applied if you first define the function as
$$
\begin{eqnarray*}
f(y,y';x)&=&\frac{\sqrt{1+y'^2}}{\sqrt{x}}.
\end{eqnarray*}
$$
The equations are then
$$
\begin{eqnarray*}
\frac{d}{dx}\frac{\partial f}{\partial y'}&=&0.
\end{eqnarray*}
$$
The simplification ensued from \( f \) not having any dependence on \( y \). This yields the differential equation
$$
\begin{eqnarray}
\frac{y'}{x^{1/2}(1+y'^2)^{1/2}}&=&(2a)^{-1/2},
\end{eqnarray}
$$
because \( \partial f/\partial y' \) must be a constant, which with some
foresight we label \( (2a)^{-1/2} \). One can now solve for \( y' \),
$$
\begin{eqnarray*}
(y')^2&=&2ax(1+y'^2)\\
\nonumber
y'&=&\sqrt{\frac{x}{2a-x}},\\
\nonumber
y(t)&=&\int_0^x dx'~\frac{\sqrt{x'}dx'}{\sqrt{2a-x'}}=\int_0^x dx'~\frac{x'dx'}{\sqrt{2ax'-x'^2}}\\
\nonumber
&=&\frac{1}{2}\int_0^x\frac{(2x'-2a)dx'}{(2ax'-x'^2)^{1/2}}+a\int_0^x\frac{dx'}{\sqrt{2ax'-x'^2}}\\
\nonumber
&=&\frac{-1}{2}\int_0^{2ax-x^2}\frac{du}{\sqrt{u}}+a\int_0^x\frac{dx'}{\sqrt{a^2-(x'-a)^2}}\\
&=&-\sqrt{2ax-x^2}+a\cos^{-1}(1-x/a).
\end{eqnarray*}
$$
This turns out to be the equation for a {\it cycloid} or a {\it
brachiostone}. If you rolled a wheel of radius \( a \) down the \( y \) axis
and followed a point on the rim, it would trace out a cycloid. Here,
the constant \( a \) must be chosen to match the boundary condition,
\( y_2=y(x_2) \). You can see the textbook for more details, plus you get
a chance to work with cycloids in the exercises at the end of this
chapter.
As an example of using Lagrange multipliers for a standard
optimization formula we attempt to maximize the following function,
$$
F(x_1\cdots x_n)=-\sum_{i=1}^n x_i\ln(x_i),
$$
with respect to the \( n \) variables \( x_i \). With no constraints, each
\( x_i \) would maximize the function for
$$
\begin{eqnarray*}
\frac{d}{dx_j}~\left[-\sum_i x_i\ln(x_i)\right]&=&0\\
-\ln(x_j)-1&=&0,~~~~x_j=e^{-1}.
\end{eqnarray*}
$$
Now, we repeat the problem but with two constraints,
$$
\sum_ix_i=1~,~~~~\sum_ix_i\epsilon_i=E.
$$
Here, \( \epsilon_i \) and \( E \) are fixed constants. We go forward by
finding the extrema for
$$
\begin{eqnarray*}
G(x_1\cdots x_n)&=&F-\alpha\sum_i x_i-\beta\sum_i\epsilon_ix_i
=\sum_i \left\{-x_i\ln(x_i)-\alpha x_i-\beta\epsilon_ix_i\right\}.
\end{eqnarray*}
$$
There are two Lagrange multipliers, \( \alpha \) and \( \beta \),
corresponding to the two constraints. One then solves for the extrema
$$
\begin{eqnarray*}
\frac{d}{dx_j}G&=&0\\
&=&-\ln(x_j)-1-\alpha-\beta\epsilon_j,\\
x_j&=&\exp\left\{-1-\alpha-\beta\epsilon_j\right\}.
\end{eqnarray*}
$$
For any given \( \alpha \) and \( \beta \) this provides a solution for
constraining \( \sum_i x_i \) and \( \sum_i\epsilon_ix_i \) to some values,
just not the values of unity and \( E \) that you wish. One would then
have to search for the correct values by adjusting \( \alpha \) and
\( \beta \) until the constraint are actually matched by solving a
transcendental equation. Although this can be complicated, it is
certainly less expensive than searching over all \( N \) values of
\( x_i \). This particular example corresponds to maximizing the entropy
for a system, \( S=-\sum_i x_i\ln(x_i) \), where \( x_i \) is the probability
of the system being in a particular discrete level \( i \) that has energy
\( \epsilon_i \). One wishes to maximize the entropy subject to the
constraints that the probabilities sum to unity and the average energy
has some given value. The result that \( x_i\sim e^{-\beta\epsilon_i} \)
demonstrates the origin of the Boltzmann factor, with the inverse
temperature \( \beta=1/T \).
Lagrange multipliers also assist with the Euler-Lagrange equation. If
one breaks an interval \( x_1 < x < x_2 \) into a large number
\( n\rightarrow\infty \) points separated by \( dx \), the Euler-Lagrange
equation involves finding the \( n \) values \( y_i \) at each point so that
\( \sum_i dx f\left\{y_i,y'_i=(y_{i+1}-y_{i-1})/(2dx)\right\} \) is
maximized for some given function \( f \). If an additional auxiliary
constraint is added, also some function of the \( n \) values \( y_i \), one
can use the method of Lagrange multipliers. In the constraint can also
be written as some function of \( C(y_i,y'_i) \), then one simply adds a
term \( \lambda C(y,y') \) to the function \( f \) and uses the Euler-Lagrange
equation to find the extrema of.
$$
\begin{eqnarray}
J&=&\int_{x_1}^{x_2}dx~f\left\{y(t),y'(t),x\right\}-\lambda C\left\{y(t),y'(t),x\right\},
\end{eqnarray}
$$
the one difference being that
$$
\begin{equation}
f\left\{y(t),y'(t),x\right\}\rightarrow f\left\{y(t),y'(t),x\right\}-\lambda C\left\{y(t),y'(t),x\right\}
\label{_auto95}
\end{equation}
$$
Consider a chain of length \( L \) and mass per unit length \( \kappa \) that
hangs from point \( x=0,y=0 \) to point \( x_f,y_f \). The shape must minimize
the potential energy. Find general expressions for the shape in terms
of three constants which must be chosen to match \( y(0)=0, y(x_f)=y_f \)
and the fixed length. Equivalently, one finds the function \( y(t) \) that
provides an extrema for the integral,
One must minimize
$$
\int d\ell~\kappa gy-\lambda\int d\ell=
\int_0^{x_f} dx~\sqrt{1+y'^2}\kappa gy-\lambda \int_0^{x_f} dx\sqrt{1+y'^2}.
$$
Here \( \lambda \) is the Lagrange multiplier associated with constraining
the length of the chain. The constrained length \( L \) appears nowhere in
the expression. Instead, one solves for form of the answer, then
adjusts \( \lambda \) to give the correct length. For the purposes of the
Euler-Lagrange minimization one considers the function
$$
\begin{eqnarray}
f(y,y';x)&=&\kappa gy\sqrt{1+y'^2}-\lambda\sqrt{1+y'^2}.
\end{eqnarray}
$$
Because \( \lambda \) is an unknown constant and because minimizing a
function multiplied by a constant is the same as minimizing the
function, we can equivlently minimize the integral using the function
$$
\begin{eqnarray}
\tilde{f}(y,y';x)&=&y\sqrt{1+y'^2}-\tilde{\lambda}\sqrt{1+y'^2},\\
\nonumber
\tilde{\lambda}&\equiv&\frac{\lambda}{\kappa g}.
\end{eqnarray}
$$
The Euler-Lagrange equations then become
$$
\begin{eqnarray*}
\frac{d}{dx}\left\{
\frac{y'}{\sqrt{1+y'^2}}y-\tilde{\lambda}\frac{y'}{\sqrt{1+y'^2}}
\right\}&=&\sqrt{1+y'^2}.
\end{eqnarray*}
$$
Here, we will guess at the form of the solution,
$$
\begin{eqnarray*}
y'&=&\sinh[(x-x_0)/a],~~y=a\cosh[(x-x_0)/a]+y_0.
\end{eqnarray*}
$$
Plugging into the Euler-Lagange equations,
$$
\begin{eqnarray*}
\frac{d}{dx}\left\{(a\cosh[(x-x_0)/a]+y_0)\frac{\sinh[(x-x_0)/a]}{\cosh[(x-x_0)/a]}-\tilde{\lambda}\frac{\sinh[(x-x_0)/a]}{\cosh[(x-x_0)/a]}\right\}&=&\cosh[(x-x_0)/a],\\
\nonumber
\frac{d}{dx}\left\{(y_0-\tilde{\lambda})\tanh[(x-x_0)/a]\right\}=0.
\end{eqnarray*}
$$
This solution works if \( y_0=\tilde{\lambda} \). So the general form of
the solution is
$$
y=\tilde{\lambda}+a\cosh[(x-x_0)/a].
$$
One must find \( \tilde{\lambda} \), \( x_0 \) and \( a \) to satisfy three
conditions, \( y(x=0)=0 \), \( y(x=x_f)=y_f \) and that the length is \( L \). For
a hanging chain \( a \) is positive. A solution with negative \( a \) would
represent a maximum of the potential energy. A remarkable property of
the solution is that once you define the length and the end-point
positions \( y_1 \) and \( y_2 \), the solution does not depend on \( \kappa \) or
\( g \). Thus, the shape of the chain would be the same if you took it to
the moon. These solutions are known as catenaries.
Lagrangians represent a powerful method for solving problems that
would be nearly impossible by direct application of Newton's third
law, \( \boldsymbol{F}=m\boldsymbol{a} \). The method works well for problems where a
system is well described by a few \textit{generalized coordinates}. A
generalized coordinate might be the angle describing the position of a
pendulum. This one angle takes the place of using \( x \) and \( y \) to
describe the position of the pendulum, then applying a clumsy
constraint.
The Lagrangian equations of motion can be derived from a principle of
least action, where the action \( S \) is defined as
$$
\begin{equation}
S=\int dt~ L(q,\dot{q},t),
\label{_auto96}
\end{equation}
$$
where \( q \) is some coordinate that describes the orientation of a
system and the Lagrangian \( L \) is defined as
$$
\begin{equation}
L=T-U,
\label{_auto97}
\end{equation}
$$
the difference of the kinetic and potential energies. Minimizing the
action through the Euler-Lagrange equations gives the Lagrangian
equations of motion,
$$
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}=\frac{\partial L}{\partial q}.
\label{_auto98}
\end{equation}
$$
We begin with two simple examples, neither of which gains from the Lagrangian approach.
Consider a particle of mass \( m \) connected to a spring with stiffness \( k \). Derive the Lagrangian equations of motion.
$$
\begin{eqnarray*}
L&=&\frac{1}{2}m\dot{x}^2-\frac{1}{2}kx^2,\\
\frac{d}{dt}\frac{\partial L}{\partial \dot{x}}&=&\frac{\partial L}{\partial x},\\
m\ddot{x}&=&-kx.
\end{eqnarray*}
$$
Derive the Lagrangian equations of motion for a pendulum of mass \( m \)
and length \( \ell \).
$$
\begin{eqnarray*}
L&=&\frac{m}{2}\ell^2\dot{\theta}^2-mg\ell(1-\cos\theta),\\
\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta}}&=&\frac{\partial L}{\partial \theta},\\
m\ell^2\ddot{\theta}&=&-mg\ell\sin\theta,\\
\ddot{\theta}&=&-\frac{g}{\ell}\sin\theta,\\
\ddot{\theta}&\approx&-\frac{g}{\ell}\theta.
\end{eqnarray*}
$$
Lagrange's equations of motion can only be applied for the following conditions:
In the first line we used the fact that \( T \) does not depend on
\( x \). Continuing with taking the derivative of \( U \),
$$
\begin{eqnarray}
-\frac{\partial U}{\partial\dot{q}_j}&=&-\sum_i\frac{\partial U}{\partial x_i}\frac{\partial x_i}{\partial\dot{q}_j}=0.
\end{eqnarray}
$$
In the first line above we used the fact that \( U \) does not depend on \( \dot{x} \) then we used the second condition that \( x \) does not depend on \( \dot{q} \). Adding the two pieces together, then taking the derivative w.r.t. time,
$$
\begin{eqnarray}
\nonumber
\frac{d}{dt}\frac{\partial}{\partial\dot{q}}(T-U)&=&\sum_im\ddot{x}_i\frac{\partial x_i}{\partial q_j}
+\sum_i m\dot{x}_i\frac{\partial\dot{x}_i}{\partial q_j}.
\end{eqnarray}
$$
Now, we consider the r.h.s. of Lagrange's equations. Because the
kinetic energy depends only on \( \dot{x} \) and not \( x \), and because the
potential depends on \( x \) but not \( \dot{x} \),
$$
\begin{eqnarray}
\frac{\partial}{\partial q_j}(T-U)&=&\sum_i\frac{\partial T}{\partial\dot{x}_i}\frac{\partial\dot{x_i}}{\partial q_j}
-\sum_i\frac{\partial U}{\partial x_i}\frac{\partial x_i}{\partial q_j}\\
\nonumber
&=&\sum_i m\dot{x}_i\frac{\partial\dot{x_i}}{\partial q_j}
-\sum_i\frac{\partial U}{\partial x_i}\frac{\partial x_i}{\partial q_j}
\end{eqnarray}
$$
Using the fact that \( m\ddot{x}_i=-(\partial/\partial x_i)U \), one can
see that the bottom expressions above are identical,
$$
\begin{equation}
\frac{d}{dt}\frac{\partial}{\partial\dot{q}_i}(T-U)=\frac{\partial}{\partial q_i}(T-U).
\label{_auto99}
\end{equation}
$$
Two examples are presented here. In the first, there are two
generalized coordinates, but the two equations of motion can be
reduced to one through conservation laws (angular momentum in this
case). In the second, there is a time-dependent constraint.
Consider a cone of half angle \( \alpha \) standing on its tip at the
origin. The surface of the cone is defined as
$$
r=\sqrt{x^2+y^2}=z\tan \alpha.
$$
Find the equations of motion for a particle of mass \( m \) moving along the surface under the influence of a constant gravitational force, \( -mg\hat{z} \). For generalized coordinates use the azimuthal angle \( \phi \) and \( r \).
The kinetic energy is
$$
\begin{eqnarray*}
T&=&\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m(\dot{r}^2+\dot{z}^2)\\
&=&\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m\dot{r}^2\left(1+\cot^2\alpha\right)\\
&=&\frac{1}{2}mr^2\dot{\theta}^2+\frac{1}{2}m\dot{r}^2\csc^2\alpha.
\end{eqnarray*}
$$
The potential energy is
$$
U=mgr\cot\alpha,
$$
so Lagrange's equations give
$$
\begin{eqnarray*}
\frac{d}{dt}\left(mr^2\dot{\theta}\right)&=&0,\\
\frac{d}{dt}\left(m\csc^2\alpha \dot{r}\right)&=&mr\dot{\theta}^2-mg\cot\alpha,\\
\ddot{r}&=&r\dot{\theta}^2\sin^2\alpha-g\cos\alpha\sin\alpha
\end{eqnarray*}
$$
The first equation is a statement of the conservation of angular
momentum with \( L=mr^2\dot{\theta} \), so the second equation can also be
expressed as
$$
\ddot{r}=\frac{L^2\sin^2\alpha}{m^2r^3}-g\sin\alpha\cos\alpha.
$$
A bead slides along a wire bent in the shape of a parabola,
$$
z=\frac{1}{2}kr^2,~~r^2=x^2+y^2.
$$
Also, the parabolic wire is rotating about the \( z \) axis with angular
velocity \( \omega \). Derive the equations of motion. Are there any
stable configurations?
Using the fact that
$$
\dot{z}=\dot{r}\frac{\partial z}{\partial r}=kr\dot{r},
$$
the kinetic and potential energies are
$$
\begin{eqnarray*}
T&=&\frac{1}{2}m\left(\dot{r}^2+\dot{z}^2+r^2\omega^2\right)\\
&=&\frac{1}{2}m\left(\dot{r}^2+(kr\dot{r})^2+r^2\omega^2\right),\\
U&=&mgkr^2/2.
\end{eqnarray*}
$$
The equations of motion are then
$$
\begin{eqnarray*}
\frac{d}{dt}\left\{m\dot{r}(1+k^2r^2)\right\}&=&-mgkr+mk^2\dot{r}^2r+m\omega^2r,\\
\ddot{r}&=&\frac{-gkr+\omega^2r-k^2\dot{r}^2r}{1+k^2r^2}
\end{eqnarray*}
$$
For a stable configuration, there needs to be a solution with
\( \dot{r}=0 \) and \( \ddot{r}=0 \). This can only happen at \( r=0 \), and then
for the acceleration to be inward for small deviations of \( r \) one
needs to have \( gk>\omega^2 \). If \( \omega^2>gk \) the bead will move
outward indefinitely.
Two examples are provided for solving for normal modes. These are
solutions with multiple generalized coordinates, where the motion is
that of simple harmonic motion. However, the motion is only simple for
a particular set of coordinates \( q_1 \) and \( q_2 \),
$$
\begin{eqnarray}
q_1&=&A\cos(\omega_1 t),\\
\nonumber
q_2&=&B\cos(\omega_2 t),
\end{eqnarray}
$$
while it is not necessarily simple in other coordinates. For example
if \( x=q_1+q_2 \), and \( y=q_1-q_2 \), the \( x \) and \( y \) motions will contain
mixtures of multiple frequencies. For many problems, or in the limit
of small vibrations about a minimum, there is some coordinate system
where the motion is simple. These are normal modes. Characterizing the
normal modes involves finding the frequencies, \( \omega_i \), and the
coordinate system where the motion is simple for each coordinate. This
involves finding the direction, or the linear combination of \( x_i \)
that form the coordinates \( q_i \) in which the motion is that of a
single oscillator in each coordinate.
For a first example, we consider a system of springs, where we write
the Lagrangian, then find the normal modes. For the second example, a
double pendulum is considered. In this case, one must first make a
small angle expansion before finding the modes. In principle, problems
could have the same number of normal modes a degrees of freedom. For
example, a system of 7 particles moving in three dimensions has 21
degrees of freedom. However, some of the degrees of freedom do not
have oscillatory behavior. For example, for a rigid body in free
space, the angles describing the orientation evolve, but do not
oscillate. Also, the center-of-mass coordinates of a system of
particles isolated from outside particles moves at constant
velocity. One can also describe these as normal modes, but acknowledge
that their characteristic frequency is zero, as there are no restoring
forces.
Consider two springs, whose relaxed lengths are \( \ell \), connected to three masses as depicted in the figure here. Describe the two normal modes of the motion. We can write the Lagrangian as
$$
\begin{eqnarray*}
\mathcal{L}&=&\frac{m}{2}\dot{x}_1^2+m\dot{x}_2^2+\frac{m}{2}\dot{x}_3^2
-\frac{k}{2}(x_2-x_1-\ell)^2-\frac{k}{2}(x_3-x_2-\ell)^2.
\end{eqnarray*}
$$
There are three coordinates, thus there are three equations of motion,
$$
\begin{eqnarray*}
m\ddot{x}_1&=&-k(x_1-x_2+\ell)\\
2m\ddot{x}_2&=&-k(x_2-x_1-\ell)-k(x_2-x_3+\ell)\\
&=&-k(2x_2-x_1-x_3)\\
m\ddot{x}_3&=&-k(x_3-x_2+\ell).
\end{eqnarray*}
$$
This is a bit complicated because the center-of-mass motion does not easily separate from the three equations. Instead, choose the following coordinates,
$$
\begin{eqnarray*}
X&=&\frac{x_1+2x_2+x_3}{4},\\
q_1&=&x_1-x_2+\ell,\\
q_3&=&x_3-x_2-\ell.
\end{eqnarray*}
$$
In these coordinates the potential energy only involves two coordinates,
$$
\begin{eqnarray*}
U&=&\frac{k}{2}(q_1^2+q_3^2).
\end{eqnarray*}
$$
To express the kinetic energy express \( x_1, x_2 \) and \( x_3 \) in terms of
\( X \), \( q_1 \) and \( q_3 \),
$$
\begin{eqnarray*}
x_1&=&(3q_1-q_3-4\ell+4X)/4,\\
x_2&=&(4X-q_1-q_3)/4,\\
x_3&=&(3q_3-q_1+4\ell+4X)/4.
\end{eqnarray*}
$$
The kinetic energy and Lagrangian are them
$$
\begin{eqnarray*}
T&=&\frac{m}{2}\frac{1}{16}(3\dot{q}_1-\dot{q}_3+4\dot{X})^2
+m\frac{1}{16}(4\dot{X}-\dot{q}_1-\dot{q}_3)^2
+\frac{m}{2}\frac{1}{16}(3\dot{q}_3-\dot{q}_1+4\dot{X})^2\\
&=&\frac{3m}{8}(\dot{q}_1^2+\dot{q}_3^2)-\frac{m}{4}\dot{q}_1\dot{q}_3
+2m\dot{X}^2,\\
\mathcal{L}&=&\frac{3m}{8}(\dot{q}_1^2+\dot{q}_3^2)-\frac{m}{4}\dot{q}_1\dot{q}_3
+2m\dot{X}^2-\frac{k}{2}q_1^2-\frac{k}{2}q_3^2.
\end{eqnarray*}
$$
The three equations of motion are then,
$$
\begin{eqnarray*}
\frac{3}{4}m\ddot{q}_1-\frac{1}{4}m\ddot{q}_3&=&-kq_1,\\
\frac{3}{4}m\ddot{q}_3-\frac{1}{4}m\ddot{q}_1&=&-kq_3,\\
4M\ddot{X}&=&0.
\end{eqnarray*}
$$
The last equation simply states that the center-of-mass velocity is
fixed. One could obtain the same result by summing the equations of
motion for \( x_1 \), \( 2x_2 \) and \( x_3 \) above. The second two equations are
more complicated. To solve them, we assume a form
$$
\begin{eqnarray*}
q_1&=&Ae^{i\omega t},\\
q_3&=&Be^{i\omega t},
\end{eqnarray*}
$$
Because this is a linear equation, we can multiply the solution by a
constant and it will still be a solution. Thus, we can set \( B=1 \), then
solve for \( A \), effectively solving for \( A/B \). Putting this guess into
the equations of motion,
$$
\begin{eqnarray*}
-\frac{3}{4}\frac{A}{B}\omega^2+\frac{1}{4}\omega^2&=&-\omega_0^2\frac{A}{B},\\
-\frac{3}{4}\omega^2+\frac{1}{4}\frac{A}{B}\omega^2&=&-\omega_0^2.
\end{eqnarray*}
$$
This is two equations and two unknowns, \( \omega^2 \) and
\( A/B \). Substituting for \( A/B \) gives a quadratic equation,
$$
\begin{eqnarray*}
\omega^4-3\omega_0^2\omega^2+2\omega_0^4&=&0,\\
\omega_0^2&\equiv&k/m.
\end{eqnarray*}
$$
The two solutions are
$$
\begin{eqnarray*}
(1)~~\omega&=&\omega_0,~~~A=-B,\\
(2)~~\omega&=&\omega_0\sqrt{2},~~~A=B.
\end{eqnarray*}
$$
The first solution corresponds to the two outer masses moving in
opposite directions, in sync, with the middle mass fixed. The second
solution has both outer masses moving in the same direction, but with
the center mass moving opposite. These two solutions are referred to
as normal modes, and are characterized by their frequency and by the
linear combinations of coordinates that oscillate together. In
general, the solution is a linear combination of normal modes, which
usually results in a chaotic looking motion. However, once the
solution is expressed in terms of the normal modes, each of which
oscillates independently in a simple manner, one can better understand
the motion. Further, the frequencies of these modes represent the
natural resonant frequencies of the system. This is important in the
construction of many structures, such as bridges or vehicles.
Consider a double pendulum confined to the \( x-y \) plane, where \( y \) is
vertical. A mass \( m \) is connected to the ceiling with a massless
string of length \( \ell \). A second mass \( m \) hangs from the first mass
with an identical massless string of the same length. Using \( \theta_1 \)
and \( \theta_2 \) to describe the orientations of the strings relative to
the vertical axis, find the Lagrangian and derive the equations of
motion, both for arbitrary angles and in the small-angle
approximation. Finally, express the equations of motion in the limit
of small oscillations.
The kinetic and potential energies are:
$$
\begin{eqnarray*}
T&=&\frac{1}{2}m\ell^2\dot{\theta}_1^2
+\frac{1}{2}m\left\{(\ell\dot{\theta}_1\cos\theta_1+\ell\dot{\theta}_2\cos\theta_2)^2
+(\ell\dot{\theta}_1\sin\theta_1+\ell\dot{\theta}_2\sin\theta_2)^2\right\}\\
&=&\frac{1}{2}m\ell^2\left\{2\dot{\theta}_1^2+\dot{\theta}_2^2+2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2)
\right\},\\
U&=&mg\ell(1-\cos\theta_1)+mg\left[\ell(1-\cos\theta_1)+\ell(1-\cos\theta_2)\right]\\
&=&mg\ell(3-2\cos\theta_1-\cos\theta_2)
\end{eqnarray*}
$$
Lagrange's equations for \( \theta_1 \) lead to
$$
\begin{eqnarray*}
m\ell^2\frac{d}{dt}\left\{2\dot{\theta}_1+\dot{\theta}_2\cos(\theta_1-\theta_2)\right\}&=&
-m\ell^2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)
-2mg\ell\sin\theta_1,\\
2\ddot{\theta}_1+\ddot{\theta}_2\cos(\theta_1-\theta_2)+\dot{\theta}_2^2\sin(\theta_1-\theta_2)
&=&-2\omega_0^2\sin\theta_1,\\
\omega_0^2&\equiv& g/\ell,
\end{eqnarray*}
$$
and the equations for \( \theta_2 \) are
$$
\begin{eqnarray*}
m\ell^2\frac{d}{dt}\left\{\dot{\theta}_2+\dot{\theta}_1\cos(\theta_1-\theta_2)\right\}&=&
m\ell^2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-mg\ell\sin\theta_2,\\
\ddot{\theta}_2+\ddot{\theta_1}\cos(\theta_1-\theta_2)&=&
-\omega_0^2\sin\theta_2.
\end{eqnarray*}
$$
For small oscillations, one can only consider terms linear in \( \theta_1 \) and \( \theta_2 \) or their derivatives,
$$
\begin{eqnarray}
2\ddot{\theta}_1+\ddot{\theta}_2&=&-2\omega_0^2\theta_1,\\
\nonumber
\ddot{\theta}_1+\ddot{\theta}_2&=&-\omega_0^2\theta_2.
\end{eqnarray}
$$
To find the solutions, assume they are of the form
\( \theta_1=Ae^{i\omega t}, \theta_2=Be^{i\omega t} \). Solve for \( \omega \)
and \( A/B \), noting that \( B \) is arbitrary.
Plug in the desired form and find
$$
\begin{eqnarray*}
e^{i\omega t}(-2\omega^2A-\omega^2B)&=&e^{i\omega t}(-2\omega_0^2A),\\
e^{i\omega t}(-\omega^2A-\omega^2B)&=&e^{i\omega t}(-\omega_0^2B).
\end{eqnarray*}
$$
We can treat \( B \) as arbitrary and set it to unity. When we find \( A \),
it is the same as \( A/B \) for arbitrary \( B \). This gives the equations
$$
\begin{eqnarray*}
2\omega^2A+\omega^2&=&2\omega_0^2A,\\
\omega^2A+\omega^2&=&\omega_0^2.
\end{eqnarray*}
$$
This is two equations and two unknowns. Solving them leads to a
quadratic equation with solutions
$$
\begin{eqnarray*}
A/B&=&\pm\frac{1}{\sqrt{2}},\\
\omega^2&=&\frac{\omega_0^2}{1\pm 1/\sqrt{2}}.
\end{eqnarray*}
$$
Again, these two solutions are the normal modes, and the general
solution is a sum of the two solutions, with two arbitrary
constants. For the angles \( \theta_1 \) and \( \theta_2 \) are:
$$
\begin{eqnarray*}
\theta_1&=&\frac{A_+}{\sqrt{2}}e^{i\omega_+t}, ~\theta_2=A_+e^{i\omega_+t},\\
\theta_1&=&\frac{-A_-}{\sqrt{2}}e^{i\omega_-t}, ~\theta_2=A_-e^{i\omega_-t},\\
\omega_{\pm}&=&\omega_0\sqrt{\frac{1}{1\pm 1/\sqrt{2}}}.
\end{eqnarray*}
$$
One can also express the solution in vector notation, with the vectors
having arbitrary amplitudes \( A_+ \) and \( A_- \),
$$
\begin{eqnarray*}
\theta_+&=&\left(\begin{array}{c}
\frac{1}{\sqrt{2}}\\ 1\end{array}\right)A_+e^{i\omega_+t},\\
\theta_-&=&\left(\begin{array}{c}
\frac{-1}{\sqrt{2}}\\ 1\end{array}\right)A_-e^{i\omega_-t}.
\end{eqnarray*}
$$
Here, the upper/lower components of the vector describe
\( \theta_1/\theta_2 \) respectively.
These problems can be treated as linear algebra exercises. Linear
algebra is not used in this course, but nonetheless we describe how
this works for the curious student. In the limit of small vibrations,
the equations of motion can be expressed in the form,
$$
\begin{eqnarray*}
M\ddot{q}&=&-Kq,
\end{eqnarray*}
$$
a form that looks like the spring equation. However, \( q \) is an
$n-$dimensional vector and \( M \) and \( k \) are \( n\times n \) matrices. In
the double pendulum example, the dimensionality is 2 and the \( q \)
refers to the \( \theta_1 \) and \( \theta_2 \), and the matrices for \( M \) and
\( K \) can be read off (add ref).
$$
\begin{eqnarray*}
M&=&\left(\begin{array}{cc}
2&1\\
1&1\end{array}\right)~,\hspace*{40pt} K=\left(\begin{array}{cc}
2\omega_0^2&0\\
0&\omega_0^2\end{array}\right).
\end{eqnarray*}
$$
Multiplying both sides of the equation by the inverse matrix \( M^{-1} \),
$$
\begin{eqnarray*}
\ddot{q}&=&-\left(M^{-1}K\right)q.
\end{eqnarray*}
$$
Here,
$$
\begin{eqnarray*}
M^{-1}&=&\left(\begin{array}{cc}
1&-1\\
-1&2\end{array}\right),\\
M^{-1}K&=&\left(\begin{array}{cc}
2&-1\\
-2&2\end{array}\right)\omega_0^2.
\end{eqnarray*}
$$
One can find a transformation, basically a rotation, that transforms
to a frame where \( M^{-1}K \) is diagonal. In this coordinate system the
diagonal components of \( M^{-1}K \) represent the squared frequencies of
the normal modes,
$$
M^{-1}K\rightarrow -\left(\begin{array}{cc}
\omega_+^2&0\\
0&\omega_-^2\end{array}\right)~,
$$
and are known as "eigen" frequencies. The corresponding unit vectors,
$$
\left(\begin{array}{c}
1\\0\end{array}\right)~{\rm and}~\left(\begin{array}{c}
0\\1\end{array}\right)~{\rm in~the~new~coordinate~system},
$$
can be rotated back into the original frame, and become the solutions
for the normal modes. These are then called "eigenvectors", which
are the same as the normal modes. Finding the eigenfrequencies is
performed by realizing that the determinant of a matrix is unchanged
by the rotation between coordinate systems. Writing the equations of
motion as an eigenvalue problem,
$$
\begin{eqnarray}
\left[A-\lambda_i\mathbb{1}\right]u_i&=&0,~~~A\equiv M^{-1}K,~\lambda_i\equiv \omega^2_i.
\end{eqnarray}
$$
In the coordinate system where \( M^{-1}K \) is diagonal, and the forms
for \( u_i \) are simple this requires that in that system, the diagonal
elements of \( M^{-1}K \) are the eigenvalues, \( \omega_i^2 \). For each
\( \omega^2_i \), the determinant \( |A-\lambda_i\mathbb{1}| \) must
vanish. This is then true in any coordinate system,
$$
\begin{eqnarray}
{\rm det}\left[A-\lambda\mathbb{1}\right]&=&0,
\end{eqnarray}
$$
which for a \( 2\times 2 \) matrix becomes
$$
\begin{eqnarray}
\left|
\begin{array}{cc}
A_{11}-\lambda&A_{12}\\
A_{21}&A_{22}-\lambda
\end{array}
\right|&=&0,\\
A_{11}A_{22}-\lambda A_{11}-\lambda A_{22}+\lambda^2-A_{21}A_{12}&=&0.
\end{eqnarray}
$$
One can solve a quadratic equation for \( \lambda \), which gives two
eigenvalues corresponding to \( \omega_+^2 \) and \( \omega_-^2 \) found
above. Choosing one of the eigenvalues, one can insert one of the
eigenvalues \( \lambda_i \) into the eigenvalue problem and solve for \( u_i \),
then choose the other eigenvalue and solve for the other corresponding
vector.
If this were a 3-dimensional set of equations, the determinant would
include terms like \( \lambda^3 \) and would become a cubic equation with
three eigenvalues. One would then solve for three eigenvectors. If one
has a system with dimensionality \( n>2 \), one usually resorts to solving
the problem numerically due to the messiness of the algebra. The main
programming languages all have packages which readily diagonalize
matrices and find eigenvectors and eigenvalues.
Energy is conserved only when the Lagrangian has no explicit
dependence on time, i.e. \( L(q,\dot{q}) \), not \( L(q,\dot{q},t) \). To show
this, we first define the Hamiltonian,
$$
\begin{eqnarray}
H&=&\sum_i\left(\dot{q}_i\frac{\partial L}{\partial\dot{q}_i}\right)-L.
\end{eqnarray}
$$
After showing that \( H \) is conserved, i.e. \( (d/dt)H=0 \), we then show
that \( H \) can be identified with the total energy, \( H=T+V \).
One can see that \( H \) is conserved by applying first using the chain
rule for \( (d/dt)H \), then applying Lagrange's
equations,
$$
\begin{eqnarray}
\frac{d}{dt}H&=&\sum_i\left\{\ddot{q}_i\frac{\partial L}{\partial\dot{q}_i}+\dot{q}_i\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{q}_i}\right)-\frac{\partial L}{\partial\dot{q}_i}\ddot{q}_i-\frac{\partial L}{\partial q_i}\dot{q}_i\right\}\\
\nonumber
&=&\sum_i\left\{\ddot{q}_i\frac{\partial L}{\partial\dot{q}_i}+\dot{q}_i\frac{\partial L}{\partial q_i}-\frac{\partial L}{\partial\dot{q}_i}\ddot{q}_i-\frac{\partial L}{\partial q_i}\dot{q}_i\right\}\\
\nonumber
&=&0.
\end{eqnarray}
$$
These steps assumed that \( L \) had no explicit time dependence, i.e. \( L \)
is a function of \( q \) and \( \dot{q} \), but not of \( t \).
Next, we show that \( L \) can be identified with the energy. Because \( V \) does not depend on \( \dot{q} \),
$$
\begin{equation}
H=\sum_i\frac{\partial T}{\partial\dot{q}_i}\dot{q}_i-T+V.
\label{_auto100}
\end{equation}
$$
If the kinetic energy has a purely quadratic form in terms of \( \dot{q} \),
$$
\begin{equation}
T=\sum_{ij}A_{ij}(q)\dot{q}_i\dot{q}_j,
\label{_auto101}
\end{equation}
$$
the Hamiltonian becomes
$$
\begin{eqnarray}
H&=&\sum_{ij}2A_{ij}(q)\dot{q}_i\dot{q}_j-\sum_{ij}A_{ij}(q)\dot{q}_i\dot{q}_j+V\\
\nonumber
&=&T+V.
\end{eqnarray}
$$
The proof that \( H \) equals the energy hinged on the fact that the
kinetic energy was quadratic in \( \dot{q} \). This can be attributed to
time-reversal symmetry. Because the Cartesian coordinates \( x_i \) do not
depend on \( \dot{q}_i \) or on time, \( \dot{x}_i=(\partial x_i/\partial
q_j)\dot{q}_j \). Thus, the kinetic energy, \( T=m\dot{x}_i^2/2 \), should
be proportional to two powers of \( \dot{q} \), which validates the
assumption above.
Here, energy conservation is predicated on the Lagrangian not having
an explicit time dependence. Without an explicit time dependence the
equations of motion are unchanged if one translates a fixed amount in
time because the physics does not depend on when the clock starts. In
contrast, the absolute time becomes relevant if there is an explicit
time dependence. In fact, conservation laws can usually be associated
with symmetries. In this case the translation symmetry in time leads
to energy conservation.
For another example of how symmetry leads to conservation laws,
consider a Lagrangian for a particle of mass \( m \) moving in a
two-dimensional plane where the generalized coordinates are the radius
\( r \) and the angle \( \theta \). The kinetic energy would be
$$
\begin{equation}
T=\frac{1}{2}m\left\{\dot{r}^2+r^2\dot{\theta}^2\right\},
\label{_auto102}
\end{equation}
$$
and if the potential energy \( V(r) \) depends only on the radius \( r \) and
not on the angle, Lagrange's equations become
$$
\begin{eqnarray}
\frac{d}{dt}(m\dot{r})&=&-\frac{\partial V}{\partial r}+m\dot{\theta}^2r,\\
\nonumber
\frac{d}{dt}(mr^2\dot{\theta})&=&0.
\end{eqnarray}
$$
The second equation implies that \( mr^2\dot{\theta} \) is a
constant. Indeed, it is the angular momentum which is conserved for a
radial force. Here, the conservation of angular momentum is associated
with the independence of the physics to changes in \( \theta \), or in
other words, rotational invariance. Once one knows the fact that
\( L=mr^2\dot{\theta} \) is conserved, it can be inserted into the
equations of motion for \( \dot{r} \),
$$
\begin{equation}
m\ddot{r}=-\frac{\partial V}{\partial r}+\frac{L^2}{mr^3}.
\label{_auto103}
\end{equation}
$$
This is related to Emmy Noether's theorem
Simply stated, if the Lagrangian \( L \) is independent of \( q_i \), one can
see that the quantity \( \partial L/\partial\dot{q}_i \) is conserved,
$$
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial\dot{q}_i}=0.
\label{_auto104}
\end{equation}
$$
Another easy example is in Cartesian coordinates where the potential
depends only on \( x \) and \( y \) but not on \( z \). In that case, there is a
translational symmetry. From the last equation, this translates
into conservation of the momentum in the \( z \) direction.
Consider a pair of particles of mass \( m_1 \) and \( m_2 \) where the potential is of the form
$$
U(\boldsymbol{r}_1,\boldsymbol{r}_2)=V_a(|m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2|/(m_1+m_2))+V_b(|\boldsymbol{r}_1-\boldsymbol{r}_2|).
$$
Using symmetry arguments alone, are there any conserved components of
the momentum? or the angular momentum??
There is no translational invariance, hence there are no conserved
components of the momentum. However, there is rotational invariance
about any axis that goes through the origin. Hence, there is angular
momentum conservation in all three directions. Symmetry arguments are
great ways to recognize the existence of conserved quantities, but
actually expressing them in terms of coordinates can be tricky. For
instance, you may need to write the Lagrangian in terms of angles.
Variational Calculus and Lagrangian Formalism
Constrained motion
Example: brachiostone
Maximizing a Function
Example
Lagrangians
Proving Lagrange's Equations of Motion from Newton's Laws
Going forward with the proof, consider \( x_i(q_1,q_2\cdots,t) \) and look
at the l.h.s. of Lagrange's equations of motion.
$$
\begin{eqnarray}
\frac{\partial T}{\partial\dot{q}_j}&=&\sum_i\frac{\partial T}{\partial\dot{x}_i}\frac{\partial\dot{x}_i}{\partial\dot{q_j}}
+\sum_i\frac{\partial T}{\partial x_i}\frac{\partial x_i}{\partial\dot{q_j}}\\
\nonumber
&=&\sum_i m\dot{x}_i\frac{\partial \dot{x}_i}{\partial\dot{q_j}}\\
\nonumber
&=&\sum_i m\dot{x}_i\frac{(\delta x_i/\delta t)|_{{\rm fixed~}q_{j'\ne j}}}{\delta q_j/\delta t}\\
\nonumber
&=&\sum_im\dot{x}_i\frac{\delta{x}_i|_{{\rm fixed~}q_{j'\ne j}}}{\delta q_j}\\
\nonumber
&=&\sum_i m\dot{x}_i\frac{\partial x_i}{\partial q_j}.
\end{eqnarray}
$$
Lagrangian Examples
Small Vibrations and Normal Modes
Conservation Laws