\documentclass[11pt, fleqn]{article}
%File creation date: Nov. 19, 2005
%second version (Wei): April, 8, 2006

\usepackage{amsmath,amsfonts}
\usepackage{graphicx}
\usepackage{color}
\usepackage{latexsym}
\usepackage{amssymb}

\setlength{\topmargin}{-.5in}      %topmargin is 1/2 inch (note negative value)
\setlength{\oddsidemargin}{0in}    %left margin is 1 inch on right-hand pages
\setlength{\evensidemargin}{0in}   %same for left-hand pages in twosided document
\setlength{\textwidth}{6.5in}      %leaves 1 inch for the right margin
\setlength{\textheight}{9in}       %9 inches reserved for the text.

\newcommand\comment[1]{}
%\newcommand\comment[1]{#1}
%\definecolor{light}{gray}{0.93}
\newcommand{\nn}{\nonumber}
\definecolor{light}{gray}{0.93}
\newcommand{\mybox}[1]{\begin{center}\colorbox{light}{\begin{minipage}[b]{5.2in}#1\end{minipage}}\end{center}}

\begin{document}

\title{\bf Computing Response Function from Linear Response Theory}


\author{Keonwook Kang, Chris Weinberger and Wei Cai \\
Department of Mechanical Engineering, Stanford University \\
Stanford, CA 94305-4040}

\title{A Short Essay on Variational Calculus}

%\date{January 13, 2005} % start date
%\date{April 8, 2006} % second version
\date{May 3, 2006} %third (final) version

\maketitle

\tableofcontents

\section{Definition of a functional}

The calculus of variations deals with {\it functionals}, which are
functions of a function, to put it simply.  For example, the
calculus of variations can be used to find an unknown function
that minimizes or maximizes a functional.
%
Let us begin our discussion with the definition of a functional.

We are familiar with the definition of a {\it function}.  It is a
mapping from one number (or a set of numbers) to another value.
For example,
\begin{equation}
 f(x) = x^2
\end{equation}
is a function, which maps $x=2$ to $f(x)=4$, and $x=3$ to
$f(x)=9$, etc.
%
On the other hand, a functional is a mapping from a function (or a
set of functions) to a value.  For example,
\begin{equation}
 I[y(x)] = \int_0^\pi \left[ y(x) \right]^2 \, dx
\end{equation}
is a functional.  Notice that we use square brackets for a
functional, to signify the fact that its argument is a function.
%
When $y(x)=x$,
\begin{equation}
 I[y(x)] = \int_0^\pi x^2 \, dx = \frac{\pi^3}{3}
\end{equation}
When $y(x)=\sin x$,
\begin{equation}
 I[y(x)] = \int_0^\pi \sin^2 x \, dx = \frac{\pi}{2}
\end{equation}
Therefore, functional $I[y(x)]$ maps $y(x)=x$ to $\pi^3/3$ and
maps $y(x)=\sin x$ to $\pi/2$.
%
Because an integral maps a function to a number, a functional
usually involves an integral.  The following form of functional
often appears in the calculus of variations,
\begin{equation}
 I[y(x)] = \int_{x_1}^{x_2} F(x,y(x),y'(x)) \, dx \label{eq:Iyx}
\end{equation}
%
where $y'(x)\equiv dy(x)/dx$.
%
Many engineering problems can be formulated into a variational
form, in which we need to find a function $y(x)$ that minimizes
(or maximizes) this functional, subject to certain boundary
conditions, such as $y(x_1)=a$ and $y(x_2)=b$.
%
Notice that while $I[y(x)]$ is a functional, $F(x,y,y')$ is simply
a function that maps three numbers into one number.





\section{First variation}

For a function $f(x)$, its {\it differential}, $df$, is how much
$f$ changes if its argument, $x$, changes by an infinitesimal
amount $dx$.  For example, when $f(x)=x^2$,
\begin{equation}
 d f = 2 x \, dx
\end{equation}
If $x$ is at a minimum (or maximum) of $f(x)$, then $df$ should be
zero for an infinitesimal change of $x$.  In this example, the
minimum occurs at $x=0$.

For a functional $I[y(x)]$, the corresponding term is its (first)
{\it variation}, $\delta I$.  $\delta I$ is how much $I$ changes
if its argument, the function $y(x)$, changes by an infinitesimal
amount $\delta y(x)$.
%
This is illustrated in Fig.~\ref{fig:firstvariation}, where
function $y(x)$ is shown in bold line.  If the function undergoes
a small change,
\begin{equation}
 \tilde{y}(x) = y(x) + \delta y(x)
\end{equation}
where $\delta y(x)$ is a small, continuous function, then the
variation of the functional is,
\begin{equation}
 \delta I = I[\tilde{y}(x)] - I[y(x)]
\end{equation}
Our task is obtain the explicit expression of $\delta I$.

\begin{figure}[h]
\begin{center}
\includegraphics[width=2.5in]{Figs/firstvariation}
\caption{Function $y(x)$ (bold line) and a small variation from
it, $\tilde{y}(x)=y(x) + \delta y(x)$.  $\delta y(x)$ is a small,
continuous function which satisfies the boundary conditions
$\delta y(x_1)=\delta y(x_2) = 0$. } \label{fig:firstvariation}
\end{center}
\end{figure}

Suppose that the functional $I[y(x)]$ takes the form of
Eq.(\ref{eq:Iyx}).  To obtain $\delta I$, let us define
\begin{equation}
 \delta y(x) = \epsilon \phi(x)
\end{equation}
where $\phi(x)$ is an arbitrary continuous function whose value is
of the order $1$ and $\epsilon$ is an infinitesimal number.  We
also require $\phi(x)$ to be sufficiently smooth so that
$\tilde{y}'(x)$ is well defined.\footnote{Another restriction on
$\phi(x)$ is that $\phi(x)$ and $\phi'(x)$ have the same order of
smallness.  In this case, the variation in $y$ is said to be
\emph{weak}~\cite{fox}.}
%
We can now Taylor expand $I[\tilde{y}(x)]$ in terms of $\epsilon$.
\begin{eqnarray}
I[\tilde{y}(x)]
  &=& \int _{x_1} ^{x_2}
    F(x, y(x)+\epsilon\phi(x), y^\prime(x) +\epsilon \phi'(x)) \, dx \nonumber \\
  &=& \int _{x_1} ^{x_2} \left[
    F(x, y(x), y^\prime(x)) +
    \frac{\partial F}{\partial y} \,\epsilon \phi(x) +
    \frac{\partial F}{\partial y'} \,\epsilon \phi'(x) +
    \mathcal{O}(\epsilon^2) \right]\, dx \nonumber\\
  &=& I[y(x)] + \epsilon\,\int _{x_1} ^{x_2} \left[
    \frac{\partial F}{\partial y} \, \phi(x) +
    \frac{\partial F}{\partial y'} \, \phi'(x) \right]\, dx
    + \mathcal{O}(\epsilon^2)  \label{eq:taylorI}
%
%&=& \int _{x_1} ^{x_2} \left [ F(x,y,y') + \frac{\partial F}{\partial y} \delta y + \frac{\partial F}{\partial y'} \delta y' \right. + {} \nonumber \\
%&& {} \left. \frac{1}{2!} \frac{\partial^2 F}{\partial y^2} (\delta y)^2 + \frac{1}{2!} \frac{\partial^2 F}{\partial y'^2} (\delta y')^2 + \frac{1}{2!} \frac{\partial^2 F}{\partial y \partial y'} \delta y \delta y' + \dots \right ] \nonumber \\
%&=& \int _{x_1} ^{x_2} \left [ F(x,y,y') + \epsilon \left( \frac{\partial F}{\partial y} \phi + \frac{\partial F}{\partial y'} \phi ' \right) \right. + {} \nonumber \\
%&& {} \left. \frac{1}{2!} \epsilon^{2} \left( \frac{\partial^2
%F}{\partial y^2} \phi^2 + \frac{\partial^2 F}{\partial y'^2}
%\phi'^2 + \frac{\partial^2 F}{\partial y \partial y'} \phi \phi'
%\right) + O(\epsilon^3) \right ] \nonumber
\end{eqnarray}
%
Keeping terms only up to the first order, we have,
\begin{eqnarray}
\delta I &=& \epsilon\, \int _{x_1} ^{x_2}
             \left[\frac{\partial F}{\partial y} \,\phi(x)
             + \frac{\partial F}{\partial y'} \,\phi'(x) \right] dx \nonumber \\
         &=&  \int _{x_1} ^{x_2}
         \left[\frac{\partial F}{\partial y} \,\delta y(x)
             + \frac{\partial F}{\partial y'} \,\delta y'(x) \right] dx
\end{eqnarray}
%
Usually, it is more desirable to express $\delta I$ solely in
terms of $\delta y(x)$, and get rid of $\delta y'(x)$ in the
expression.  This can be done by integration-by-parts.
%
\begin{eqnarray}
\delta I &=& \epsilon\, \int _{x_1} ^{x_2}
             \frac{\partial F}{\partial y} \, \phi(x)\, dx
             + \epsilon\, \int _{x_1} ^{x_2}
               \frac{\partial F}{\partial y'} \, d\phi(x) \nonumber \\
         &=&  \epsilon\, \int _{x_1} ^{x_2}
             \frac{\partial F}{\partial y} \,\phi(x)\, dx
             + \epsilon\, \left[\frac{\partial F}{\partial y'} \phi(x)
             \right]_{x_1}^{x_2}
             - \epsilon\, \int _{x_1} ^{x_2}
               \frac{d}{dx} \left(\frac{\partial F}{\partial y'}\right)
              \phi(x)  \, dx \label{eq:intepart}
\end{eqnarray}
If we constrain the value of function $y(x)$ at $x_1$ and $x_2$ to
be constant, i.e., apply boundary conditions
$\phi(x_1)=\phi(x_2)=0$, then,\footnote{The same conclusion can be
made if $y(x)$ satisfies periodic boundary conditions,
$y(x_1)=y(x_2)$, $y'(x_1)=y'(x_2)$.}
\begin{eqnarray}
\delta I &=& \epsilon\, \int _{x_1} ^{x_2}
             \left[ \frac{\partial F}{\partial y}
             -  \frac{d}{dx} \left(\frac{\partial F}{\partial y'}\right)
             \right] \phi(x) \, dx \nonumber \\
         &=& \int _{x_1} ^{x_2}
             \left[ \frac{\partial F}{\partial y}
             -  \frac{d}{dx} \left(\frac{\partial F}{\partial y'}\right)
             \right] \delta y(x) \, dx
\end{eqnarray}
%
Therefore, if $y(x)$ is a minimum (or a maximum) of functional
$I[y(x)]$, subjected to the constraint that $y(x_1)$ and $y(x_2)$
are fixed, then $y(x)$ must satisfy the differential equation,
\begin{equation}
 \frac{\partial F}{\partial y}
             -  \frac{d}{dx} \left(\frac{\partial F}{\partial
             y'}\right)  \label{eq:eulerlag}
 = 0
\end{equation}
for $x\in(x_1,x_2)$.  This is called the {\it Euler-Lagrange
equation}.

\mybox{ %
{\it Example.} \ \
%
What is the shortest path between two points $A$ and $B$, as shown
in in Fig.\ref{fig:firstvariation}?  Of course, the answer is a
straight line.  Here we will derive this result from variational
calculus.
%
Given a curve $y(x)$ connecting points $A$ and $B$, the length of
the curve is,
\begin{equation}
 L[y(x)] = \int _{x_1} ^{x_2} F(x,y(x),y'(x)) \, dx
  = \int _{x_1} ^{x_2} \sqrt{ 1 + (y'(x))^2 } \, dx
\end{equation}
%
Therefore, $L$ is a functional of $y(x)$.  The minimum of
functional $L$ must satisfy the Euler-Lagrange equation
(\ref{eq:eulerlag}).  Notice that in this case $F$ has no explicit
dependence on $x$ or $y(x)$, so that,
\begin{displaymath}
  \frac{\partial F}{\partial y} = 0 \qquad \textrm{and} \qquad
  \frac{\partial F}{\partial y'} = \frac{y'(x)}{\sqrt{1 + (y'(x))^2 }}
\end{displaymath}
The Euler-Lagrange equation becomes
\begin{displaymath}
  -\frac{d}{dx}\left(\frac{y'(x)}{\sqrt{1 + (y'(x))^2 }}\right) = 0
\end{displaymath}
The solution is that $y'(x)$=const, hence the shortest path is a
straight line.
%
}





\section{Essential and natural boundary conditions}

In the above derivation of the Euler-Lagrange equation, we assumed
that $\delta y(x_1)=\delta y(x_2)=0$.  This makes the second term
in Eq.~(\ref{eq:intepart}) vanish.
%
In this case, the function $y(x)$ satisfies the boundary
condition,
\begin{eqnarray}
 y(x_1) &=& y_1  \\
 y(x_2) &=& y_2
\end{eqnarray}
where $y_1$ and $y_2$ are constants.  This is called the {\it
essential} (or geometric) boundary condition.
%
This is the appropriate boundary condition for some applications,
such as the example discussed above.
%
However, in other applications, we may need to apply other types
of boundary conditions to the function $y(x)$.

If we still want the second term in Eq.~(\ref{eq:intepart}) to
vanish (so that we obtain the familiar Euler-Lagrange equation),
but allowing $\delta y(x_1)$ and $\delta y(x_2)$ to be non-zero,
then we need to have,
\begin{eqnarray}
 \left. \frac{\partial F}{\partial y'} \right|_{x=x_1} &=& 0 \\
 \left. \frac{\partial F}{\partial y'} \right|_{x=x_2} &=& 0
\end{eqnarray}
This is called a {\it natural} boundary condition.
%
A system may also have a {\it natural} boundary condition at one
end ($x=x_1$) and an {\it essential} boundary condition at the
other end ($x=x_2$).

If the functional involves higher derivatives of $y(x)$, such as,
\begin{equation}
 I[y(x)] = \int_{x_1}^{x_2} F(x,y(x),y'(x),y''(x)) \, dx
\end{equation}
then the minimum (or maximum) of $I[y(x)]$ satisfies the following
Euler-Lagrange equation,
\begin{equation}
  \frac{d ^2}{d x^2}\left(\frac{\partial F}{\partial y''}\right) -
  \frac{d}{d x}\left(\frac{\partial F}{\partial y'}\right) +
  \frac{\partial F}{\partial y} = 0
\end{equation}
provided that $y(x)$ is subjected to the following essential or
natural boundary conditions,
\begin{center}
\begin{tabular}{|l|l|}
  \hline
  essential  & natural  \\
  \hline \hline
  $y(x_1)=y_1$   & $\frac{\partial F}{\partial y'}-\frac{d}{dx}
    \left(\frac{\partial F}{\partial y''}\right)= 0$ at $x=x_1$  \\
  \hline
  $y'(x_1)=y'_1$   & $\frac{\partial F}{\partial y''}= 0$ at $x=x_1$  \\
  \hline
  \hline
  $y(x_2)=y_2$   & $\frac{\partial F}{\partial y'}-\frac{d}{dx}
    \left(\frac{\partial F}{\partial y''}\right)= 0$ at $x=x_2$  \\
  \hline
  $y'(x_2)=y'_2$   & $\frac{\partial F}{\partial y''}= 0$ at $x=x_2$  \\
  \hline
\end{tabular}
\end{center}
%
In general, when the functional $I$ contains higher derivatives of
$y(x)$,
\begin{equation}
 I[y(x)] = \int_{x_1}^{x_2} F(x,y,y',y'',\cdots,y^{(n)}) \, dx
\end{equation}
the Euler-Lagrange equation is
\begin{displaymath}
  \sum_{i=0}^{n} (-1)^i \frac{d^{i}}{dx^{i}}
    \left( \frac{\partial F}{\partial y^{(i)}} \right) = 0
  \end{displaymath}
when appropriate boundary condition is applied.



\section{Second variation}

When $y(x)$ satisfies the Euler-Lagrange equation, the first
variation $\delta I$ vanishes and $y(x)$ is an extremum of the
functional $I[y(x)]$.
%
To find out whether $y(x)$ is a minimum or a maximum of $I[y(x)]$,
we need to look at the second variation of $I$.
%
Let us continue the Taylor expansion in Eq.~(\ref{eq:taylorI}) to
the second order,
\begin{eqnarray}
I[\tilde{y}(x)]
  &=& \int _{x_1} ^{x_2}
    F(x, y(x)+\epsilon\phi(x), y^\prime(x) +\epsilon \phi'(x)) \, dx \nonumber \\
  &=& I[y(x)] + \epsilon\,\int _{x_1} ^{x_2} \left[
    \frac{\partial F}{\partial y} \, \phi(x) +
    \frac{\partial F}{\partial y'} \, \phi'(x) \right]\, dx
    \nonumber\\
   && + \frac{\epsilon^2}{2}\,\int _{x_1} ^{x_2} \left[
    \frac{\partial^2 F}{\partial y^2}\, \phi^2(x) +
    \frac{\partial^2 F}{\partial y'^2}\, \phi'^2(x) +
    \frac{\partial^2 F}{\partial y\partial y'}\, \phi(x)\phi'(x)
    \right] \, dx
    + \mathcal{O}(\epsilon^3)  \label{eq:taylorI2}
\end{eqnarray}
%
The second variation of $I$ is the term that is of the second
order of $\epsilon$,
\begin{eqnarray}
\delta ^2 I &=& \frac{\epsilon^2}{2}\,\int _{x_1} ^{x_2} \left[
    \frac{\partial^2 F}{\partial y^2}\, \phi^2(x) +
    \frac{\partial^2 F}{\partial y'^2}\, \phi'^2(x) +
    \frac{\partial^2 F}{\partial y\partial y'}\, \phi(x)\phi'(x)
    \right] \, dx \nonumber\\
   &=& \frac{1}{2}\,\int _{x_1} ^{x_2} \left[
       \frac{\partial^2 F}{\partial y^2}\, \delta y^2(x) +
       \frac{\partial^2 F}{\partial y'^2}\, \delta y'^2(x) +
       \frac{\partial^2 F}{\partial y\partial y'}\, \delta y(x)\delta y'(x)
       \right] \, dx
\end{eqnarray}
When the first variation $\delta I$ vanishes, $y(x)$ is a maximum
if the second variation $\delta^2 I > 0$ (for arbitrary
$\phi(x)$), or a minimum if $\delta^2 I < 0$.  If $\delta^2 I =
0$, we will have to examine the higher order variations of $I$ to
determine whether $y(x)$ is a minimum or a maximum.
%
For more discussions see~\cite{fox,gel:fom}.




\section{Functional derivative}

For a function with multiple arguments, $f(x_1,x_2,\cdots,x_n)$,
if the differential $df$ can be written as,
\begin{equation}
 df = \sum_{i=1}^n g_i(x_1,\cdots,x_n)\,dx_i
\end{equation}
then $g_i(x_1,\cdots,x_n)$ is called the (partial) {\it
derivative} of $f$ with respect to $x_i$, for $i=1,\cdots,n$,
\begin{equation}
 \frac{\partial f}{\partial x_i} = g_i(x_1,\cdots,x_n)
\end{equation}
Notice that the derivative of a function $f(x_1,\cdots,x_n)$ with
respect to $x_i$ is another function $g_i(x_1,\cdots,x_n)$.

\vspace{0.1in}

Similarly, if the (first) variation of a functional
\begin{equation}
 I[y(x)] = \int_a^b F(y(x),y'(x)) \, dx
\end{equation}
can be written as,
\begin{equation}
 \delta I = \int_{a}^{b} g(x) \, \delta y(x) \, dx
\end{equation}
then the {\it functional derivative} of $I$ is
\begin{equation}
 \frac{\delta I}{\delta y(x)} = g(x)
\end{equation}
The functional derivative of a functional $I[y(x)]$ is a function
$g(x)$.
%
In the following, we discuss the connection between the (ordinary)
derivative and the functional derivative.

Imagine that we approximate function $y(x)$ by a piece-wise linear
curve that passes through a set of points $(x_i,y_i)$,
$i=0,\cdots,n+1$, where $x_i=i\Delta x$, $y_i=y(x_i)$, $x_0=a$,
and $x_{n+1}=b$, as shown in Fig.~\ref{fig:discretevar}.  Let
$y(x)$ subject to essential boundary conditions, so that
$y_0=y(a)$ and $y_{n+1}=y(b)$ are fixed.
%
Then the functional $I[y(x)]$ can be approximated by a sum,
\begin{eqnarray}
 I[y(x)] &=& \int_a^b F(y(x),y'(x)) \, dx  \nonumber\\
     &=& \lim_{n\to \infty} \sum _{i=0}^n
         F\!\left(y_i, \frac{y_{i+1}-y_i}{\Delta x}\right)
         \, \Delta x \nonumber\\
     &\equiv& \lim_{n\to \infty} I_n(y_1,\cdots,y_n)
\end{eqnarray}
where
\begin{eqnarray}
I_n(y_1, y_2, \cdots, y_n) &\equiv& \sum _{i=0} ^n F\!\left(
%      y_i,\frac{y_{i+1}-y_i}{\Delta x}  \right)\Delta x
      y_i,z_i  \right)\Delta x \\
z_i &\equiv& \frac{y_{i+1}-y_i}{\Delta x}
\end{eqnarray}

\begin{figure}[thb]
\begin{center}
\includegraphics[scale = 1]{Figs/discretevar}
\end{center}
\caption{A piece-wise linear approximation of function $y(x)$.
$x_i=i\Delta x$, $y_i=y(x_i)$, $x_0=a$, $x_{n+1}=b$.  Consider
that one of the node, $(x_k,y_k)$ moves vertically by $\delta
y_k$, the variational derivative $\delta I/\delta y$ at $x=x_k$ is
the change of $I$ divided by the change of area $\delta y\Delta x$
(shade region) in the limit of $\Delta x\to 0$.}
\label{fig:discretevar}
\end{figure}


Therefore, functional $I[y(x)]$ is the limit of an ordinary
function $I_n(y_1,\cdots,y_n)$ when the number of its arguments,
$n$, goes to infinity.
%
The partial derivative of $I_n$ with respect to one of its
arguments is,
%
\begin{eqnarray}
\frac{\partial I_n}{\partial y_k} &=& \left(
   \left.\frac{\partial F}{\partial y}\right|_{y=y_k}
 + \frac{1}{\Delta x}\, \left.\frac{\partial F}{\partial y'}\right|_{y'=z_{k-1}}
 - \frac{1}{\Delta x}\, \left.\frac{\partial F}{\partial y'}\right|_{y'=z_{k}}
   \right) \, \Delta x  \nonumber\\
 &\approx&  \left[
   \left.\frac{\partial F}{\partial y}\right|_{y=y_k}
  - \left. \frac{d}{dx} \left( \frac{\partial F}{\partial y'}
    \right) \right|_{x=x_k}
   \right] \, \Delta x
\end{eqnarray}
%
Compare this with the definition of the functional derivative, we
have,
\begin{eqnarray}
\frac{\delta I}{\delta y(x)} = \frac{\partial F}{\partial y} -
\frac{d}{dx} \left( \frac{\partial F}{\partial y'} \right) =
\lim_{\Delta x \to 0} \,
    \frac{\partial I_n}{\partial y_k}\, \frac{1}{\Delta x}
\end{eqnarray}
%
Therefore, the functional derivative $\delta I/\delta y(x)$ can be
interpreted as the following. Imagine a small and localized
perturbation of $y(x)$ near point $x$, $\delta I/\delta y(x)$ is
the change of $I$ divided by the change of area under the curve
$y(x)$.







\section{Examples}

\subsection{Lagrange's equation of classical mechanics}

The equations of motion for particles in classical mechanics can
be obtained from Hamilton's principle or the {\it principle of
least action}.  Let $q$ be the coordinate of a particle, $t$ be
the time, so that $q(t)$ is the trajectory of this particle.
Suppose that the particle's position at times $t_1$ and $t_2$ are
known,
\begin{eqnarray}
 q(t_1) &=& q_1 \\
 q(t_2) &=& q_2
\end{eqnarray}
i.e., $q(t)$ is subjected to essential boundary conditions.  The
question is: what trajectory $q(t)$ would the particle take to go
through points $q_1$ and $q_2$ exactly at times $t_1$ and $t_2$?

The Hamilton's principle states that the real trajectory of the
particle is the one that minimizes (or maximizes) the {\it action}
functional~\cite{lan:lif},
\begin{equation}
 S[q(t)] = \int_{t_{1}}^{t_{2}} L(t,q(t),q'(t)) \, dt \label{eqn:action}
\end{equation}
where $L(t,q(t),q'(t))$ is the {\it Lagrangian} of the system,
which is usually the kinetic energy minus the potential energy,
\begin{equation}
 L(q,q') = \frac{1}{2}m (q'(t))^2 - V(q) \label{eq:lagran}
\end{equation}
The real trajectory of the particle must have $\delta S=0$, and
hence satisfies the Euler-Lagrange equation,
\begin{equation}
\frac{d}{dt} \left( \frac{\partial L}{\partial q'} \right) -
\frac{\partial L}{\partial q} = 0.
\end{equation}
When the Lagrangian takes the form of Eq.~(\ref{eq:lagran}), this
leads to the following equation of motion,
\begin{equation}
 m q''(t) = - \frac{\partial V}{\partial q}
\end{equation}
which is identical to the Newton's equation of motion $F=ma$,
where $F=- \frac{\partial V}{\partial q}$ is the force, $a=q''(t)$
is the acceleration.
%
In general, if the system has $n$ degrees of freedom, $q_1,
\cdots, q_n$, the Lagrange's equations of motion are
\begin{displaymath}
\frac{d}{dt}\left( \frac{\partial L}{\partial q'_i}\right) -
\frac{\partial L}{\partial q_i} = 0 \qquad \textrm{for all $i$ =
1,2,$\cdots$,$n$}.
\end{displaymath}
%They are Lagrange's equations of motion.


\subsection{Tight string under lateral loading}
%\footnote{This example is excerpted from Kwon's book.\cite{kwon}}}

When a system subjected to external loading reaches thermal
equilibrium, its {\it Gibbs free energy} is minimized.  This
thermodynamic principle allows us to determine the equilibrium
shape of an elastic body subjected to external loading.

For example, consider an elastic string subjected to lateral
loading $w(x)$, as shown in Fig.~\ref{fig:stringload}.  Ignoring
the entropic contribution, the Gibbs free energy of the system is
simply the enthalpy $H$, i.e. the store elastic energy minus the
work done by the external load.  Assuming the string has a line
tension $T$, the enthalpy can be written as the following
functional of the shape of the string $y(x)$, (in the limit of
$y'(x)\ll 1$)
\begin{equation}
 H[y(x)] = \int_0^L \frac{1}{2}\, T (y'(x))^2 - w(x)\,y(x) \, dx
\end{equation}
%
Therefore, the variational derivative is
\begin{equation}
 \frac{\delta H}{\delta y(x)} = -T\, y''(x) - w(x)
\end{equation}
%
Because the equilibrium shape of the string minimizes the
functional $H[y(x)]$, it satisfies the following Euler-Lagrange
equation,
\begin{equation}
 T\, y''(x) + w(x) = 0
\end{equation}

\begin{figure}[hbt]
\begin{center}
\includegraphics[width=3in]{Figs/stringload}
\end{center}
\caption{Distributed load $w(x)$ is applied to the string
  hanging tightly between two walls. The deflection $y(x)$ and its derivative
  is assumed to be very small because the string is tight. }
\label{fig:stringload}
\end{figure}


\subsection{Phase field equation for microstructure evolution}

Many systems evolves in the direction that minimizes its Gibbs
free energy.  This observation can be used to construct evolution
equations of a wide range of systems, such as the phase field
model of material microstructure~\cite{csd}.
%
For example, consider a simple model of crystal growth by
solidification from the melt.  The system may be described by a
phase field, $\phi(x)$, such that $\phi(x)\approx 1$ if the
material at $x$ is solid and $\phi(x)\approx -1$ if it is liquid.
%
The free energy $F$ of the system is a functional of the phase
field $\phi(x)$, which can be written as,
\begin{equation}
 F[\phi(x)] = \int_0^L U(\phi(x)-1)^2 (\phi(x)+1)^2 + \epsilon\,
 (\phi'(x))^2 \, dx
\end{equation}
The first term in the integrand accounts for the system's
preference in either one of the two states, $\phi=\pm1$.  The
second term in the integrand introduces a free energy penalty for
the liquid-solid interface.
%
A simple evolution equation for $\phi(x)$ is that
$\dot{\phi}(x)\equiv\frac{d}{dt}\phi(x)$ follows the steepest
descent direction, i.e. the direction of the variational
derivative,
\begin{equation}
 \dot{\phi}(x) = - k \frac{\delta F}{\delta \phi(x)}
  = -k \left[ \, 4 U \phi(x)\, (\phi(x)^2-1)
            -\epsilon\, \phi''(x) \right]
\end{equation}
where $k$ is a rate constant.




\subsection{Schr\"odinger's equation of quantum mechanics}

An electron in quantum mechanics is described by its wave function
$\psi(x)$.  For a single electron moving in a potential $V(x)$,
the Hamiltonian operator is,
\begin{equation}
 \hat{H}  = -\frac{\hbar^2}{2m} \frac{\partial^2}{\partial x^2} + V(x)
\end{equation}
This means that,
\begin{equation}
 \hat{H} \psi(x) = -\frac{\hbar^2}{2m}
 \frac{\partial^2 \psi(x) }{\partial x^2} + V(x) \psi(x)
\end{equation}
%
The energy $E$ of the electron is a functional of its wave
function $\psi(x)$,
\begin{equation}
 E[\psi(x)] = \int_{-\infty}^{\infty} \psi^*(x) \hat{H} \psi(x) \,
 dx
\end{equation}
For simplicity, we will assume that $\psi(x)$ is real, so that its
complex conjugate $\psi^*(x)=\psi(x)$.  We will also assume that
$\psi(-\infty)=\psi(\infty)=0$, i.e. assume essential boundary
conditions.

The ground state wave function is the minimizer of functional
$E[\psi(x)]$.  However, the wave function must satisfy the
normalization condition,
\begin{equation}
 \int_{-\infty}^{\infty} \psi^*(x) \psi(x) \, dx  = 1
\end{equation}
Therefore, the ground state wave function satisfies the following
Euler-Lagrange equation,
\begin{equation}
 \frac{\delta E'}{\delta \psi(x)} = 0
\end{equation}
where
\begin{equation}
 E'[\psi(x)] = E[\psi(x)] - \lambda \int_{-\infty}^{\infty} \psi^*(x)
 \psi(x) \, dx
\end{equation}
and $\lambda$ is the Lagrange multiplier.
%
The explicit expression of the variational derivative is,
\begin{equation}
 \frac{\delta E'}{\delta \psi(x)} =
   2\left[ -\frac{\hbar^2}{2m}\frac{\partial^2 \psi(x) }{\partial x^2}
    + V(x)\,\psi(x) - \lambda\, \psi(x) \right]
\end{equation}
Therefore, the ground state wave function satisfies the following
eigen equation,
\begin{equation}
 \hat{H} \psi(x) = \lambda \psi(x)
\end{equation}
We can further show that $\lambda=E[\psi(x)]$ in this case, so
that the constant $\lambda$ is usually written as $E$. This leads
to the time-independent Schr\"odinger's equation,
\begin{equation}
 \left[ -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x)
 \right] \psi(x) = E \, \psi(x)
\end{equation}


%how to do Hartree-Fock, effective potential, can be derived from
% variational calculus.



\begin{thebibliography}{99}
  \bibitem{sch} R. S. Schechter, {\it The Variational Method in
  Engineering}, (McGrow-Hill Book Company, 1967).
  \bibitem{fox} C. Fox, {\it An Introduction to the Calculus of
  Variations}, (Oxford University Press, 1950).
  \bibitem{gel:fom} I. M. Gelfand and S. V. Fomin, {\it Calculus
  of Variations}, (Dover, New York, 2000).
  \bibitem{bru} B. van Brunt, {\it The Calculus of Variations},
  (Springer, 2004).
  \bibitem{lan:lif} L. D. Landau and E. M. Lifshitz,
  {\it Mechanics}, (Pergamon Press, 1969).
  \bibitem{lanczos} C. Lanczos, {\it The Variational Principles of
  Mechanics}, (Dover, New York, 1970).h 9
  \bibitem{csd} V. V. Bulatov and W. Cai, {\it Computer Simulation
  of Dislocations}, (Oxford University Press, to be published,
  2006),   chapter 11.
  %\bibitem{kwon} T. H. Kwon, ``Introduction to Finite Element
  %Method'', Unpublished
  %\bibitem{dav} A. Davol, Course notes on ``Finite Element Method''
  %\bibitem{math} MathWorld
  %http://mathworld.wolfram.com/FunctionalDerivative.html
  %\bibitem{math} J. P. Lewis's note on ``Variational Calculus and
  %Discrete Alternatives'' (http://www.idiom.com/~zilla/Work/Notes/varcalce.pdf)
\end{thebibliography}
\end{document}






















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% old text, no longer used
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\comment{

%\section{Calculus of Variations}

From a mathematical point of view, the Calculus of Variations is
the branch of mathematics that deals with finding an unknown
function which minimizes or maximizes a functional. A functional
is a mapping from a set of functions to the real space, so it can
be thought of as a function that has functions for arguments.
%in the integrand of an integral which is to be made to have some
%extreme value, either maximum or minimum.
In ordinary differential calculus, the derivative is defined in
terms of the infinitesimal change of a function with respect to
its argument.  The first derivative test can be used to determine
extreme values of the function, and the second derivative test can
differentiate between local minima, maxima, and points of
inflection.  However, in Calculus of Variations, the problems
involve changes of the functional with respect to its function.
Similarly, the first variation determines extremals of the
functional and the second variation can determine if they are
minimums, maximums, or other extremals.  For example, we can think
of the following statement as one of the simplest variational
problems:\\
%\begin{tabular}{|l|}
%  \hline

Find \footnote{This is sometimes referred to as the
\emph{Variational Form} of the problem.} a function $y(x)$ which
minimizes (or maximizes) a functional $I(y)$ such that
\begin{equation}
I(y) =\int _{x_1} ^{x_2} F[x, y(x), y^\prime(x) ] dx.
\label{eqn:functional}
\end{equation}
%subject to the boundary conditions $F(x_1)=a$ and $F(x_2)=b$.
subject to the boundary conditions $y(x_1)=a$ and $y(x_2)=b$.
%  \hline
%\end{tabular}


\subsection{Variation operator $\delta$}
The variation symbol $\delta$ is used to express an infinitesimal
change or a small fluctuation of a function. If $y(x)$ is the
exact solution of the statement (\ref{eqn:functional}) and
$\tilde{y}(x)$ is an approximate, then the variation of $y(x)$ can
be expressed using $\delta$ operator as
\begin{displaymath}
  \delta y = \tilde{y}(x) - y(x) = \epsilon \phi(x)
\end{displaymath}
where $\epsilon$, the variable parameter, is an arbitrarily very
small constant independent of $x$ and $y$.
%where $\epsilon$, the variable parameter, is a constant independent of
%$x$ and $y$ and has a very small value($\ll 1$).
The function $\phi(x)$ can be any arbitrary but admissible
function.  Admissability usually requires that the function is
sufficiently continuous and differentiable such that equation
\ref{eqn:functional} is well defined.  An admissible function must
also satisfy certain boundary conditions. For this example the
boundary is fixed which requires that $\phi(x_1)=\phi(x_2)=0$.
Boundary conditions and their effects will be discussed more in
detail later.




\subsection{Properties of the operator $\delta$}
\flushleft
\begin{enumerate}
  \item There is no variation in $x$, because the function $y(x)$ has
  an infinitesimal change at fixed $x$. Or
    \begin{equation}
      \delta x = 0
    \end{equation}
  \item The $\delta$ operator is commutative with differential and
  integral operators.
    \begin{eqnarray}
      \delta  \frac{dy}{dx} &=& \frac{d}{dx}\delta y \\
      \delta \int y dx &=& \int \delta y dx
    \end{eqnarray}
  \item Chain rule
    \begin{eqnarray}
      \delta F[x,y(x),y^\prime(x)] &=& F(x, \tilde{y}(x), \tilde{y}^\prime (x))
                 - F(x, y(x), y^\prime(x)) \nonumber \\
               &=& F(x, y(x) + \epsilon \phi(x), y^\prime(x) + \epsilon \phi(x)^\prime)
                 - F(x, y(x), y^\prime(x)) \nonumber \\
               &=& \frac{\partial F}{\partial y}\delta y +
                 \frac{\partial F}{\partial y^\prime}\delta y^\prime
    \end{eqnarray}
  \item The $\delta$ operator is linear. For any real numbers $a$ and
  $b$
    \begin{eqnarray}
      \delta (aF + bG)[x,y,y^\prime,\cdots] &=&
                        (a\tilde{F} + b\tilde{G}) - (aF + bG) \nonumber \\
           &=& a(\tilde{F}-F) + b(\tilde{G}-G) \nonumber \\
           &=& a\delta F[x,y,y^\prime,\cdots] + b\delta G[x, y, y^\prime,\cdots]
    \end{eqnarray}
      %where $J = \int _{x_1} ^{x_2} F[x, y(x), y^\prime (x)] dx$ and
      %$K = \int _{x_1} ^{x_2} G[x, y(x), y^\prime (x)] dx$.
  \item Product rule
    \begin{equation}
      \delta (FG) = F \delta G + G \delta F
    \end{equation}
  \item Since $\delta y$ and $\delta y'$ are functions of $x$,
    \begin{equation}
      \frac{\partial (\delta y)}{\partial y} = 0 \qquad
      \frac{\partial (\delta y')}{\partial y} = 0 \qquad
      \frac{\partial (\delta y)}{\partial y'} = 0 \qquad
      \frac{\partial (\delta y')}{\partial y'} = 0
    \end{equation}
  \item $\delta (F^n) = nF^{(n-1)}\delta F$

\end{enumerate}

}


\comment{

The Eqn.(\ref{eqn:func_deriv}) is actually the left-hand side of
the Euler-Lagrange Eqn.(\ref{eqn:strong}). The denominator
$\partial y_k \Delta x$ in the Eqn.(\ref{eqn:func_deriv}) means
the area change in $x$-$y$ plane according to the variation in
$y$. (See Fig.\ref{fig:fig3}.) Generally, the functional
derivative $\delta I(y(x))/ \delta y(x)$ can be said as the
limiting case of the the variation in $I$ over the area
change($\delta y \Delta x$) due to the variation in $y$ around a
fixed $x$.

}
