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:

(3.42.1.1)¶ $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:

$\delta L = 0\,,$

which yields:

(3.42.1.2)¶ $\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 (3.42.1.2) vanishes and we obtain:

(3.42.1.3)¶ $\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:

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

Now we substitute (3.42.1.4) into (3.42.1.2) and obtain:

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

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

(3.42.1.6)¶ $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 (3.42.1.2) and thus also in (3.42.1.5). This implies that (3.42.1.6) 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 (3.42.1.1) becomes:

(3.42.1.7)¶ $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 (3.42.1.7) 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:

(3.42.1.8)¶ $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 (3.42.1.7)). Then from calculus of variations:

(3.42.1.9)¶ $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:

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

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

The definition of the space $U_0$ in (3.42.1.10) is derived from the definition of the space $U$ in (3.42.1.8).

To compute the variation of $L$ , we substitute (3.42.1.9) into (3.42.1.7), differentiate with respect to $\varepsilon$ and then set $\varepsilon=0$ using (3.17.3.2):

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

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

(3.42.1.11)¶ $\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 (3.42.1.11) holds for all $\delta u \in U_0$ . From (3.42.1.11) one obtains (as in (3.42.1.2)):

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

The governing equation (3.42.1.4) is the same:

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

The boundary term (3.42.1.5) becomes (see (3.42.1.12)):

$(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 (3.42.1.8), 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 (3.42.1.10) 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 (3.42.1.11).

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

Define the functional $L[u]$ in (3.42.1.7) and the space $U$ for the trial functions $u\in U$ in (3.42.1.8).
Define the weak form (3.42.1.11) and the two spaces $U$ and $U_0$ , where $u\in U$ and $\delta u \in U_0$ .
Define the strong form (3.42.1.13) 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 (3.42.1.11):

(3.42.1.14)¶ $\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 (3.42.1.14) 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 (3.42.1.14) 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:

(3.42.2.1)¶ $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:

(3.42.2.2)¶ $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:

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

Substituting (3.42.2.3) into (3.42.2.2) we obtain:

(3.42.2.4)¶ $[r^2 R' \delta R]_0^\infty = 0$

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

(3.42.2.5)¶ $\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 (3.42.2.5) vanishes. Or we allow $\delta R(\infty)$ to vary, and then (3.42.2.5) implies $R'(\infty) = 0$ .