3.42. 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.

3.42.1. Poisson Equation

The Lagrangian for Poisson equation is:

([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:

\delta L = 0\,,

which yields:

(\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 ( vanishes and we obtain:

(\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:

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

Now we substitute ( into ( and obtain:

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

Thus ( implies both ( and ( The equation ( holds for any \delta u(x) (generally not zero at the boundary) and thus it implies:

('(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 ( and thus also in ( This implies that ( does not hold at the Dirichlet part of the boundary and we have to set the value u(x) there directly.


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

([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 ( 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:

( := \{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 ( Then from calculus of variations:

( = 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:

\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:

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

The definition of the space U_0 in ( is derived from the definition of the space U in (

To compute the variation of L, we substitute ( into (, differentiate with respect to \varepsilon and then set \varepsilon=0 using (

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

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

(\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 ( holds for all \delta u \in U_0. From ( one obtains (as in (

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

The governing equation ( is the same:

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

The boundary term ( becomes (see (

(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 (, 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 ( 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 (


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 ( and the space U for the trial functions u\in U in (
  2. Define the 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 ( 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.

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:

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 (

(\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 ( holds for all w\in U_0, where

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:

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 ( as:

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

Strong Formulation

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

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

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

3.42.2. 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:

([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:

( &= \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:

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

Substituting ( into ( we obtain:

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

And we can see that ( implies both the equation ( and the boundary term ( The boundary term is zero for r=0, so it reduces to:

(\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 ( vanishes. Or we allow \delta R(\infty) to vary, and then ( 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 ( has no surface term.