===============================
Variational Formulation of PDEs
===============================

Not every equation allows a variational formulation (e.g., Navier-Stokes or
Euler equations do not have such a formulation), but many equations have one,
and we explain how it works on several examples.

Poisson Equation
================

The Lagrangian for Poisson equation is:

.. math::
    :label: poisson0

    L[u] = \int_a^b \left[ \half u'^2(x) - f(x) u(x) \right] \d x
        - [g(x) u(x)]_a^b\,.

Important note: technically, as we will see below, this imposes the Neumann
boundary condition and 1D Poisson equation with two Neumann boundary conditions
does not have a unique solution. At least one Dirichlet boundary condition is
needed for a unique solution. For example with $u(a) = u_0$ and $u'(b) = g$ the
boundary term becomes just $-gu(b)$. However, for simplicity, we will show the
derivation with two Neumann boundary conditions first and we will discuss how
to impose the Dirichlet boundary condition later.

The variational formulation is:

.. math::

    \delta L = 0\,,

which yields:

.. math::
    :label: poisson1

    \delta L &= \int_a^b \left[ u'(x) \delta u'(x)
        - f(x) \delta u(x) \right] \d x
        - [g(x) \delta u(x)]_a^b

             &= \int_a^b \left[-u''(x) \delta u(x)
        - f(x) \delta u(x) \right] \d x
        + [u'(x) \delta u(x)]_a^b
        - [g(x) \delta u(x)]_a^b

             &= \int_a^b \left[-u''(x) - f(x)\right] \delta u(x) \d x
        + [(u'(x) - g(x)) \delta u(x)]_a^b

             &= 0\,,

where we applied integration by parts.
This equation holds for any $\delta u(x)$, and in particular it holds for
$\delta u(x) = 0$ at the boundary (i.e., for $\delta u(a) = 0$ and
$\delta u(b) = 0$). Then the boundary term in :eq:`poisson1` vanishes and we
obtain:

.. math::
    :label: poisson2

    \int_a^b \left[-u''(x) - f(x)\right] \delta u(x) \d x = 0\,,

This equation holds for any $\delta u(x)$ that is zero at the boundary, and
thus it implies:

.. math::
    :label: poisson3

    u''(x) + f(x) = 0\,.

Now we substitute :eq:`poisson3` into :eq:`poisson1` and obtain:

.. math::
    :label: poisson3_boundary

    [(u'(x) - g(x)) \delta u(x)]_a^b = 0\,.

Thus :eq:`poisson1` implies both :eq:`poisson3` and :eq:`poisson3_boundary`.
The equation :eq:`poisson3_boundary` holds for any $\delta u(x)$ (generally
not zero at the boundary) and thus it implies:

.. math::
    :label: poisson3_boundary2

    u'(x) - g(x) = 0

at the boundary. Thus $g(x)$ imposes the Neumann boundary condition, i.e.,
the value of the derivative $u'(x) = g(x)$ at the boundary. This condition
is imposed variationally.

To impose a Dirichlet boundary condition, we want to impose the value of
$u(x)=u_0(x)$ at the boundary for some constant $u_0(x)$. As such, $u(x)$
is not allowed to vary at that part of the boundary, which means that the
variation $\delta u(x) = 0$ at the boundary. So we restrict the variation
$\delta u(x)$ to be zero at the Dirichlet part of the boundary in
:eq:`poisson1` and thus also in :eq:`poisson3_boundary`. This implies that
:eq:`poisson3_boundary2` does not hold at the Dirichlet part of the boundary
and we have to set the value $u(x)$ there directly.

Example
-------

As a particular example, let $u(a) = u_0$ and $u'(b) = g$. Then the Lagrangian
:eq:`poisson0` becomes:

.. math::
    :label: poisson_example_L

    L[u] = \int_a^b \left[ \half u'^2(x) - f(x) u(x) \right] \d x - g u(b)\,.

We can explicitly define the space $U$ of all trial functions $u \in U$ that
one can choose (admissible) and substitute in :eq:`poisson_example_L` as
follows. We have to impose the Dirichlet condition $u(a) = u_0$ on the space
itself, and we also have to choose how smooth functions we want. For finite
element applications one typically chooses $H^1$ (i.e., values and first
derivatives are from $L^2$) and we obtain:

.. math::
    :label: poisson_example_UU

    U := \{u : u \in H^1(a,b), u(a)=u_0\}

Now we derive what space the variation $\delta u(x)$ belongs to. Let
$u_\textrm{min}$ be the solution (the extremum of the functional
:eq:`poisson_example_L`). Then from calculus of variations:

.. math::
    :label: poisson_example_u

    u = u_\textrm{min} + \varepsilon \delta u(x)

Here $u$ is called the trial function and $\delta u(x)$ is called the test
function. Both $u$ and $u_\textrm{min}$ are from the space $U$. Thus we can
compute:

.. math::

    \delta u(a) = {u(a) - u_\textrm{min}(a) \over \varepsilon}
    = {u_0 - u_0 \over \varepsilon} = 0\,.

In addition, both $u,u_\textrm{min}\in H^1(a,b)$, so also their difference
$u(x) - u_\textrm{min}(x)$ and thus also $\delta u(x)={u(x) - u_\textrm{min}(x)
\over \varepsilon}$ is from $H^1(a,b)$. There are no other conditions ($u(b)$
and $u_\textrm{min}(b)$ are generally different, so in general $\delta u(b) \ne
0$) and so $\delta u(x) \in U_0$ where the space $U_0$ is:

.. math::
    :label: poisson_example_U0

    U_0 := \{w : w \in H^1(a,b), w(a)=0\}\,.

The definition of the space $U_0$ in :eq:`poisson_example_U0` is derived from
the definition of the space $U$ in :eq:`poisson_example_UU`.

To compute the variation of $L$, we substitute :eq:`poisson_example_u` into
:eq:`poisson_example_L`, differentiate with respect to $\varepsilon$ and then
set $\varepsilon=0$ using :eq:`functional_deriv`:

.. math::

    \delta L[u] = \left.{\d\over\d\varepsilon}L[u_\textrm{min}+\varepsilon
        \delta u] \right|_{\varepsilon=0}

as was done in :eq:`poisson1` and one obtains the weak form (below we drop the
label $\textrm{min}$ from $u_\textrm{min}$ and just use $u$):

.. math::
    :label: poisson_example_weak_form

    \delta L[u] =
    \int_a^b \left[ u'(x) \delta u'(x) - f(x) \delta u(x) \right] \d x
        - g \delta u(b) = 0\,.

The task is to find such function $u\in U$ so that
:eq:`poisson_example_weak_form` holds for all $\delta u \in U_0$.  From
:eq:`poisson_example_weak_form` one obtains (as in :eq:`poisson1`):

.. math::
    :label: poisson_example_2

    \int_a^b \left[-u''(x) - f(x)\right] \delta u(x) \d x
        + (u'(b) - g) \delta u(b) = 0\,.

The governing equation :eq:`poisson3` is the same:

.. math::
    :label: poisson_example_strong

    u''(x) + f(x) = 0\,.

The boundary term :eq:`poisson3_boundary` becomes (see
:eq:`poisson_example_2`):

.. math::

    (u'(b) - g) \delta u(b) = 0\,.

Which implies $u'(b) = g$.

The Dirichlet boundary condition is part of the definition of the function
space :eq:`poisson_example_UU`, so all trial functions $u$ that one can choose
(admissible) and substitute in $L[u]$ must lie in $U$. From the derivation of
the space $U_0$ in :eq:`poisson_example_U0` we can see that since the value of
$u(a)$ is fixed, we always have $\delta u(a) = 0$; on the other hand, since
$u(b)$ is not fixed, in general we have $\delta u(b) \ne 0$.

The Neumann boundary condition is imposed variationally due to the surface term
in the weak form :eq:`poisson_example_weak_form`.

Summary
~~~~~~~

We have shown above that there are three equivalent formulations which
fully and uniquely determine the solution and boundary conditions (both
Dirichlet and Neumann):

1. Define the functional $L[u]$ in :eq:`poisson_example_L` and the space $U$
   for the trial functions $u\in U$ in :eq:`poisson_example_UU`.
2. Define the weak form :eq:`poisson_example_weak_form` and the two spaces $U$
   and $U_0$, where $u\in U$ and $\delta u \in U_0$.
3. Define the strong form :eq:`poisson_example_strong` and the boundary
   conditions $u(a) = u_0$ and $u'(b) = g$.

Let us write down the three formulations in detail.

Variational Formulation
~~~~~~~~~~~~~~~~~~~~~~~

The variational formulation is the formulation 1. above.

.. math::

    L[u] = \int_a^b \left[ \half u'^2(x) - f(x) u(x) \right] \d x - g u(b)\,.

The task is to find such $u\in U$ that extremizes this functional
($\delta L[u] = 0$), where:

.. math::

    U := \{u : u \in H^1(a,b), u(a)=u_0\}\,.

Weak Formulation
~~~~~~~~~~~~~~~~

Weak formulation is the formulation 2. above, and it is customary to
write $w(x) \equiv \delta u(x)$ in the weak form
:eq:`poisson_example_weak_form`:

.. math::
    :label: poisson_example_weak_form2

    \int_a^b \left[ u'(x) w'(x) - f(x) w(x) \right] \d x - g w(b) = 0\,.

The task is to find such $u\in U$ so that :eq:`poisson_example_weak_form2`
holds for all $w\in U_0$, where

.. math::

    U   & := \{u : u \in H^1(a,b), u(a)=u_0\}\,,

    U_0 & := \{w : w \in H^1(a,b), w(a)=0\}\,.

We can also define:

.. math::

    a(u,w) &= \int_a^b u'(x) w'(x) \d x\,,

    b(w)   &= \int_a^b f(x) w(x) \d x + g w(b)

and write :eq:`poisson_example_weak_form2` as:

.. math::

    a(u,w) = b(w)\,.

Strong Formulation
~~~~~~~~~~~~~~~~~~

Strong formulation is the formulation 3. above. We are solving the equation:

.. math::

    u''(x) + f(x) = 0

subject to boundary conditions $u(a) = u_0$ and $u'(b) = g$.

Radial Schrödinger Equation
===========================

The derivation is similar as for the Poisson equation, except that we have
$g(x) = 0$ based on physical reasoning (that we cannot set the derivative
to a given value, or, alternatively, that we require the operator to be
self-adjoint).

The Lagrangian for the radial Schrödinger equation is:

.. math::
    :label: schr_radial0

    L[R] = \int_0^\infty \left[\half R'^2(r)
        + \left(V(r) + {l(l+1)\over 2 r^2}\right) R^2(r) \right] r^2 \,\d r\,.

We minimize the Lagrangian subject to the normalization condition
$N[R] = \int_0^\infty R^2(r) r^2\, \d r = 1$ as follows:

.. math::
    :label: schr_radial1

    0 &= \delta (L - \epsilon (N-1))

    &= \delta \int_0^\infty \left[ \half r^2 R'^2
    + (r^2 V + \half l(l+1)) R^2 - \epsilon r^2R^2 \right] \,\d r =

    &= 2\int_0^\infty \left[ \half r^2 R'(\delta R)'
    + (r^2 V + \half l(l+1)) R\delta R - \epsilon r^2 R\delta R \right]
      \,\d r =

    &= 2\int_0^\infty \left[ -\half (r^2 R')'
    + (r^2 V + \half l(l+1)) R - \epsilon r^2 R\right]\delta R \,\d r
      + [r^2 R' \delta R]_0^\infty

This equation holds for any $\delta R(r)$, and so it also holds when we
restrict $\delta R(r) = 0$ on the boundary and the boundary term vanishes. Then
it implies the radial Schrödinger equation:

.. math::
    :label: schr_radial2

    -\half (r^2 R'(r))' + (r^2 V(r) + \half l(l+1)) R(r) = \epsilon r^2 R(r)

Substituting :eq:`schr_radial2` into :eq:`schr_radial1` we obtain:

.. math::
    :label: schr_radial_boundary

      [r^2 R' \delta R]_0^\infty = 0

And we can see that :eq:`schr_radial1` implies both the equation
:eq:`schr_radial2` and the boundary term :eq:`schr_radial_boundary`. The
boundary term is zero for $r=0$, so it reduces to:

.. math::
    :label: schr_radial_boundary2

    \lim_{r\to\infty} r^2 R'(r) \delta R(r) = 0

We can see that there is no natural condition at $r=0$, and for $r=\infty$ we
only have two possible options. Either we impose $\delta R(\infty) = 0$ and
obtain the Dirichlet condition and the boundary term
:eq:`schr_radial_boundary2` vanishes. Or we allow $\delta R(\infty)$ to vary,
and then :eq:`schr_radial_boundary2` implies $R'(\infty) = 0$.

Unlike for the Poisson equation we are not allowed to set $R'(\infty)$ to
anything other than zero, and that's why :eq:`schr_radial0` has no surface
term.