\documentstyle[maple]{acmconf}
%%
%% <SHAREFILE=analysis/calcvar/calcvar.tex >
%% <DESCRIBE>
%% calcvar.tex     A paper documenting the design of the calcvar package.
%%              AUTHOR: Robert M. Corless, rcorless@uwo.ca
%% </DESCRIBE>
%% <UPDATE=R4 >

\title{Sufficiency Analysis for the Calculus of Variations}
\author{Robert M. Corless\\
Dept. Applied Mathematics\\
University of Western Ontario\\
London, CANADA N6A 5B7\\
rcorless@uwo.ca
}
\begin{document}
\maketitle
\begin{abstract}
Many of the computations in the calculus of variations are algebraic in
nature: computing the Euler-Lagrange equations and solving them, for
example.  But deciding whether or not the computed extremals provide
minima, maxima, or neither is an essentially analytic problem, and one
that has not been previously implemented in a computer algebra package.  I
describe here a partial attempt to do so in Maple, and detail some
successes and some failures of the attempt.  I believe the successes are
interesting, whilst some of 
the failures point to areas where Maple could be improved.
\end{abstract}
\section{Introduction}
This paper describes a Maple package for the calculus of variations,
which package is a completely revised version of
an Applied Math 511 
project by Yanan Wu, which was later presented as a
poster at ISSAC '92, Berkeley~\cite{Wu}.  
The present package is intended to {\sl help\/} the user
solve certain classes of problems from the calculus of variations, and in
particular is not intended as a `black box' for solving problems {\sl
for\/} the
user.  The package implements some necessary conditions, but most of the
interest for this context is from the partial implementation of some
sufficient conditions.
\section{Introduction to the Classical Theory}
This section will summarize the most elementary portions of the classical
theory of the calculus of variations, dealing with functions of one
independent variable.
The principal references used to prepare this section were 
\cite{Bliss,Clarke,Weinstock} and \cite{Young}.

\subsection{The Euler-Lagrange Equations}
The so-called `basic problem' of the calculus of variations is to find a
function $y(x)$ such that the functional
\[
  I(y) = \int_a^b f(t,y(t),y'(t))\, dt
  \]
is minimized over all smooth (or piecewise smooth) functions $y(t)$
satisfying certain boundary conditions---for example, $y(a) = A$ and $y(b)
= B$.  

% If $\hat y(x)$ is the smooth function at
% which the desired minimum of $I(y)$ occurs, and if $I(\hat y(x) + \e \eta(x))$ is
% defined for some arbitrary smooth function $eta(x)$ with $\eta(a) = 0$ and 
% $\eta(b) = 0$, for small enough $\e$, then
% \[
%   I(\hat y + \e \eta) = \int_a^b f(x,\hat y+\e\eta,\hat y'+\e\eta')\,dx
%   \]
% is now a function of $\e$, which must have a minimum at $\e = 0$.  In that
% case, if $I(\e)$ is smooth enough, we must have 
% \[
%   \frac{dI}{d\e}|_{\e=0} = \int_a^b f_y(x,\hat y,\hat y')\eta(x) 
%   + f_{y'}(x,\hat y,\hat y')\eta'(x)\,dx = 0\>.
%   \]
% If we integrate the second term by parts we get, using $\eta(a)= 0$
% and $\eta(b) = 0$,
% \[
%   \int_a^b \Biggl( f_y(x,\hat y,\hat y') - \frac{d}{dx} f_{y'}(x,\hat y,\hat y')
%            \Biggr)\eta(x)\,dx = 0\>.
%   \]
% One can then argue that since $\eta(x)$ was arbitrary and $\hat y$ is
% smooth, we must have the
% quantity in brackets identically zero.  This gives the {\sl
% Euler-Lagrange\/} equations:
If a smooth function $y$ exists which minimizes $I(y)$, and $I(y)$ is not
too `strange', then $y$ must satisfy the Euler-Lagrange equations:

\begin{equation}
  f_y(t, y, y') - \frac{d}{dt} f_{y'}(t, y, y')  = 0\>. \label{EuLag}
\end{equation}
This typically gives a second-order (nonlinear) boundary-value problem to
solve for the unknown function $ y(t)$, although for certain classes of
problems (where $f$ is independent of $t$, or independent of $y$, or
independent of $y'$) we can automatically reduce the problem to a
first-order differential equation (or first-order system of differential
equations if there is more than one dependent variable).

So the first problem for our Maple package is, given a functional $I$, 
find the Euler-Lagrange equations, and solve them if possible (such
solutions are called {\sl extremals\/} of the problem).

We list here several examples of this basic problem, together with some
refinements.
\begin{enumerate}
\item Geodesics.  Find the path from $(a,A)$ to $(b,B)$ which gives the
shortest distance
\begin{enumerate}
\item in the plane;  hence
\begin{equation}
  J(y) = \int_a^b \sqrt{ 1 + (y')^2}\,dx \label{plane}
\end{equation}
or, in the parametric case,
\begin{equation}
  J(x,y) = \int_{t_0}^{t_1} \sqrt{ \dot x^2 + \dot y^2 } \,dt
\label{parplane}
\end{equation}
\item on a sphere; hence
\begin{equation}
 J(\theta) = \int_{\phi_0}^{\phi_1} \sqrt{1 + \sin^2\phi
(d\theta/d\phi)^2}\,d\phi\>. \label{sphere}
\end{equation}
Now, in fact, the minimizing curves for both these examples can be found
by direct methods, without recourse to the classical theory of the
calculus of variations.  Nonetheless they provide useful examples as we
will see.
\end{enumerate}
\item The Brachistochrone.  This famous problem asks for the shape of a
wire joining the origin to a point $(a,A)$ along which a frictionless bead
can slide from the origin to $(a,A)$ in the least possible time.  If we
suppose that gravity acts in the positive $y$-direction, and that the
bead starts its journey with an initial velocity $v_0 = \sqrt{2gy_0}$ for
some $y_0$ (i.e. it is as if the particle fell from rest at $y_0$), then
the time of transit down a wire with shape $y(x)$ can be expressed as
\begin{equation}
T(y) = \int_0^a \sqrt{\frac{1 + (dy/dx)^2}{2g(y-y_0)} }\,dx 
\label{brachist}
\end{equation}
\item Queen Dido's Problem.  Here we are asked to 
find the shape of a boundary enclosing the
maximum area, for a given length of boundary.  
The classical formulation of this problem is that Queen Dido was granted
as much land as she could enclose in an oxhide; cutting the oxhide into
narrow strips and making a leather rope, she marked off a rather larger
area than was considered possible, and thus founded the city of Carthage.
In modern terms, the integral we are to minimize is, by the theory of
Lagrange multipliers,
\begin{equation}
A(x,y) = \int_{t_0}^{t_1} \frac{1}{2}(y\dot x - x\dot y) + \lambda
\sqrt{\dot x^2 + \dot y^2}\,dt 
\label{dido}
\end{equation}
\item Define the {\sl defect\/}~\cite{Defect} $\delta(y(x))$ by $\delta(y(x)) = y'(x) -
f(y(x))$, where $f(y(x))$ is some given function.
We are interested in the minimum 2-norm interpolant: that is, the 
function $y(x)$ with minimum 2-norm of the defect between the points
$y(t_k)$ and $y(t_{k+1})$, which come from a discrete numerical method
for solving $y' = f(y)$---that is, we wish to use the calculus of
variations to analyze the performance of the numerical method.  In that
case, the functional we wish to minimize is
\begin{equation}
I(y) = \int_{t_k}^{t_{k+1}} \left[y'(t) - f(y(t))\right]^T
\left[y'(t) - f(y(t))\right]\,dt
\end{equation}

\end{enumerate}
This is by no means an exhaustive list of problems in the calculus of
variations. Sturm-Liouville
problems can be cast as variational problems, 
Fermat's principle and  Hamilton's principle of least action lead to very
important classes of variational problems, and there are
many others~\cite{JJ}.

\subsection{Computation of extremals with Maple}
For example, here is how the package computes and solves the
Euler-Lagrange equations in the geodesic problem and the brachistochrone
problem.  We start with a small utility program which sets up some aliases
to make the output more readable.
\begin{maple}
> depends({x,y},t);
               I, x'', x', x, y'', y', y

> f := (dx^2 + dy^2)^(1/2);
                         2     2 1/2
                 f := (x'  + y' )

> Euler_Lagrange( f, t, [x,y]);
     y' (2 x' x'' + 2 y' y'')         y''
{1/2 ------------------------ - --------------,
             2     2 3/2           2     2 1/2
          (x'  + y' )           (x'  + y' )

        x' (2 x' x'' + 2 y' y'')         x''       
    1/2 ------------------------ - --------------, 
                2     2 3/2           2     2 1/2  
             (x'  + y' )           (x'  + y' )     


        x'                      y'
  -------------- = K[1],  -------------- = K[2]}
     2     2 1/2              2     2 1/2
  (x'  + y' )              (x'  + y' )
\end{maple}
Note that there are two first integrals here.
\begin{maple}
> first_integrals := select(has,",K);
                           x'               
 first_integrals := {-------------- = K[1], 
                        2     2 1/2         
                     (x'  + y' )            

		             y'
		       -------------- = K[2]}
		          2     2 1/2
		       (x'  + y' )

\end{maple}
Maple cannot solve these differential equations as given---they have an
infinite number of solutions, each a parameterization of the straight line
joining the boundary points (this is obvious to a human on dividing the
two first integrals to get $dy/dx = const$).  This is the first example of
where this package points out a possible area of improvement in computer
algebra systems:  I would like CAS to be able to solve differential
equations in {\sl parametric form\/}: however, I am unaware of any
existing theory to be implemented.  Most parametric solutions that I have
seen were derived on an {\sl ad-hoc\/} basis.

The next example, the brachistochrone, also raises this point, but in a
different way:
\begin{maple}
> f := (1 + dy^2)^(1/2)/y^(1/2);
                      2 1/2
               (1 + y' )
          f := ------------
                    1/2
                   y

> Euler_Lagrange(f,t,y);
              2 1/2          2          
       (1 + y' )           y'  y''      
{- 1/2 ------------ + ----------------- 
            3/2              2 3/2  1/2 
           y          (1 + y' )    y    

                2
              y'                    y''
 + 1/2 -----------------  -  -----------------,
              2 1/2  3/2            2 1/2  1/2
       (1 + y' )    y        (1 + y' )    y

         2 1/2            2
  (1 + y' )             y'
  ------------ - ----------------- = K[1]}
       1/2              2 1/2  1/2
      y          (1 + y' )    y

\end{maple}
One can verify that solutions of this are cycloids.  Maple
can find an implicit representation of the cycloids:
\begin{maple}
> select(has,",K);
               2 1/2            2
        (1 + y' )             y'
       {------------ - ----------------- = K[1]}
             1/2              2 1/2  1/2
            y          (1 + y' )    y

> dsolve(",y);
                
/  y      2\1/2 
|----- - y |    
|    2     |    
\K[1]      /    

                      2
         arcsin(2 K[1]  y - 1)
   - 1/2 --------------------- + t = _C1,
                     2
                 K[1]
                  
  /  y      2\1/2 
- |----- - y |    
  |    2     |    
  \K[1]      /    

                     2
        arcsin(2 K[1]  y - 1)
  + 1/2 --------------------- + t = _C1
                    2
                K[1]
\end{maple}
This is a success, of sorts.  The usual parametric representation of a
cycloid, $x = a(\theta - \sin \theta)$, $y = a(1-\cos\theta)$, is far
superior to this implicit representation, and it would be very nice if
Maple could find such a parametric solution of a given differential
equation.  Indeed, it would be nice to be able to find parametric
solutions of {\sl algebraic\/} equations, but again I know of no theory
to do this in a reasonably general way.  In \cite{Bif} a Gr\"obner basis
is used to help find the equilibria of a differential equation system:
however, the solution of that Gr\"obner basis proceeded by hand to a
parametric solution, which again was far superior to the mechanically
generated solution to the polynomial system.




\section{Sufficient Conditions}
It is possible for an extremal to provide a `global minimum', 
`weak local minimum', `strong local minimum', or analogously provide
a `maximum'; it is also possible that an extremal provides neither a
minimum nor a maximum but instead provides a sort of saddle point; 
it is possible that there will be {\sl
no\/} extremal joining $y(a)=A$ to $y(b)=B$; it is possible that even if
$y(x)$ provides neither minimum nor maximum it will {\sl look\/} to simple
tests as if it does; and, finally, it may be that the true minimum is
provided by a {\sl non\/}-smooth function.

Much of the research on the calculus of variations in the early part of
this century focussed on whether or not we could provide tools
for automatically deciding whether our extremal provides a global minimum,
a strong local minimum, a weak local minimum, a saddle, a weak local
maximum, a strong local maximum, or a global maximum.  
We describe some of the simplest such tools here.
\subsection{Convexity}
The best and simplest test for global minima (or maxima) is convexity
(concavity).  If we can show that our functional is convex (concave), 
then we are
done, and any extremals will automatically provide global minima
(maxima)~\cite{Clarke}.
The convexity test is itself very simple:  we examine the Hessian matrix
$H(f)$ of second partial derivatives 
and test to see if it is positive-, negative-, or in-definite (actually
semi-definite will do, as we will see).
This can be done by rational means, simply by doing Gaussian elimination
without pivoting (if pivoting is necessary then the matrix is indefinite).
It turns out that testing a matrix for definiteness is already part of
Maple's linear algebra package, and so this test is trivial to implement
in Maple.

Here we give the results of this convexity test for the geodesic
problem.  For reasons of space, we do not give the full
Hessian matrix of the geodesic problem, but simply look at the nonzero
submatrix:
\begin{maple}
> H := map(normal,linalg[submatrix](",3..4,3..4));
       [          2                         ]
       [        y'                y' x'     ]
       [  --------------   - -------------- ]
       [     2     2 3/2        2     2 3/2 ]
       [  (x'  + y' )        (x'  + y' )    ]
  H := [                                    ]
       [                            2       ]
       [        y' x'             x'        ]
       [ - --------------   --------------  ]
       [      2     2 3/2      2     2 3/2  ]
       [   (x'  + y' )      (x'  + y' )     ]

> linalg[definite](H,positive_semidef);
             2                         2
           y'                        x'
   - -------------- <= 0 and - -------------- <= 0
        2     2 3/2               2     2 3/2
     (x'  + y' )               (x'  + y' )

\end{maple}
Since $x$ and $y$ are real, we see that $H$ is positive semidefinite, and
hence the function is convex, and hence the extremals provide global
minima.  When this test succeeds, as it does here, it is the most powerful
of the sufficiency conditions.  In the package, one does not have to
separately compute the Hessian matrix and call \verb+definite+, but simply
call the routine \verb+Convex+ which does this for you.

\section{Legendre, Jacobi, and Weierstrass Conditions}
It is possible for extremals to give global minima even when the integrand
is not convex.  In such a case, the test of the previous section will
fail.  The following compound test for sufficiency consists of verifying
three strengthened necessary conditions, those of Legendre, Jacobi, and
Weierstrass.  The theory on which they are based is that of so-called
`fields of extremals'.  In practice, these tests can often be performed in
a mechanical way, without real reference to the field of extremals---that
is, it is often true that the inequality in the strengthened Weierstrass
condition holds for {\sl all\/} values of $p$, and not just for the
slope-field values of the theory: if this is true, then the test can be
performed in a mechanical way.  

The compound test is usually expressed as
\begin{itemize}
\item[I.] $y$ is an extremal ({\sl i.e.\/} a solution of the
Euler-Lagrange equation).
\item[${\rm II}_b'$.] The Weierstrass {\sl excess function\/} $E$ is positive,
where 
\[ 
  E(t,y,y',p) = f(t,y,p) - f(t,y,y') - (y'-p)f_{y'}(t,y,p)\>.
  \]
If $E>0$ for every arc $y$ which is different from the extremal, in a neighbourhood
of the extremal, where $p(t)$ is taken
as the slope of the particular extremal which passes through the point
$y(t)$ when the independent variable takes the value $t$, then the strengthened Weierstrass condition is said to be
satisfied.  In practice, this slope field is difficult to calculate---if,
however, the above inequality is true for {\sl all\/} $p$, then it is true
in particular for the slope field:  we will be concerned only with such
cases here.
\item[${\rm III}'$] (Legendre) At every element $(t,y,y')$ of the minimizing
arc, 
\[
  f_{y'y'} > 0\>.
  \]
This condition, strengthened to hold in a neighbourhood of the extremal,
is denoted $III_b'$.
\item[${\rm IV}'$] (Jacobi) There is no point conjugate to the initial point
anywhere on the minimizing arc between the initial point and the endpoint,
inclusive of the endpoint.  [We will define conjugacy shortly.]
\end{itemize}

If $y$ satisfies $I$, $III'$, and $IV'$, then $y$ is known to provide a
{\sl weak local minimum\/} (i.e. a minimum over all arcs close to $y$ with
derivatives also close to $y'$); if it satisfies $I$, $II_b'$, and $IV'$, $y$
is then known to provide a {\sl strong local minimum\/} (i.e. over all
arcs close enough to $y$, with no restriction on their derivatives).   
The conditions
$I$, $III_b'$, and $IV'$ also ensure that $y$ provides a strong local
minimum.

The Maple package provides tools for the computation of the Weierstrass
excess function and conjugate points: verification that the excess
function is positive {\sl is left to the user}, because that information
is essentially geometric and Maple is not particularly good at
inequalities and geometry.
The following demonstrates the Weierstrass excess function computation for
the bra\-chis\-to\-chrone prob\-lem.
\begin{maple}
> signum(y) := 1;
                    signum(y) := 1
\end{maple}
We will discuss the necessity of that statement after the example.
\begin{maple}
> Weierstrass(f,t,y,p);
         2 1/2            2 1/2
  (1 + y' )      (1 + p[1] )    
  ------------ - -------------- 
       1/2             1/2      
      y               y         


          (y' - p[1]) p[1]
      - -------------------
                 2 1/2  1/2
        (1 + p[1] )    y

> normal(");
            2 1/2          2 1/2
   - (1 + y' )    (1 + p[1] )    + 1 + p[1] y'
 - -------------------------------------------
                        2 1/2  1/2
               (1 + p[1] )    y
\end{maple}
Now, using the Cauchy-Schwartz inequality $|a| |b| \geq a\cdot b$, with
equality only if $a$ and $b$ are proportional, we can see that $E \geq 0$
(taking $a = (1,y')$ and $b = (1,p_1)$).  Further, we see that $E = 0$ if
and only if $y' = p_1$, and by uniqueness of the solutions of initial
value problems $E = 0$ only at the extremal.  It was not necessary to
compute the slopes of the field of extremals: the Weierstrass condition
holds for all $p$, not just the necessary ones.  Hence the brachistochrone
extremals provide minima over whatever region can be covered by a field of
extremals (in general this is determined by the location of conjugate
points).

I remark that the command \verb+signum(y) :=1+ was necessary here to allow
simplification of things like $\sqrt{ y^2 }$.  Maple's \verb+assume+
facility is not yet strong enough to allow us to say 
\begin{maple}
> assume(y>0)
\end{maple}
because here $y = y(t)$ (it prints as $y$ because of an alias).
Similarly, the routine \verb+Weierstrass+ assumes  $p_1$ 
is real
by saying
\begin{maple}
> signum(p[1]^2) := 1
\end{maple}
instead of \verb+assume(p[1],real)+, which is
the more natural and general syntax (but will not work, because
\verb+assume+ does not know about indexed names).  
It is hoped that these weaknesses
will be remedied for Maple V Release 4.

\subsection{Conjugate Points}
The usual way to find conjugate points is to first
find a linear differential equation, Jacobi's equation, and solve it. 
Then one looks for zeros of this solution.  However, forming Jacobi's
equation requires one to know the extremals for the problem already, and
if (as is usually the case) we know a basis for the solution space for the
Euler-Lagrange equations, then there is a simple way to find an algebraic
equation whose zeros (if any) are the conjugate points: using this
technique is much less computationally intensive.  Analytically, Jacobi's
equation is to be preferred, because using it sometimes one can determine the
existence of conjugate points without actually finding them.  But in the
present computational context, most of the time the algebraic equation
will be more useful.   See~\cite[p. 150]{Bliss} for a derivation of the
theory for this algebraic equation, which is just the determinant
\[
J(t) = \det \left[ \begin{array}{cc} 
y_a(t,a_0,b_0) & y_b(t,a_0,b_0) \\
y_a(t_1,a_0,b_0) & y_b(t_1,a_0,b_0) 
\end{array}
\right]\>,
  \]
where $t_1$ is the value of the independent variable at the initial point,
$a_0$ and $b_0$ are the values of the constants used in expressing the
extremal in terms of the basis of the solution space.  The zeros of $J(t)$
(not including $t_1$) are the points conjugate to $t_1$.

The following example shows the use of this in a Sturm-Liouville problem.
\begin{maple}
> depends({p,q,r,y},t);
   I, p'', p', p, r'', r', r, 
       q'', q', q, y'', y', y

> f := -p*dy^2 + (q-lambda*r)*y^2;
                2                   2
     f := - p y'  + (q - lambda r) y

> Euler_Lagrange(f,t,y);
    {2 (q - lambda r) y + 2 p' y' + 2 p y''}

> Weierstrass(f,t,y,P);
          2         2
    - p y'  + p P[1]  + 2 (y' - P[1]) p P[1]

> factor(");
                        2
         - p (P[1] - y')
\end{maple}
That shows that the solution to the Sturm-Liouville problem
$(py')' + qy = \lambda r y$, plus boundary conditions,
provides a strong local maximum to the functional
$\int_a^b f$.  In fact, one can show that the value
of this maximum is zero, and hence if we normalize
the eigenfunctions by $\int_a^b r y^2 dt = 1$, then the
value of the maximum of $\int_a^b -p(y')^2 + qy^2\,dt$ subject
to that normalization is just $\lambda$, the first eigenvalue
of the Sturm-Liouville problem.

The fact that $E \leq 0$ with equality only if $y'=P$ shows that a
`global' maximum occurs over the portion of the plane covered
by the field of extremals.  We cannot further investigate this
field without actually solving the Euler-Lagrange equations and
looking for conjugate points.
\begin{maple}
> f1 := subs(p=1,q=0,r=-1,f);
                2           2
      f1 := - y'  + lambda y

> Euler_Lagrange(",t,y);
                        2           2
 {2 lambda y + 2 y'', y'  + lambda y  = K[1]}

> select(has,",K);
        2           2
     {y'  + lambda y  = K[1]}

> dsolve(",y);
                   1/2            
           y lambda               
    arcsin(-----------)           
                 1/2              
             K[1]                 
  - ------------------- + t = _C1,
               1/2                
         lambda                   

                     1/2
             y lambda
      arcsin(-----------)
                   1/2
               K[1]
      ------------------- + t = _C1
                 1/2
           lambda
\end{maple}
Obviously if our boundary conditions are $y(0) = 0$ and  $y(\pi) = 0$, our
eigenfunctions are $y(t) = A \sin(nt)$, with eigenvalues $n^2 \geq 1$.

We can now look for conjugate points.

\begin{maple}
> conjugate_eqn(A*sin(t)+B*cos(t),
                [A,B],[1,0],t,0);
              sin(t)
\end{maple}

This equation is zero when $t=\pi$, which is the nearest conjugate
point to the origin in the positive $t$-direction.

The field of extremals is thus the strip $0 < t < Pi$, $-\infty<y<\infty$.

% \section{The Royal Road of Carath\'eodory}
% This test is the `weakest' sufficiency test of all---that is, this
% condition is both necessary {\sl and} sufficient.  It is also
% computational: we simply have to compute an integral, and verify an
% inequality.  This makes it a natural one to use in a computer algebra
% context.  Unfortunately, the integral calculation is often extremely
% difficult, as is the verification of the inequality.  We give only one
% example here.
\section{Concluding Remarks}
The present calculus of variations package provides a useful
functionality, and highlights some weaknesses of current computer algebra
systems.  In particular, it would be very helpful to have facilities for
solving algebraic and differential equations in parametric form, and also
for easy testing of classical inequalities.
\section*{Acknowledgements}  
I wish to thank Colin Clark and F.~H.~Clarke for introducing me to the calculus of
variations.  I wish also to thank my Applied Computer Algebra students for
many discussions.  This work was supported by NSERC.
\begin{thebibliography}{99}
\bibitem{Bliss} Gilbert Ames Bliss, {\sl Calculus of Variations}, Carus
Mathematical Monographs, MAA, 1925.
\bibitem{Clarke} Frank H. Clarke, {\sl Lecture Notes in the Calculus of
Variations and Optimal Control\/}, course notes for Math 426, University
of British Columbia, 1979.
\bibitem{SR} Robert M. Corless, {\sl Symbolic Recipes}, vol I,
Springer-Verlag, to appear 1994.
\bibitem{Bif} Robert M. Corless, ``Bifurcation in a flow-induced vibration
model'', to appear in Proc. Fields Institute Workshop on Normal Forms and
Homoclinic Chaos, W. Nagata and W. F. Langford, eds.
\bibitem{Defect} Robert M. Corless, ``Defect-Controlled Numerical Methods
and Shadowing for Chaotic Ordinary Differential Equations'', Physica D,
{\bf 60}, 1992.
\bibitem{JJ} Sir Harold Jeffreys and Bertha Swirles (Lady Jeffreys), {\sl
Methods of Mathematical Physics}, 3rd ed., Cambrige, 1956.
\bibitem{Weinstock} Robert Weinstock, {\sl Calculus of Variations\/},
Dover, 1974.
\bibitem{Wu} Yanan Wu and Robert M. Corless, ``Doing Calculus of
Variatiosn Symbolically'', poster at ISSAC '92, Berkeley.
\bibitem{Young} L. C. Young, {\sl Lectures on the Calculus of Variations
and Optimal Control\/}, Saunders, Philadelphia, 1969.
\end{thebibliography}
\end{document}


