Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part I of the course APMA0340
Introduction to Linear Algebra with Mathematica

Glossary

Preface

Many differential equations arise from problems in mechanics, electrical engineering, quantum mechanics, and other areas where conservation laws may be applied. In particular, a conservation law states that some physical quantity, which is usually energy, remains constant. In reality, a physical system is never conservative. However, mathematical models often neglect effects such as friction, electrical resistance, or temperature fluctuation if they are small enough. Therefore, we operate with idealized mathematical models that may obey conservative laws. We will see later that in many cases mathematical expressions that have no physical meaning behave conservatively. It is natural to study the effect that energy has on the solutions to the system. In particular, stability analysis can be done more efficiently upon using energy approach.

Conservative systems

When modeling a physical system that includes some form of damping — due to friction, viscosity, or dissipation — linearization will usually suffice to resolve the stability or instability of equilibria. However, when dealing with conservative systems, when damping is absent and energy is preserved, the linearization test is often inconclusive, and one must rely on more sophisticated stability criteria (see Lyapunov section). A typical example gives oscillations of ideal pendulum. In such situations, one can often exploit conservation of energy that minimizers of an energy function should be stable (but not necessarily asymptotically stable) equilibria.

By saying that energy is conserved, we mean that it remains constant as the solution evolves. Conserved quantities are also known as first integrals for the system of ordinary differential equations. Additional well-known examples include the laws of conservation of mass, and conservation of linear and angular momentum. Let us mathematically formulate the general definition.

A first integral of an autonomous system \( \dot{y} = {\bf F}({\bf y}) \) is a real-valued function I(y) that is constant on solutions of the given system.
In other words, for each solution y(t) to the differential equation \( \dot{y} = {\bf F}({\bf y}) \) , we have
\[ I({\bf y}(t)) = \mbox{constant} \qquad \mbox{for all $t$.} \]
The value of constant is fixed by the initial data since, in particular, \( c = I({\bf y}(t_0 )) = I({\bf y}_0 ) . \) Or, to rephrase this condition in another, equivalent manner, every solution to the dynamical system is constrained to move along a single level set { I(y) = c } of the first integral, namely the level set that contains the initial data y(t0).

Note first that any constant function is trivially a first integral, but this provides no useful information whatsoever about the solutions, and so is uninteresting. We will call any autonomous system that possesses a nontrivial first integral a conservative system. A physical system is conservative when it is driven by a conservative force. From mathematical point of view, we don't need to appeal to a conservative force and we can work with the corresponding first integral even if it may have no physical meaning.

Consider a mechanical system that is governed by Newton's second law,

\[ F= m\,\ddot{y} ,\qquad \ddot{y} = \text{d}^2 y/\text{d}t^2 , \]
where the force F = F(y) depends only on displacement y.

By dividing the initial equation \( F= m\,\ddot{y} \) throughout by the bothersome mass m, we can rewrite it as

\[ \ddot{y} + f(y) =0, \]
where f(y) = -F(y)/m. Now we show that this differential equation possesses a conservative law. We first multiply the equation by \( \dot{y} , \) obtaining
\[ \dot{y}\,\ddot{y} + f(y)\,\dot{y} =0. \]
Recalling the chain rule of calculus, we see that the first term on the left-hand side is
\[ \dot{y}\,\ddot{y} = \frac{\text{d}}{\text{d}t} \left[ \frac{1}{2}\, (\dot{y} )^2 \right] . \]
Likewise, if Π(y) denotes an antiderivative of f(y), we express the second term as
\[ f(y)\,\dot{y} = \frac{\text{d}}{\text{d}t} \, \Pi (y), \qquad \mbox{with} \quad \frac{\text{d}}{\text{d}y} \, \Pi (y) = f(y). \]
Using these derivative expressions, we form the differential equation
\[ \frac{\text{d}}{\text{d}t} \left[ \frac{1}{2}\, (\dot{y} )^2 + \Pi (y) \right] =0. \]
Recognizing in the expression \( (1/2)(\dot{y})^2 , \) the kinetic energy, and in Π (y), the potential energy of the system, we see that the total energy is a constant (denoted by K):
\[ E(y, \dot{y}) \equiv \frac{1}{2}\, (\dot{y} )^2 + \Pi (y) = K . \]
The above equation is the underlying conservative law. For our mechanical system, if y(t) represents a displacement, then the term \( (1/2) (\dot{y})^2 /) is kinetic energy per unit mass, and Π (y) is potential energy per unit mass.

To verify that the real-valued function I(y) is first integral, we just need to differentiate it with respect to t and invoking the chain rule leads to the basic condition

\[ 0 = \frac{\text d}{{\text d}t}\, I({\bf y}) = \nabla I({\bf y}) \cdot {\bf F} ({\bf y}) . \]
The final expression can be identified as the directional derivative of I(y) with respect to the vector field v = F(y) hat specifies the differential equation. Expanding the gradient, we obtain
\[ F_1 \left( y_1 , y_2 , \ldots , y_n \right) \frac{\partial I}{\partial y_1} + \cdots F_n \left( y_1 , y_2 , \ldots , y_n \right) \frac{\partial I}{\partial y_n} = 0 , \]
which shows that the first integral satisfies the first order differential equation. We don't have an appropriate method to solve this first order partial differential equation. However, for two dimensional case, we may try.

Let us consider the following system of autonomous equations

\[ \frac{{\text d}x}{{\text d}t} = f(x,y) , \qquad \frac{{\text d}y}{{\text d}t} = g(x,y) . \]
Any first integral I(x,y) of the above system must satisfy the linear partial differential equation
\[ f(x,y)\,\frac{\partial I}{\partial x} + g(x,y) \,\frac{\partial I}{\partial y} = 0 . \]
This first order partial differential equation can be solved by the method of characteristics. Consider the auxiliary first order scalar ordinary differential equation
\[ \frac{{\text d}y}{{\text d}x} =\frac{{\text d}y/{\text d}t}{{\text d}x/{\text d}t} = \frac{g(x,y)}{f(x,y)} . \]
Theorem: Let I(y) be a first integral for the autonomous system \( \dot{y} = {\bf F}({\bf y}) . \) If y* is a strict local extremum — minimum or maximum — of I, then y* is a stable equilibrium point for the system.

The energy functions typically do not have local maxima. Indeed, physical energy is the sum of kinetic and potential contributions. While potential energy can admit maxima, e.g., the pendulum at the top of its arc, these are only unstable saddle points for the full energy function.

Generally speaking, it is a tough problem to find a first integral for a system of differential equations. For a conservative system, typically the total energy is a first integral.

Example: Consider the equations of motions of n vortices of strength γi with positions ri = (xi,yi) in the plane

\[ \dot{x}_i = - \sum_{j\ne j} \gamma_j \,\frac{y_i - y_j}{\left\vert {\bf r}_i - {\bf r}_j \right\vert^2} , \qquad \dot{y}_i = \sum_{j\ne j} \gamma_j \, \frac{x_i - x_j}{\left\vert {\bf r}_i - {\bf r}_j \right\vert^2} \]
Upon introducing the Poisson brackets
\[ \left\{ f, g \right\} = \sum_{i=1}^n \frac{1}{\gamma_i} \left( \frac{\partial f}{\partial x_i}\, \frac{\partial g}{\partial y_i} - \frac{\partial f}{\partial y_i}\, \frac{\partial g}{\partial x_i} \right) \]
and the Hamiltonian
\[ H = - \frac{1}{2} \gamma_i \gamma_j \,\ln \left\vert {\bf r}_i - {\bf r}_j \right\vert , \]
we can rewrite the equations of motion as
\[ \dot{x}_i = \left\{ x_i , H \right\} , \qquad \dot{y}_i = \left\{ y_i , H \right\} . \]

We show that the following quantiries are conservative:

\[ P_x = \sum_{i=1}^n \gamma_i y_i \qquad\mbox{and} \qquad P_y = - \sum_{i=1}^n \gamma_i x_i . \]
Indeed, Poisson brackets give
\[ \dot{P}_x = \frac{{\text d} P_x}{{\text d}t} = \left\{ P_x , H \right\} , \qquad \dot{P}_y = \frac{{\text d} P_y}{{\text d}t}= \left\{ P_y , H \right\} , \qquad \]
Direct calculations show that
\[ \dot{P}_x = \sum_k \gamma_k \dot{y}_k = \sum_{i\ne k} \gamma_i \gamma_k \, \frac{x_k - x_i}{\left\vert {\bf r}_k - {\bf r}_i \right\vert^2} = 0. \]
It is easy to see that the result is zero since the summand is antisymmetric in i and k. Similarly, one can prve that \( \dot{P}_y = 0 , \) too. Thus Px and Py are integrals of motion. As argued above,
\[ \left\{ P_x , H \right\} = \left\{ P_y , H \right\} = 0 \]
and
\[ \left\{ P_x , P_y \right\} = \sum_k \frac{1}{\gamma_k} \left( \frac{\partial P_x}{\partial x_k} \,\frac{\partial P_y}{\partial y_k} - \frac{\partial P_x}{\partial y_k} \,\frac{\partial P_y}{\partial x_k} \right) = \sum_k \gamma_k . \]

For the system of two vortices the equations of motion read

\[ \begin{split} \dot{x}_1 &= - \gamma_2 \,\frac{y_1 - y_2}{\left\vert {\bf r}_1 - {\bf r}_2 \right\vert^2} , \qquad \dot{x}_2 &= - \gamma_1 \, \frac{y_2 - y_1}{\left\vert {\bf r}_2 - {\bf r}_1 \right\vert^2} , \\ \dot{y}_1 &= \gamma_2 \,\frac{x_1 - x_2}{\left\vert {\bf r}_1 - {\bf r}_2 \right\vert^2} , \qquad \dot{y}_2 &= \gamma_1 \, \frac{x_2 - x_1}{\left\vert {\bf r}_2 - {\bf r}_1 \right\vert^2} . \end{split} \]
The integrals of motion or two vortices take the form
\[ P_x = \gamma_1 y_1 + \gamma_2 y_2 , \qquad P_y = -\gamma_1 x_1 - \gamma_2 x_2 . \]
They can be interpreted as coordinates of the center of mass of the two vortices. The relative motion of the vortices is described by the variables
\[ X = x_1 - x_2 , \qquad Y= y_1 - y_2 . \]
Ror X and Y one obtains the equations of motion
\[ \dot{X} = - \left( \gamma_1 + \gamma_2 \right) \frac{Y}{R^2} , \qquad \dot{Y} = \left( \gamma_1 + \gamma_2 \right) \frac{X}{R^2} \]
with
\[ R^2 = X^2 + Y^2 . \]
Note that the distance between the vortices R is a constant of motion:
\[ \frac{{\text d}R^2}{{\text d}t} = 2X\dot{X} + 3Y\dot{Y} = -2 \left( \gamma_1 + \gamma_2 \right) \frac{XY}{R^2} + 2 \left( \gamma_1 + \gamma_2 \right) \frac{XY}{R^2} =0 . \]
The equations of motion can be rewritten as
\[ \dot{X} = -\omega_0 Y, \qquad \dot{Y} = \omega_0 X, \]
where
\[ \omega_0 = \frac{\gamma_1 + \gamma_2}{R^2} \]
is an angular velocity. Introducing a complex variable
\[ Z = X + {\bf j}Y , \]
the equation of motion can be rewritten as
\[ \dot{Z} = {\bf j}\omega_0 Z . \]
   ■

Example: The total energy

\[ E = \mbox{K} + \Pi = \frac{1}{2}\, m\,\ell^2 \dot{\theta}^2 + mg\,\ell \left( 1 - \cos \theta \right) \]
for the pendulum equation
\[ \ddot{\theta} + \frac{g}{\ell}\,\sin \theta =0 \]
is the first integral. The Lagrangian is the difference of kinetic and potential energies:
\[ {\cal L} = \mbox{K} - \Pi = \frac{1}{2}\, m\,\ell^2 \dot{\theta}^2 - mg\,\ell \left( 1 - \cos \theta \right) . \]
The conjugate momentum is
\[ p_{\theta} = \frac{\partial {\cal L}}{\partial \dot{\theta}} = m\,\ell^2 \dot{\theta} , \]
where ℓ is the length of pendulum rod and g is the acceleration due to gravity.    ■

Example: Consider the motion of a rigid rod rocking back and forth on the circular surface without slipping. The governing equation of motion can be expressed as (see the article by Wu et al.)

\[ \left( \frac{1}{12} + \frac{u^2}{16} \right) \frac{{\text d}^2 u}{{\text d} t^2} + \frac{u}{16} \left( \frac{{\text d}u}{{\text d}t} \right)^2 + \frac{g\,u}{4\,\ell}\,\cos u =0 \]
subject to the initial conditions
\[ u(0) = A, \qquad \dot{u}(0) =0 . \]
Here \( \dot{u} = \texttt{D} u = {\text d}u/{\text d}t \) is the derivative of unknown function u(t) and g, ℓ are positive constants. It is also convenient to assume that A in the initial conditions is a positive constant (considered as amplitude, in this content).

Its first integral is obtained by multiplying the governing equation by the velocity \( \dot{u} \) and integration. This yields the first integral

\[ I(u) = \int_0^t \left( - \frac{1}{2}\, \dot{u}^2 - \frac{3}{8}\, u^2 \dot{u}^2 + \frac{3g}{\ell} \left( \cos u + \sin u \right) \right) {\text d}t . \]
   
StreamPlot[{y, -(x/16*y^2 + (x/4)*Cos[x])/(1/12 + x/16)}, {x, -3,
4}, {y, -3, 3}, Mesh -> 12, AspectRatio -> 1,
MeshStyle -> {Red, Thick}, PerformanceGoal -> "Quality",
StreamPoints -> Fine]
   ■

 

Hamiltonian systems

Now we turn our attention to an important particular class of conservative systems that was invented in 1835 by the Irish mathematician William Rowan Hamilton (1805--1865). Although finding a first integral for a given system of differential equations is a very tough problem in general, Hamiltonian functions provide the first integral right away. The most important features of the Hamiltonian is its ability to reformulate not only classical mechanics (with finite degrees of freedom), but also quantum mechanics.

A function of 2n variables q = [ q1, ... , qn ] and p = [ p1, ... , pn ] is called the Hamiltonian or Hamiltonian function if the motion of the system is described by
\[ \dot{q}_i = \frac{\partial H}{\partial p_i} , \qquad\dot{p}_i = - \frac{\partial H}{\partial q_i} , \qquad i =1,2,\ldots n . \]

As we see from the system of differential equations, it is written in variables that have no direct relationship with the positions of a particle system. Therefore, these variables, usually denoted historically by q, are called the generalized coordinates. Another variables are referred to as the generalized momentum, and usually denoted by p. These variables are related to the Lagrangian via the Legendre transformation

\[ {p}_i = \frac{\partial {\cal L}}{\partial \dot{q}_i} , \qquad \]
where ℒ is the Lagrangian of the system (the difference of the kinetic energy and the potential energy. The solution of Hamilton's equations of motion will yield a trajectory in terms of positions and momenta as functions of time. Again, Hamilton's equations can be easily shown to be equivalent to Newton's equations, and, like the Lagrangian formulation, Hamilton's equations can be used to determine the equations of motion of a system in any set of coordinates.
\[ H({\bf q}, {\bf p}) = \sum_{i=1}^n p_i \cdot \dot{q}_i - {\cal L} ({\bf q}, \dot{\bf q}) . \]
If the Hamiltonian H does not explicitly depend on time t, then
\[ H \quad \mbox{is a}\quad \begin{cases} \mbox{ constant of the motion}, \\ \mbox{ conserved quantity}, \\ \mbox{ first integral of the motion}. \end{cases} \]

Let us consider the simplest case---a dynamics of a particle of mass m that is described by the Hamiltonian

\[ H(x,p) = \frac{p^2}{2m} + \Pi (x) , \]
where \( \displaystyle p = m\dot{x} \) is the particle's momentum, x is a particle position (that is usually denoted by q), and Π(x) is the potential energy. The Hamilton's equations become
\[ \frac{{\text d}x}{{\text d}t} = \frac{p}{m}, \qquad \frac{{\text d}p}{{\text d}t} = - \frac{{\text d}\Pi (x)}{{\text d}x} . \]
Since the Hamiltonian (and Lagrangian) is time independent, the energy conservation law states that H(x,p) = E. In turn,this conservation law implies that the particle's velocity \( \displaystyle \dot{x} \) can be expressed as
\[ \dot{x} (x,E) = \pm \sqrt{\frac{2}{m} \left[ E - \Pi (x) \right]} , \]
where the sign of \( \displaystyle \dot{x} \) is determined from the initial conditions. It is immediately clear that physical motion is possible only if E ≥ Π(x); points where E = Π(x) are known as turning points. Furthermore, motion is bounded for a ̄fixed value of energy E if two roots exist (labeled by x1 and x2).

The dynamical solution \( \displaystyle x(t,E) \) of the Hamilton's equations is ̄first expressed as

\[ t(x,E) = \sqrt{\frac{m}{2}} \int_{x_0}^x \frac{{\text d}s}{\sqrt{E - \Pi (s)}} , \]
where x0 is the particle's initial position Next, inversion of the above relation yields the solution x(t,E). For bounded motion in one dimension, the particle bounces back and forth between the two turning points x1 and x2. The period of oscillation T(E) is a function of energy alone
\[ T(E) = 2\int_{x_1}^{x_2} \frac{{\text d}x}{\dot{x}(x,E)} = \sqrt{2m} \int_{x_1}^{x_2} \frac{{\text d}x}{\sqrt{E - \Pi (x)}} . \]

Example: The equation of motion for a charged particle of mass m and charge e moving in an electromagnetic field represented by the electric field E and magnetic field B are

\[ \begin{split} \frac{{\text d}{\bf x}}{{\text d}t} &= {\bf v} , \\ \frac{{\text d}{\bf v}}{{\text d}t} &= \frac{c}{m} \left( {\bf E} + \frac{{\text d}{\bf x}}{{\text d}t} \times {\frac{\bf B}{c} \right) , \end{split} \]
where x denotes the position of the particle and v is its velocity.

By treating the coordinates (x;v) as generalized coordinates (i.e., δv is independent of δx), we now show that the corresponding equations of motion can be obtained as Euler--Lagrange equations from the Lagrangian

\[ L({\bf x}, \dot{\bf x}, {\bf v}, \dot{\bf v}, t) = \left(m\,{\bf v} + \frac{e}{c} \,{\bf A}({\bf x}, t) \right) \cdot \dot{\bf x} - \left( e\,\Phi ({\bf x},t) + \frac{m}{2} \left\vert {\bf v} \right\vert^2 \right) , \]
where Φ and A are the electromagnetic potentials in terms of which electric and magnetic fields are defined
\[ {\bf E} = - \nabla \Phi - c^{-1} \frac{\partial {\bf A}}{\partial t} \qquad \mbox{and} \qquad {\bf B} = \nabla \times {\bf A} . \]

First, we look at the Euler--Lagrange equation for x:

\[ \begin{split} \frac{\partial L}{\partial \dot{\bf x}} &= m\,{\bf v} + \frac{e}{c} \,{\bf A} \qquad\Longrightarrow \qquad \frac{\text d}{{\text d}t} \left( \frac{\partial L}{\partial \dot{\bf x}} \right) = m\,\dot{\bf v} + \frac{e}{c} \left( \frac{\partial {\bf A}}{\partial t} + \dot{\bf x} \cdot \nabla {\bf A} \right) , \\ \frac{\partial L}{\partial {\bf x}} &= \frac{e}{c} \,\nabla {\bf A} \cdot \dot{\bf x} - e\,\nabla \Phi , \end{split \]
because
\[ m\,\dot{\bf v} = -e \left( \nabla \Phi + c^{-1} \frac{\partial {\bf A}}{\partial t} \right) + \frac{e}{c} \, \dot{\bf x} \times \nabla \times {\bf A} = e\,{\bf E} + \frac{e}{c} \,\dot{\bf x} \times {\bf B} , \]
where we used representation of fields through potentials.

Next, we look at the Euler-Lagrange equation for v:

\[ \frac{\partial L}{\partial \dot{\bf v}} = 0 \qquad\Longrightarrow \qquad \frac{\text d}{{\text d}t} \left( \frac{\partial L}{\partial \dot{\bf v}} \right) = 0 = \frac{\partial L}{\partial {\bf v}} = m\,\dot{\bf x} - m\,{\bf v} . \]
The Lagrangian is
\[ L \left( {\bf x}, \dot{\bf x}, t \right) = \frac{m}{2} \left\vert \dot{\bf x} \right\vert^2 + \frac{e}{c} \, {\bf A}({\bf x},t) \cdot \dot{\bf x} -e\,\Phi ({\bf x},t) . \]
The energy conservation law is expressed as
\[ \frac{\partial L}{\partial t} = \frac{\text d}{{\text d} t} \left( L - \dot{\bf x} \cdot \frac{\partial L}{\partial \dot{\bf x}} -\dot{\bf v} \cdot \frac{\partial L}{\partial \dot{\bf v}} \right) , \]
from which we derive the equations of motion.

The canonical momentum p for a particle of mass m and charge e in an electric field defined as

\[ {\bf p} ({\bf x}, {\bf v}, t) = \frac{\partial L}{\partial \dot{\bf x}} = m\,{\bf v} + \frac{e}{c} \, {\bf A}({\bf x},t) . \]
The canonical Hamiltonian function H(x,p,t) is now constructed through the Legendre transformation
\begin{align*} H ({\bf x}, {\bf p}, t) &= {\bf p}\cdot \dot{\bf x}({\bf x}, {\bf p}, t) - L\left( {\bf x}, \dot{\bf x}({\bf x}, {\bf p}, t), t \right) \\ &= e\,\Phi ({\bf x},t) + \frac{1}{2m} \left\vert {\bf p} - \frac{e}{c} \, {\bf A}({\bf x},t) \right\vert^2 , \end{align*}
where v(x,p,t) was obtained by inverting p(x,v,t). Using the canonical Hamiltonian function, we get
\begin{align*} \dot{\bf x} &= \frac{\partial H}{\partial {\bf p}} = \frac{1}{m} \left( {\bf p} - \frac{e}{c} \, {\bf A}({\bf x},t) \right) , \\ \dot{\bf p} &= -\frac{\partial H}{\partial {\bf x}} = -e\,\nabla \Phi - \frac{e}{c} \, \nabla {\bf A} \cdot \dot{\bf x} , \end{align*}
from which we recover the equations of motion.    ■

Example: A spherical pendulum of length ℓ and mass m carries a positive charge e and moves under the action of a constant gravitational field (with acceleration g) and a constant magnetic ̄field B. The position vector of the bob of the pendulum is

\[ {\bf x} = \ell \left[ \cos\theta \left( \cos\phi \,\hat{x} + \sin \phi\, \hat{y} \right) + \sin \theta\,\hat{z} \right] \]
and thus its velocity \( {\bf v} = \dot{\bf x} \) is
\[ {\bf v} = \ell \dot{\theta} \left[ \sin\theta \left( \cos\phi \,\hat{x} + \sin \phi\, \hat{y} \right) - \cos \theta\,\hat{z} \right] + \ell \sin\theta \,\dot{\phi} \left( -\sin\phi \,\hat{x} + \cos \phi\, \hat{y} \right) , \]
and the kinetic energy of the pendulum becomes
\[ \mbox{K} = \frac{m\,\ell^2}{2} \left( \dot{\theta}^2 + \sin^2 \theta \,\dot{\phi}^2 \right) . \]
Because the charged pendulum moves in a magnetic ̄field \( {\bf B} = -B\,\hat{z} , \) we must include the magnetic term \( {\bf v} \cdot e\,{\bf A}/c \) in the Lagrangian. Here, the vector potential A must be evaluated at the position of the pendulum and is thus expressed as
\[ {\bf A} = - \frac{B}{2}\,\ell \,\sin\theta \left( - \sin^2 \phi\,\hat{x} + \cos\phi\,\hat{y} \right) , \]
and we find
\[ \frac{e}{c} \, {\bf A}\dot {\bf v} = - \frac{B}{2}\,\ell \,\sin\theta \, \dot{\phi} . \]
Finally, the charged pendulum is under the influence of two potential energy terms: gravitational potential energy \( \left( -mg\,\ell\,\cos\theta \right) \) and magnetic potential energy \( \left( -{\bf \mu} \cdot {\bf B} = \mu_z B \right) \) where
\[ {\bf \mu} = \frac{e}{m\,c} \left({\bf x} \times {\bf v} \right) \]
denotes the magnetic moment of a charge e moving about a magnetic field line. Here, it is easy to ̄find
\[ \mu_z B = - \frac{e\,B}{2c} \,\ell^2 \sin^2 \theta \,\dot{\phi} \]
and by combining the various terms, the Lagrangian for the system becomes
\[ L \left( \theta , \dot{\theta} , \dot{\phi} \right) = m\ell^2 \left[ \frac{1}{2}\,\dot{\theta}^2 + \sin^2 \theta \left( \frac{1}{2}\,\dot{\phi}^2 - \omega_c \dot{\phi} \right) \right] + mg\ell \,\cos \theta , \]
where the cyclotron frequency ωc is defined as
\[ \omega_c = \frac{e\,N}{mc} . \]
The Euler-Lagrange equation for θ is
\begin{align*} \frac{\partial L}{\partial \dot{\theta}} &= m\ell^2 \dot{\theta} \qquad \Longrightarrow \qquad \frac{\text d}{{\text d} t} \left( \frac{\partial L}{\partial \dot{\theta}} \right) = m\ell^2 \ddot{\theta} , \\ \frac{\partial L}{\partial \theta} &= - mg\ell \,\sin \theta + m\ell^2 \left( \dot{\phi}^2 - 2\,\omega_c \dot{\phi} \right) \sin\theta \,\cos \theta , \end{align*}
or
\[ \ddot{\theta} + \frac{g}{\ell}\,\sin\theta = \left( \dot{\phi}^2 - 2\,\omega_c \dot{\pi} \right) \sin\theta \,\cos\theta . \]
The Euler-Lagrange equation for φ immediately leads to a constant of the motion for the system since the Lagrangian is independent of the azimuthal angle φ and hence
\[ p_{\phi} = \frac{\partial L}{\partial \dot{\phi}} = m\ell^2 \sin^2 \theta \left( \dot{\phi} - \omega_c \right) \]
is a constant of the motion, i.e., the Euler-Lagrange equation for φ states that \( \dot{p}_{\phi} = 0 . \)

Since \( p_{\phi} \) is a constant of the motion, we can use it the rewrite \( \dot{\phi} \) in the Euler-Lagrange equation for θ as

\[ \dot{\phi} = \omega_c + \frac{p_{\phi}}{m\ell^2 \sin^2 \theta} , \]
and thus
\[ \dot{\phi}^2 - 2\,\omega_c \dot{\phi} = \left( \frac{p_{\phi}}{m\ell^2 \sin^2 \theta} \right)^2 - \omega_c^2 , \]
so that
\[ \ddot{\theta} + \frac{g}{\ell}\,\sin\theta = \sin\theta \,\cos\theta \left[ \left( \frac{p_{\phi}}{m\ell^2 \sin^2 \theta} \right)^2 - \omega_c^2 \right] . \]

The Hamiltonian for the system is obtained through the Legendre transformation

\begin{align*} H &= \dot{\theta}\,\frac{\partial L}{\partial \dot{\theta}} + \dot{\phi} \, \frac{\partial L}{\partial \dot{\phi}} -L , \\ &= \frac{p_{\theta}^2}{2m\,\ell^2} + \frac{1}{2m\,\ell^2 \sin^2 \theta} \left( p_{\phi} + m\ell^2 \omega_c \,\sin^2 \theta \right)^2 - mg\ell\,\cos\theta . \end{align*}
The Hamilton's equations for variables (θ,pθ) are
\begin{align*} \dot{\theta} &= \frac{\partial H}{\partial p_{\theta}} = \frac{p_{\theta}}{m\ell^2} , \\ \dot{p}_{\theta} &= - \frac{\partial H}{\partial \theta} = -mg\ell \, \sin\theta + m\ell^2 \sin\theta \,\cos \theta \left[ \left( \frac{p_{\phi}}{m\ell^2 \sin^2 \theta} \right)^2 - \omega_c^2 \right] , \end{align*}
while the Hamilton's equations for variables (φ,pφ) are
\begin{align*} \dot{\phi} &= \left\{ \phi , H \right\} = \frac{\partial H}{\partial p_{\phi}} = \frac{p_{\phi}}{m\ell^2 \sin^2 \theta} + \omega_c^2 , \\ \dot{p}_{\phi} &= \left\{ p_{\phi} , H \right\} = - \frac{\partial H}{\partial \phi} = 0 . \end{align*}
It is readily checked that these Hamilton equations lead to the same equations as the Euler--Lagrange equations for variables θ and φ.

It turns out that a most useful application of the Hamiltonian formalism resides in the use of the constants of the motion to plot Hamiltonian orbits in phase space. Indeed, for the problem considered here, a Hamiltonian orbit is expressed in the form pθ(θ, E, pφ), that is each orbit is labeled by values of the two constants of motion E (the total energy) and pφ the azimuthal canonical momentum (actually an angular momentum):

\[ p_{\theta} = \pm \sqrt{2\,m\ell^2 \left( E + mg\ell \,\cos\theta \right) - \frac{1}{\sin^2 \theta} \left( p_{\phi} + m\ell^2 \sin^2 \theta \,\omega_c \right)^2} . \]
Hence, for charged pendulum of given mass m and charge e with a given cyclotron frequency ωc (and g), we can completely determine the motion of the system once initial conditions are known (from which E and pφ can be calculated).    ■

Another subclass of conservative Hamiltonian systems was first introduced in the work of Siméon Dennis Poisson (1781--1840) in the early nineteenth century. His research was concerned with methods for solving the conservative systems arising in classical mechanics. The Vlasov--Poisson system is an important non-linear transport equation, used to describe the evolution of particles under their self-consistent electric or gravitational field. The existence of classical solutions is limited to dimensions ≤3 under strong assumptions on the initial data, while weak solutions are known to exist under milder conditions.

A Poisson system is a first order system of ordinary differential equations of the form
\[ \dot{\bf y} = {\bf J} \,\nabla H({\bf y}) , \qquad \mbox{where} \quad {\bf J}^{\mathrm T} = -{\bf J} \]
is a constant skew-symmetric matrix. The real-valued function H(y) is known as the Hamiltonian function for the system.
The Hamiltonian function is automatically a first integral for the Poisson system because
\[ \frac{{\text d}H}{{\text d}t} = \nabla H \cdot \frac{{\text d} {\bf y}}{{\text d}t} = \nabla H \cdot {\bf F} = \nabla H^{\mathrm T} {\bf J}\,\nabla H = 0 . \]
The simplest example is when \( {\bf J} = \begin{pmatrix} 0&1 \\ -1&0 \end{pmatrix} \) is a 2×2 matrix. In this case, we set y = [ q(t), p(t) ]T, and the Poisson system takes the form
\[ \frac{{\text d}q}{{\text d}t} = \frac{\partial H}{\partial p} , \qquad \frac{{\text d}p}{{\text d}t} = -\frac{\partial H}{\partial q} . \]
The most important case, generalizing the two-dimensional version is
\[ {\bf J} = \begin{pmatrix} {\bf O}& {\bf I} \\ - {\bf I} & {\bf O} \end{pmatrix} , \]
where O denotes the n×n zero matrix and I stands for the n×n identity matrix. In this case the variables are split into two vector components:
\[ {\bf y} = \left[ {\bf q}, {\bf p} \right]^{\mathrm T} = \left[ q_1 , q_2 , \ldots , q_n , p_1 , p_2 , \ldots , p_n \right]^{\mathrm T} . \]
We have in this case that
\[ \nabla H \left( {\bf q}, {\bf p} \right) = \left( \frac{\partial H}{\partial q_1} , \ldots , \frac{\partial H}{\partial q_n} , \frac{\partial H}{\partial p_1} , \ldots , \frac{\partial H}{\partial p_n} \right)^{\mathrm T} , \]
and equations of motion become
\[ \frac{{\text d}q_1}{{\text d} t} = \frac{\partial H}{\partial p_1} , \ldots , \frac{{\text d}q_n}{{\text d} t} = \frac{\partial H}{\partial p_n} , \frac{{\text d}p_1}{{\text d} t} = -\frac{\partial H}{\partial q_1} , \ldots , \frac{{\text d}p_n}{{\text d} t} = -\frac{\partial H}{\partial q_n} . \]
In physical applications, the p variables typically represent momenta, while the q variables represent positions of the objects in the system; this will be shown in the following pendulum example.

Example: Consider the undamped pendulum equation

\[ \ddot{\theta} + \omega_0^2 \sin\theta =0 , \qquad \mbox{where} \quad \omega_0^2 = g/\ell . \]
Let us introduce the position variable q = θ representing the angular coordinate of the pendulum. Let \( p = m\,\dot{\theta} \) be its angular momentum. The pendulum energy function can be rewritten in terms of the position and momentum variables:
\[ H(p, q) = \frac{p^2}{2m} + \omega_0^2 \left( 1− \cos q\right) . \]
The resulting Hamiltonian system becomes
\[ \frac{{\text d}p}{{\text d}t} = \omega_0^2 \sin q . \]
If we convert back to y = q and velocity \( v = \dot{\theta} = p/m , \) we immediately recover the phase plane form
\[ \frac{{\text d}y}{{\text d}t} = v, \qquad \frac{{\text d}v}{{\text d}t} = - \omega_0^2 \sin y . \]
of the pendulum equation. Thus, we recover the constancy of the energy first integral directly from the general Hamiltonian framework.    ■

Conservative mechanical systems (with n particles) can always be represented in Hamiltonian form

\[ H\left( {\bf p}, {\bf q} \right) = \frac{p_1^2}{2\,m_1} + \cdots + \frac{p_n^2}{2\, m_n} + \Pi \left( q_1 , \ldots , q_n \right) , \]
where each mi > 0 is a positive constant, and Π is the potential energy. The canonical system has the form
\[ \frac{{\text d}q_1}{{\text d}t} = \frac{p_1}{m_1} , \ldots , \frac{{\text d}q_n}{{\text d}t} = \frac{p_n}{m_n} , \frac{{\text d}p_1}{{\text d}t} = - \frac{\partial \Pi}{\partial q_1} , \ldots , \frac{{\text d}p_n}{{\text d}t} = - \frac{\partial \Pi}{\partial q_n} , \]
which we can rewrite in vectorial form
\[ {\bf M} \,\frac{{\text d} {\bf q}}{{\text d} t} = {\bf p} , \qquad \frac{{\text d} {\bf p}}{{\text d} t} = - \nabla \Pi ({\bf q}) , \]
where M is the diagonal matrix: M = diag{ m1, ... , mn }. Eliminating the momentum variable p, we obtain the Newtonian form
\[ {\bf M} \,\frac{{\text d}^2 {\bf q}}{{\text d} t^2} = - \nabla \Pi ({\bf q}) . \]
The vector q represents the position vector for the mechanical system, while Π(q) is the potential energy. The constants mi are masses, and each \( p_i = m\,\dot{q}_i \) represents the momentum variable associated with the position variable qi. The Hamiltonian function is the total energy, the first term being the kinetic energy because the term
\[ \frac{p_i^2}{2\,m_i} = \frac{m_i}{2} \left( \frac{{\text d}q_i}{{\text d}t} \right)^2 \]
is precisely one half mass times velocity squared.

Example: For computing the motion of two bodies (planet and sun) which attract each other, we choose one of the bodies (sun) as the center of our coordinate system; the motion will then stay in a plane and we can use two-dimensional coordinates q = (q1,q2) for the position of the second body. Newton’s laws, with a suitable normalization, then yield the following differential equations

\[ \ddot{q}_1 = - \frac{q_1}{\left( q_1^2 + q_2^2 \right)^{3/2}} , \qquad \ddot{q}_2 = - \frac{q_2}{\left( q_1^2 + q_2^2 \right)^{3/2}} , \qquad \]
This is equivalent to a Hamiltonian system with the Hamiltonian
\[ H(q_1 , q_2 , p_1 , p_2 ) = \frac{1}{2} \left( p_1^2 + p_2^2 \right) - \left( q_1^2 + q_2^2 \right)^{-1/2} . \]
The Lagrangian in polar coordinates becomes
\[ L = \frac{1}{2}\,m \left( \dot{r}^2 + r^2 \dot{\theta}^2 \right) + \frac{GM}{r} \]
The coordinate θ is cyclic and therefore
\[ \ell = \frac{\partial L}{\partial \dot{\theta}} = mr^2 \dot{\theta} \]
is conserved, but this quantity does not explicitly occur in the Lagrangian. However, the Hamiltonian in polar coordinates is
\[ H = \frac{1}{2m} \left( p_r^2 + \frac{p_{\theta}^2}{r^2} \right) - \frac{GM}{r} , \]
where the gravitational constant G = 6.674×10−11 m³3⋅kg−1⋅s−2. Hence pθ = ℓ is constant and we immediately have
\[ H = \frac{1}{2m} \left( p_r^2 + \frac{\ell^2}{r^2} \right) - \frac{GM}{r} , \]
The planet moves in elliptic orbits with the sun at one of the foci (Kepler’s first law). With initial values
\[ q_1 (0) = 1 e , \quad q_2 (0) =0, \quad \dot{q}_1 (0) =0 , \quad \dot{q}_2 (0) = \sqrt{\frac{1+e}{1-e}} , \]
the solution is an ellipse with eccentricity e (0 ≤ e ≤ 1), a =1, \( b = \sqrt{1 - e^2} , \) d = 1 - e², and the period 2π. The total energy is H0 = -1/2, and the angular momentum is L0 = b.    ■

Example: While at Princeton in 1962, the French mathematician and astronomer Michel Hénon and the American astrophysicist Carl Heiles worked on the non-linear motion of a star around a galactic center with the motion restricted to a plane. In 1964 they published an article titled "The applicability of the third integral of motion: Some numerical experiments" in The Astronomical Journal, doi: 10.1086/109234. Their original idea was to find a third integral of motion in a galactic dynamics. For that purpose they took a simplified two-dimensional nonlinear axi-symmetric potential (known now as the Hénon–Heiles potential):

\[ V(x,y) = \frac{1}{2} \left( x^2 + y^2 \right) + \lambda \left( x^2 y - \frac{y^3}{3} \right) . \]
This leads to the polynomial time‐independent Hamiltonian in two degrees of freedom
\[ H({\bf q},{\bf p}) = \frac{1}{2} \left( p_1^2 + p_2^2 \right) + \frac{1}{2} \left( q_1^2 + q_2^2 \right) + q_1^2 q_2 - \frac{1}{3}\,q_2^3 \]
generates a Hamiltonian differential equations
\[ \begin{split} \dot{x} &= p_x , \\ \dot{p}_x &= -x - 2\lambda xy , \\ \dot{y} &= p_y , \\ \dot{p}_y &= -y - \lambda \left( x^2 - y^2 \right) , \end{split} \]
that can have chaotic solutions. Their numerical studies showed that both chaotic and nonchaotic motion existed for this system. For H < 0.11 (approximately),the motion appeared to be non chaotic. For H > 0.11, the motion was in general chaotic, although for certain initial conditions some non chaotic motion still existed at the higher energies.

Upon plotting their numerical outputs, Hénon and Heiles discovered that there exists another conserved quantity in addition to the total energy H. This nontrivial conserved quantity is a first integral, which is valid over some domain of parameter space, or a configurational invariant, which is defined for some specific value of a parameter.    ■

Example: Let's use the Hamiltonian method to find the equations of motion of a particle of mass m constrained to move on the surface of a cylinder defined by x² + y² = R². The particle is subject to a force directed toward the origin and proportional to the distance of the particle from the origin: F = - kr.

Introducing cylindrical coordinates (r,θ,z), we express the potential energy as

\[ \Pi = \frac{1}{2}\,k\,r^2 = \frac{1}{2}\,k \left( x^2 + y^2 + z^2 \right) = \frac{1}{2}\,k \left( R^2 + z^2 \right) . \]
The linear velocity in cylindrical coordinates is
\[ v^2 = \dot{r}^2 + r^2 \dot{\theta}^2 + \dot{z}^2 . \]
Since the radial variable is constant, we have
\[ v^2 = R^2 \dot{\theta}^2 + \dot{z}^2 . \]
The kinetic energy becomes
\[ \mbox{K} = \frac{m}{2} \left( v^2 \right) = \frac{m}{2} \left( R^2 \dot{\theta}^2 + \dot{z}^2 \right) . \]
The Lagrangian is
\[ {\cal L} = \mbox{K} - \Pi = \frac{m}{2} \left( R^2 \dot{\theta}^2 + \dot{z}^2 \right) - \frac{1}{2}\,k \left( R^2 + z^2 \right) . \]
Our generalized coordinates are θ, z, and the generalized momenta are
\[ \begin{split} p_{\theta} &= \frac{\partial {\cal L}}{\partial \dot{\theta}} = m\,R^2 \dot{\theta} , \\ p_{z} &= \frac{\partial {\cal L}}{\partial \dot{z}} = m\,\dot{z} . \end{split} \]
Because the system is conservative and the equations of transformation between rectangular and cylindrical coordinates do not explicitly involve the time, the Hamiltonian H is just the total energy expressed in terms of the variables θ, pθ, pz:
\[ H\left( z, p_{\theta}, p_z \right) = \mbox{K} + \Pi = \frac{m}{2} \left( R^2 \dot{\theta}^2 + \dot{z}^2 \right) + \frac{1}{2}\,k \left( R^2 + z^2 \right) = \frac{p_{\theta}^2}{2m\,R^2} + \frac{p_z^2}{2m} + \frac{1}{2}\, k\,z^2 , \]
where we have noticed that there is no dependence on θ, and we suppress the constant term k R² /2. The canonical equations of motion are therefore:
\[ \begin{split} \dot{p}_{\theta} &= - \frac{\partial H}{\partial \theta} = 0 , \\ \dot{p}_z &= \frac{\partial H}{\partial z} = -k\,z , \\ \dot{theta} &= \frac{\partial H}{\partial p_{\theta}} = \frac{p_{\theta}}{m\,R^2} , \\ \dot{z} &= \frac{\partial H}{\partial p_{z}} = \frac{p_z}{m} . \end{split} \]
From the Legendre transformation, we have
\[ \begin{split} p_{\theta} &= \frac{\partial {\cal L}}{\partial \dot{\theta}} = \frac{\partial}{\partial \dot{\theta}} \,\frac{m}{2} \left( R^2 \dot{\theta}^2 + \dot{z}^2 \right) = m\,R^2 \dot{\theta} , \\ p_z &= \frac{\partial {\cal L}}{\partial \dot{z}} = \frac{\partial}{\partial \dot{z}} \,\frac{m}{2} \left( R^2 \dot{\theta}^2 + \dot{z}^2 \right) = m\,\dot{z} . \end{split} \]
Hence,
\[ \dot{\theta} = \frac{p_{\theta}}{m\,R^2} \qquad\mbox{and}\qquad \dot{z} = \frac{p_z}{m} . \]
We see immediately that pθ is a constant--conservation of angular momentum about the z (symmetry)axis. We can get the equation of motion for z from \( \dot{p}_z = -k\,z = m\,\ddot{z} , \) so
\[ \ddot{z} + \omega^2 z =0 , \]
so the motion in the z direction is simple harmonic.

Of course these equations could have been found from the Lagrangian method.    ■

Example: Let us consider the problem of finding a curve on a vertical circular cylinder, a curve that evolves from a given point A and terminates on the vertical line that is the generator of the cylinder. Under the influence of gravity, a particle should move on this curve in the least time. We assume that the particle starts from rest, that the friction force between the particle and the cylinder is negligible, and that the terminal vertical line is located by a given angle θ1, as showen on the figure.

The particle has two degrees of freedom and we select the cylindrical coordinates (θ,z) as the generalized coordinates. Due to the fact that the motion of the particle is conservative, the minimum-time curve, called the brachistochrone, can be determined by minimizing the functional T. The kinetic energy and potential energies are given by

\[ \text K} = \frac{m}{2} \left( r^2 \dot{\theta}^2 + \dot{z}^2 \right) , \qquad \Pi = -mgz , \]
where r denotes the radius of the cylinder. The total energy of the particle is E = L + Π = 0. Taking the angle θ as an independent variable and consider the time functional
\[ T = \int_0^{\theta_1} \sqrt{\frac{r^2 + \left( z' \right)^2}{2g\,z}}\,{\text d}theta, \qquad z(0) = 0 , \]
where z(θ1) is not given, \( z' = {\text d}z/{\text d}\theta . \) The location of the terminal line CD is specified by the angle θ1, however the position of the point B on this line is not given. Thus, a boundary condition is missing, although the upper bound in the integral T is given. This is a typical case, and we will demonstrate that an appropriate boundary condition can be supplied in the course of the variational analysis. This is one of the characteristic features of variational calculus.    ■

 

  1. Chen, Y.M., Liu, J.K., 2007. A new method based on the harmonic balance method for nonlinear oscillators. Physics Letters A, 368:371–378. [doi:10.1016/j.physleta.2007.04.025]
  2. Olver, P.J., Nonlinear ordinary differential equations, University of Minnesota
  3. Papadakis, K.E., Goudas, C.L., and Katsiaris, G.A., The general solution of the Hénon--Heile problem, Astrophysics and Space Science, 2005, Vol. 295, Issue 3, pp. 375--396.
  4. Wu, B. S., Lim, C. W. and He, L. H., A new method for approximate analytical solutions to nonlinear oscillations of nonnatural systems, Nonlinear Dynamics, Vol. 32, 1–13, 2003.

 

Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions