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