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
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.
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(t_{0}).
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,
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
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):
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
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
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
This first order partial differential equation can be solved by the method of
characteristics. Consider the auxiliary first order scalar ordinary
differential equation
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 r_{i} =
(x_{i},y_{i}) in the plane
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
P_{x} and P_{y} are integrals of motion.
As argued above,
\[
\left\{ P_x , H \right\} = \left\{ P_y , H \right\} = 0
\]
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.)
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
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 = [ q_{1}, ... ,
q_{n} ] and p = [ p_{1}, ... ,
p_{n} ] 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
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.
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
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
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 x_{1} and x_{2}).
The dynamical solution \( \displaystyle x(t,E) \)
of the Hamilton's equations is ̄first expressed as
where x_{0} 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 x_{1} and x_{2}.
The period of oscillation T(E) is a function of energy alone
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
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
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
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
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
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
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):
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
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
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.
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:
The vector q represents the position vector for the mechanical system,
while Π(q) is the potential energy. The constants
m_{i} are masses, and each
\( p_i = m\,\dot{q}_i \)
represents the momentum variable associated with the position variable
q_{i}. The Hamiltonian function is the total energy, the first
term being the kinetic energy because the term
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 = (q_{1},q_{2}) for the
position of the second body. Newton’s laws, with a suitable normalization,
then yield the following differential equations
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 H_{0} = -1/2, and the angular momentum is
L_{0} = 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):
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
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_{θ}, p_{z}:
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
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
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.
■
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]
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