Variational formulations with global basis functions¶

The finite element method is a very flexible approach for solving partial differential equations. Its two most attractive features are the ease of handling domains of complex shape in two and three dimensions and the variety of polynomials (with different properties and orders) that are available. The latter feature typically leads to errors proportional to $h^{d+1}$, where $h$ is the element length and $d$ is the polynomial degree. When the solution is sufficiently smooth, the ability to use larger $d$ creates methods that are much more computationally efficient than standard finite difference methods (and equally efficient finite difference methods are technically much harder to construct).

However, before we attack finite element methods, with localized basis functions, it can be easier from a pedagogical point of view to study approximations by global functions because the mathematics in this case gets simpler.

Basic principles for approximating differential equations¶

The finite element method is usually applied for discretization in space, and therefore spatial problems will be our focus in the coming sections. Extensions to time-dependent problems usually employs finite difference approximations in time.

The coming sections address at how global basis functions and the least squares, Galerkin, and collocation principles can be used to solve differential equations.

Differential equation models¶

Let us consider an abstract differential equation for a function $u(x)$ of one variable, written as

$$\begin{equation} \mathcal{L}(u) = 0,\quad x\in\Omega{\thinspace .} \label{_auto1} \tag{1} \end{equation}$$

Here are a few examples on possible choices of $\mathcal{L}(u)$, of increasing complexity:

$$\begin{equation} \mathcal{L}(u) = \frac{d^2u}{dx^2} - f(x), \label{fem:deq:1D:L1} \tag{2} \end{equation}$$

$$\begin{equation} \mathcal{L}(u) = \frac{d}{dx}\left({\alpha}(x)\frac{du}{dx}\right) + f(x), \label{fem:deq:1D:L2} \tag{3} \end{equation}$$

$$\begin{equation} \mathcal{L}(u) = \frac{d}{dx}\left({\alpha}(u)\frac{du}{dx}\right) - au + f(x), \label{fem:deq:1D:L3} \tag{4} \end{equation}$$

$$\begin{equation} \mathcal{L}(u) = \frac{d}{dx}\left({\alpha}(u)\frac{du}{dx}\right) + f(u,x) \label{fem:deq:1D:L4} \tag{5} {\thinspace .} \end{equation}$$

Both ${\alpha}(x)$ and $f(x)$ are considered as specified functions, while $a$ is a prescribed parameter. Differential equations corresponding to (2)-(3) arise in diffusion phenomena, such as stationary (time-independent) transport of heat in solids and flow of viscous fluids between flat plates. The form (4) arises when transient diffusion or wave phenomena are discretized in time by finite differences. The equation (5) appears in chemical models when diffusion of a substance is combined with chemical reactions. Also in biology, (5) plays an important role, both for spreading of species and in models involving generation and propagation of electrical signals.

Let $\Omega =[0,L]$ be the domain in one space dimension. In addition to the differential equation, $u$ must fulfill boundary conditions at the boundaries of the domain, $x=0$ and $x=L$. When $\mathcal{L}$ contains up to second-order derivatives, as in the examples above, we need one boundary condition at each of the (two) boundary points, here abstractly specified as

$$\begin{equation} \mathcal{B}_0(u)=0,\ x=0,\quad \mathcal{B}_1(u)=0,\ x=L \label{_auto2} \tag{6} \end{equation}$$

There are three common choices of boundary conditions:

$$\begin{equation} \mathcal{B}_i(u) = u - g,\quad \hbox{Dirichlet condition} \label{_auto3} \tag{7} \end{equation}$$

$$\begin{equation} \mathcal{B}_i(u) = -{\alpha} \frac{du}{dx} - g,\quad \hbox{Neumann condition} \label{_auto4} \tag{8} \end{equation}$$

$$\begin{equation} \mathcal{B}_i(u) = -{\alpha} \frac{du}{dx} - H(u-g),\quad \hbox{Robin condition} \label{_auto5} \tag{9} \end{equation}$$

Here, $g$ and $H$ are specified quantities.

From now on we shall use ${u_{\small\mbox{e}}}(x)$ as symbol for the exact solution, fulfilling

$$\begin{equation} \mathcal{L}({u_{\small\mbox{e}}})=0,\quad x\in\Omega, \label{_auto6} \tag{10} \end{equation}$$

while $u(x)$ is our notation for an approximate solution of the differential equation.

Remark on notation.

In the literature about the finite element method, it is common to use $u$ as the exact solution and $u_h$ as the approximate solution, where $h$ is a discretization parameter. However, the vast part of the present text is about the approximate solutions, and having a subscript $h$ attached all the time is cumbersome. Of equal importance is the close correspondence between implementation and mathematics that we strive to achieve in this text: when it is natural to use u and not u_h in code, we let the mathematical notation be dictated by the code's preferred notation. In the relatively few cases where we need to work with the exact solution of the PDE problem we call it ${u_{\small\mbox{e}}}$ in mathematics and u_e in the code (the function for computing u_e is named u_exact).

Simple model problems and their solutions¶

A common model problem used much in the forthcoming examples is

$$\begin{equation} -u''(x) = f(x),\quad x\in\Omega=[0,L],\quad u(0)=0,\ u(L)=D {\thinspace .} \label{fem:deq:1D:model1} \tag{11} \end{equation}$$

A closely related problem with a different boundary condition at $x=0$ reads

$$\begin{equation} -u''(x) = f(x),\quad x\in\Omega=[0,L],\quad u'(0)=C,\ u(L)=D{\thinspace .} \label{fem:deq:1D:model2} \tag{12} \end{equation}$$

A third variant has a variable coefficient,

$$\begin{equation} -({\alpha}(x)u'(x))' = f(x),\quad x\in\Omega=[0,L],\quad u'(0)=C,\ u(L)=D{\thinspace .} \label{fem:deq:1D:model3} \tag{13} \end{equation}$$

The solution $u$ to the model problem (11) can be determined as

$$\begin{align*} u'(x) &= -\int_0^x f(x) + c_0, \\ u(x) &= \int_0^x u'(x) + c_1, \end{align*}$$

where $c_0$ and $c_1$ are determined by the boundary conditions such that $u'(0) = C$ and $u(L) = D$.

Computing the solution is easily done using sympy. Some common code is defined first:

The following function computes the solution symbolically for the model problem (11):

Calling model1(2, L, D) results in the solution

$$\begin{equation} u(x) = \frac{1}{L}x \left(D + L^{2} - L x\right) \label{fem:deq:1D:model1:sol} \tag{14} \end{equation}$$

The model problem (12) can be solved by

to yield

$$\begin{equation} u(x) = - x^{2} + C x - C L + D + L^{2}, \label{fem:deq:1D:model2:sol} \tag{15} \end{equation}$$

if $f(x)=2$. Model (13) requires a bit more involved code,

With $f(x)=0$ and ${\alpha}(x)=1+x^2$ we get

$$u(x) = \frac{C \tan^{-1}\left (L \right ) - C \tan^{-1}\left (x \right ) + D \tan^{-1}\left (x \right )}{\tan^{-1}\left (L \right )}$$

Forming the residual¶

The fundamental idea is to seek an approximate solution $u$ in some space $V$,

$$V = \hbox{span}\{ {\psi}_0(x),\ldots,{\psi}_N(x)\},$$

which means that $u$ can always be expressed as a linear combination of the basis functions $\left\{ {{\psi}}_j \right\}_{j\in{\mathcal{I}_s}}$, with ${\mathcal{I}_s}$ as the index set $\{0,\ldots,N\}$:

$$u(x) = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x){\thinspace .}$$

The coefficients $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ are unknowns to be computed.

(Later, we will see that if we specify boundary values of $u$ different from zero, we must look for an approximate solution $u(x) = B(x) + \sum_{j} c_j{\psi}_j(x)$, where $\sum_{j}c_j{\psi}_j\in V$ and $B(x)$ is some function for incorporating the right boundary values. Because of $B(x)$, $u$ will not necessarily lie in $V$. This modification does not imply any difficulties.)

We need principles for deriving $N+1$ equations to determine the $N+1$ unknowns $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$. When approximating a given function $f$ by $u=\sum_jc_j{\varphi}_j$, a key idea is to minimize the square norm of the approximation error $e=u-f$ or (equivalently) demand that $e$ is orthogonal to $V$. Working with $e$ is not so useful here since the approximation error in our case is $e={u_{\small\mbox{e}}} - u$ and ${u_{\small\mbox{e}}}$ is unknown. The only general indicator we have on the quality of the approximate solution is to what degree $u$ fulfills the differential equation. Inserting $u=\sum_j c_j {\psi}_j$ into $\mathcal{L}(u)$ reveals that the result is not zero, because $u$ in general is an approximation and not identical to ${u_{\small\mbox{e}}}$. The nonzero result,

$$\begin{equation} R = \mathcal{L}(u) = \mathcal{L}(\sum_j c_j {\psi}_j), \label{_auto7} \tag{16} \end{equation}$$

is called the residual and measures the error in fulfilling the governing equation.

Various principles for determining $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ try to minimize $R$ in some sense. Note that $R$ varies with $x$ and the $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ parameters. We may write this dependence explicitly as

$$\begin{equation} R = R(x; c_0, \ldots, c_N){\thinspace .} \label{_auto8} \tag{17} \end{equation}$$

Below, we present three principles for making $R$ small: a least squares method, a projection or Galerkin method, and a collocation or interpolation method.

The least squares method¶

The least-squares method aims to find $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$ such that the square norm of the residual

$$\begin{equation} ||R|| = (R, R) = \int_{\Omega} R^2 {\, \mathrm{d}x} \label{_auto9} \tag{18} \end{equation}$$

is minimized. By introducing an inner product of two functions $f$ and $g$ on $\Omega$ as

$$\begin{equation} (f,g) = \int_{\Omega} f(x)g(x) {\, \mathrm{d}x}, \label{_auto10} \tag{19} \end{equation}$$

the least-squares method can be defined as

$$\begin{equation} \min_{c_0,\ldots,c_N} E = (R,R){\thinspace .} \label{_auto11} \tag{20} \end{equation}$$

Differentiating with respect to the free parameters $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$ gives the $N+1$ equations

$$\begin{equation} \int_{\Omega} 2R\frac{\partial R}{\partial c_i} {\, \mathrm{d}x} = 0\quad \Leftrightarrow\quad (R,\frac{\partial R}{\partial c_i})=0,\quad i\in{\mathcal{I}_s}{\thinspace .} \label{fem:deq:1D:LS:eq1} \tag{21} \end{equation}$$

The Galerkin method¶

The least-squares principle is equivalent to demanding the error to be orthogonal to the space $V$ when approximating a function $f$ by $u\in V$. With a differential equation we do not know the true error so we must instead require the residual $R$ to be orthogonal to $V$. This idea implies seeking $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$ such that

$$\begin{equation} (R,v)=0,\quad \forall v\in V{\thinspace .} \label{fem:deq:1D:Galerkin0} \tag{22} \end{equation}$$

This is the Galerkin method for differential equations.

The above abstract statement can be made concrete by choosing a concrete basis. For example, the statement is equivalent to $R$ being orthogonal to the $N+1$ basis functions $\{{\psi}_i\}$ spanning $V$ (and this is the most convenient way to express (22):

$$\begin{equation} (R,{\psi}_i)=0,\quad i\in{\mathcal{I}_s}, \label{fem:deq:1D:Galerkin} \tag{23} \end{equation}$$

resulting in $N+1$ equations for determining $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$.

The method of weighted residuals¶

A generalization of the Galerkin method is to demand that $R$ is orthogonal to some space $W$, but not necessarily the same space as $V$ where we seek the unknown function. This generalization is called the method of weighted residuals:

$$\begin{equation} (R,v)=0,\quad \forall v\in W{\thinspace .} \label{fem:deq:1D:WRM0} \tag{24} \end{equation}$$

If $\{w_0,\ldots,w_N\}$ is a basis for $W$, we can equivalently express the method of weighted residuals as

$$\begin{equation} (R,w_i)=0,\quad i\in{\mathcal{I}_s}{\thinspace .} \label{fem:deq:1D:WRM} \tag{25} \end{equation}$$

The result is $N+1$ equations for $\left\{ {c}_i \right\}_{i\in{\mathcal{I}_s}}$.

The least-squares method can also be viewed as a weighted residual method with $w_i = \partial R/\partial c_i$.

Variational formulation of the continuous problem.

Statements like (22), (23), (24), or (25)) are known as weak formulations or variational formulations. These equations are in this text primarily used for a numerical approximation $u\in V$, where $V$ is a finite-dimensional space with dimension $N+1$. However, we may also let the exact solution ${u_{\small\mbox{e}}}$ fulfill a variational formulation $(\mathcal{L}({u_{\small\mbox{e}}}),v)=0$ $\forall v\in V$, but the exact solution lies in general in a space with infinite dimensions (because an infinite number of parameters are needed to specify the solution). The variational formulation for ${u_{\small\mbox{e}}}$ in an infinite-dimensional space $V$ is a mathematical way of stating the problem and acts as an alternative to the usual (strong) formulation of a differential equation with initial and/or boundary conditions.

Much of the literature on finite element methods takes a differential equation problem and first transforms it to a variational formulation in an infinite-dimensional space $V$, before searching for an approximate solution in a finite-dimensional subspace of $V$. However, we prefer the more intuitive approach with an approximate solution $u$ in a finite-dimensional space $V$ inserted in the differential equation, and then the resulting residual is demanded to be orthogonal to $V$.

Remark on terminology.

The terms weak or variational formulations often refer to a statement like (22) or (24) after integration by parts has been performed (the integration by parts technique is explained in the section Integration by parts). The result after integration by parts is what is obtained after taking the first variation of a minimization problem (see the section Variational problems and minimization of functionals). However, in this text we use variational formulation as a common term for formulations which, in contrast to the differential equation $R=0$, instead demand that an average of $R$ is zero: $(R,v)=0$ for all $v$ in some space.

Test and trial functions¶

In the context of the Galerkin method and the method of weighted residuals it is common to use the name trial function for the approximate $u = \sum_j c_j {\psi}_j$.

The space containing the trial function is known as the trial space. The function $v$ entering the orthogonality requirement in the Galerkin method and the method of weighted residuals is called test function, and so are the ${\psi}_i$ or $w_i$ functions that are used as weights in the inner products with the residual. The space where the test functions comes from is naturally called the test space.

We see that in the method of weighted residuals the test and trial spaces are different and so are the test and trial functions. In the Galerkin method the test and trial spaces are the same (so far).

The collocation method¶

The idea of the collocation method is to demand that $R$ vanishes at $N+1$ selected points $x_{0},\ldots,x_{N}$ in $\Omega$:

$$\begin{equation} R(x_{i}; c_0,\ldots,c_N)=0,\quad i\in{\mathcal{I}_s}{\thinspace .} \label{fem:deq:1D:collocation} \tag{26} \end{equation}$$

The collocation method can also be viewed as a method of weighted residuals with Dirac delta functions as weighting functions. Let $\delta (x-x_{i})$ be the Dirac delta function centered around $x=x_{i}$ with the properties that $\delta (x-x_{i})=0$ for $x\neq x_{i}$ and

$$\begin{equation} \int_{\Omega} f(x)\delta (x-x_{i}) {\, \mathrm{d}x} = f(x_{i}),\quad x_{i}\in\Omega{\thinspace .} \label{fem:deq:1D:Dirac} \tag{27} \end{equation}$$

Intuitively, we may think of $\delta (x-x_{i})$ as a very peak-shaped function around $x=x_{i}$ with an integral $\int_{-\infty}^\infty \delta(x-x_{i})dx$ that evaluates to unity. Mathematically, it can be shown that $\delta (x-x_{i})$ is the limit of a Gaussian function centered at $x=x_{i}$ with a standard deviation that approaches zero. Using this latter model, we can roughly visualize delta functions as done in Figure. Because of (27), we can let $w_i=\delta(x-x_{i})$ be weighting functions in the method of weighted residuals, and (25) becomes equivalent to (26).

Approximation of delta functions by narrow Gaussian functions.

The subdomain collocation method¶

The idea of this approach is to demand the integral of $R$ to vanish over $N+1$ subdomains $\Omega_i$ of $\Omega$:

$$\begin{equation} \int_{\Omega_i} R\, {\, \mathrm{d}x}=0,\quad i\in{\mathcal{I}_s}{\thinspace .} \label{_auto12} \tag{28} \end{equation}$$

This statement can also be expressed as a weighted residual method

$$\begin{equation} \int_{\Omega} Rw_i\, {\, \mathrm{d}x}=0,\quad i\in{\mathcal{I}_s}, \label{_auto13} \tag{29} \end{equation}$$

where $w_i=1$ for $x\in\Omega_i$ and $w_i=0$ otherwise.

Examples on using the principles¶

Let us now apply global basis functions to illustrate the different principles for making the residual $R$ small.

The model problem¶

We consider the differential equation problem

$$\begin{equation} -u''(x) = f(x),\quad x\in\Omega=[0,L],\quad u(0)=0,\ u(L)=0 {\thinspace .} \label{fem:deq:1D:model1b} \tag{30} \end{equation}$$

Basis functions¶

Our choice of basis functions ${\psi}_i$ for $V$ is

$$\begin{equation} {\psi}_i(x) = {\sin\left((i+1)\pi\frac{x}{L}\right)},\quad i\in{\mathcal{I}_s}{\thinspace .} \label{fem:deq:1D:ex:sines:psi} \tag{31} \end{equation}$$

An important property of these functions is that ${\psi}_i(0)={\psi}_i(L)=0$, which means that the boundary conditions on $u$ are fulfilled:

$$u(0) = \sum_jc_j{\psi}_j(0) = 0,\quad u(L) = \sum_jc_j{\psi}_j(L) =0 {\thinspace .}$$

Another nice property is that the chosen sine functions are orthogonal on $\Omega$:

$$\begin{equation} \int\limits_0^L {\sin\left((i+1)\pi\frac{x}{L}\right)}{\sin\left((j+1)\pi\frac{x}{L}\right)}\, {\, \mathrm{d}x} = \left\lbrace \begin{array}{ll} \frac{1}{2} L & i=j \\ 0, & i\neq j \end{array}\right. \label{_auto14} \tag{32} \end{equation}$$

provided $i$ and $j$ are integers.

The residual¶

We can readily calculate the following explicit expression for the residual:

$$R(x;c_0, \ldots, c_N) = u''(x) + f(x),\nonumber$$

$$= \frac{d^2}{dx^2}\left(\sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x)\right) + f(x),\nonumber$$

$$\begin{equation} = \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j''(x) + f(x){\thinspace .} \label{fem:deq:1D:ex:sines:res} \tag{33} \end{equation}$$

The least squares method¶

The equations (21) in the least squares method require an expression for $\partial R/\partial c_i$. We have

$$\begin{equation} \frac{\partial R}{\partial c_i} = \frac{\partial}{\partial c_i} \left(\sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j''(x) + f(x)\right) = \sum_{j\in{\mathcal{I}_s}} \frac{\partial c_j}{\partial c_i}{\psi}_j''(x) = {\psi}_i''(x){\thinspace .} \label{_auto15} \tag{34} \end{equation}$$

The governing equations for the unknown parameters $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$ are then

$$\begin{equation} (\sum_j c_j {\psi}_j'' + f,{\psi}_i'')=0,\quad i\in{\mathcal{I}_s}, \label{_auto16} \tag{35} \end{equation}$$

which can be rearranged as

$$\begin{equation} \sum_{j\in{\mathcal{I}_s}}({\psi}_i'',{\psi}_j'')c_j = -(f,{\psi}_i''),\quad i\in{\mathcal{I}_s}{\thinspace .} \label{_auto17} \tag{36} \end{equation}$$

This is nothing but a linear system

$$\sum_{j\in{\mathcal{I}_s}}A_{i,j}c_j = b_i,\quad i\in{\mathcal{I}_s}{\thinspace .}$$

The entries in the coefficient matrix are given by

$$\begin{align*} A_{i,j} &= ({\psi}_i'',{\psi}_j'')\nonumber\\ & = \pi^4(i+1)^2(j+1)^2L^{-4}\int_0^L {\sin\left((i+1)\pi\frac{x}{L}\right)}{\sin\left((j+1)\pi\frac{x}{L}\right)}\, {\, \mathrm{d}x} \end{align*}$$

The orthogonality of the sine functions simplify the coefficient matrix:

$$\begin{equation} A_{i,j} = \left\lbrace \begin{array}{ll} {1\over2}L^{-3}\pi^4(i+1)^4 & i=j \\ 0, & i\neq j \end{array}\right. \label{_auto18} \tag{37} \end{equation}$$

The right-hand side reads

$$\begin{equation} b_i = -(f,{\psi}_i'') = (i+1)^2\pi^2L^{-2}\int_0^Lf(x){\sin\left((i+1)\pi\frac{x}{L}\right)}\, {\, \mathrm{d}x} \label{_auto19} \tag{38} \end{equation}$$

Since the coefficient matrix is diagonal we can easily solve for

$$\begin{equation} c_i = \frac{2L}{\pi^2(i+1)^2}\int_0^Lf(x){\sin\left((i+1)\pi\frac{x}{L}\right)}\, {\, \mathrm{d}x}{\thinspace .} \label{fem:deq:1D:ex:sines:solution} \tag{39} \end{equation}$$

With the special choice of $f(x)=2$, the coefficients can be calculated in sympy by

The answer becomes

$$c_i = 4 \frac{L^{2} \left(\left(-1\right)^{i} + 1\right)}{\pi^{3} \left(i^{3} + 3 i^{2} + 3 i + 1\right)}$$

Now, $1+(-1)^i=0$ for $i$ odd, so only the coefficients with even index are nonzero. Introducing $i=2k$ for $k=0,\ldots,N/2$ to count the relevant indices (for $N$ odd, $k$ goes to $(N-1)/2$), we get the solution

$$\begin{equation} u(x) = \sum_{k=0}^{N/2} \frac{8L^2}{\pi^3(2k+1)^3}{\sin\left((2k+1)\pi\frac{x}{L}\right)}{\thinspace .} \label{_auto20} \tag{40} \end{equation}$$

The coefficients decay very fast: $c_2 = c_0/27$, $c_4=c_0/125$. The solution will therefore be dominated by the first term,

$$u(x) \approx \frac{8L^2}{\pi^3}\sin\left(\pi\frac{x}{L}\right){\thinspace .}$$

The Galerkin method¶

The Galerkin principle (22) applied to (30) consists of inserting our special residual (33) in (22)

$$(u''+f,v)=0,\quad \forall v\in V,$$

or

$$\begin{equation} (u'',v) = -(f,v),\quad\forall v\in V{\thinspace .} \label{_auto21} \tag{41} \end{equation}$$

This is the variational formulation, based on the Galerkin principle, of our differential equation. The $\forall v\in V$ requirement is equivalent to demanding the equation $(u'',v) = -(f,v)$ to be fulfilled for all basis functions $v={\psi}_i$, $i\in{\mathcal{I}_s}$, see (22) and (23). We therefore have

$$\begin{equation} (\sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j'', {\psi}_i)=-(f,{\psi}_i),\quad i\in{\mathcal{I}_s}{\thinspace .} \label{_auto22} \tag{42} \end{equation}$$

This equation can be rearranged to a form that explicitly shows that we get a linear system for the unknowns $\left\{ {c}_j \right\}_{j\in{\mathcal{I}_s}}$:

$$\begin{equation} \sum_{j\in{\mathcal{I}_s}} ({\psi}_i,{\psi}_j'')c_j = (f, {\psi}_i),\quad i\in{\mathcal{I}_s}{\thinspace .} \label{_auto23} \tag{43} \end{equation}$$

For the particular choice of the basis functions (31) we get in fact the same linear system as in the least squares method because ${\psi}''= -(i+1)^2\pi^2L^{-2}{\psi}$. Consequently, the solution $u(x)$ becomes identical to the one produced by the least squares method.

The collocation method¶

For the collocation method (26) we need to decide upon a set of $N+1$ collocation points in $\Omega$. A simple choice is to use uniformly spaced points: $x_{i}=i\Delta x$, where $\Delta x = L/N$ in our case ($N\geq 1$). However, these points lead to at least two rows in the matrix consisting of zeros (since ${\psi}_i(x_{0})=0$ and ${\psi}_i(x_{N})=0$), thereby making the matrix singular and non-invertible. This forces us to choose some other collocation points, e.g., random points or points uniformly distributed in the interior of $\Omega$. Demanding the residual to vanish at these points leads, in our model problem (30), to the equations

$$\begin{equation} -\sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j''(x_{i}) = f(x_{i}),\quad i\in{\mathcal{I}_s}, \label{_auto24} \tag{44} \end{equation}$$

which is seen to be a linear system with entries

$$A_{i,j}=-{\psi}_j''(x_{i})= (j+1)^2\pi^2L^{-2}\sin\left((j+1)\pi \frac{x_i}{L}\right),$$

in the coefficient matrix and entries $b_i=2$ for the right-hand side (when $f(x)=2$).

The special case of $N=0$ can sometimes be of interest. A natural choice is then the midpoint $x_{0}=L/2$ of the domain, resulting in $A_{0,0} = -{\psi}_0''(x_{0}) = \pi^2L^{-2}$, $f(x_0)=2$, and hence $c_0=2L^2/\pi^2$.

Comparison¶

In the present model problem, with $f(x)=2$, the exact solution is $u(x)=x(L-x)$, while for $N=0$ the Galerkin and least squares method result in $u(x)=8L^2\pi^{-3}\sin (\pi x/L)$ and the collocation method leads to $u(x)=2L^2\pi^{-2}\sin (\pi x/L)$. We can quickly use sympy to verify that the maximum error occurs at the midpoint $x=L/2$ and find what the errors are. First we set up the error expressions:

If the derivative of the errors vanish at $x=L/2$, the errors reach their maximum values here (the errors vanish at the boundary points).

Finally, we can compute the maximum error at $x=L/2$ and evaluate the expressions numerically with three decimals:

The error in the collocation method is about 6 times larger than the error in the Galerkin or least squares method.

Integration by parts¶

A problem arises if we want to apply popular finite element functions to solve our model problem (30) by the standard least squares, Galerkin, or collocation methods: the piecewise polynomials ${\psi}_i(x)$ have discontinuous derivatives at the cell boundaries which makes it problematic to compute the second-order derivative. This fact actually makes the least squares and collocation methods less suitable for finite element approximation of the unknown function. (By rewriting the equation $-u''=f$ as a system of two first-order equations, $u'=v$ and $-v'=f$, the least squares method can be applied. Also, differentiating discontinuous functions can actually be handled by distribution theory in mathematics.) The Galerkin method and the method of weighted residuals can, however, be applied together with finite element basis functions if we use integration by parts as a means for transforming a second-order derivative to a first-order one.

Consider the model problem (30) and its Galerkin formulation

$$-(u'',v) = (f,v)\quad\forall v\in V{\thinspace .}$$

Using integration by parts in the Galerkin method, we can "move" a derivative of $u$ onto $v$:

$$\int_0^L u''(x)v(x) {\, \mathrm{d}x} = - \int_0^Lu'(x)v'(x){\, \mathrm{d}x} + [vu']_0^L\nonumber$$

$$\begin{equation} = - \int_0^Lu'(x)v'(x) {\, \mathrm{d}x} + u'(L)v(L) - u'(0)v(0){\thinspace .} \label{fem:deq:1D:intbyparts} \tag{45} \end{equation}$$

Usually, one integrates the problem at the stage where the $u$ and $v$ functions enter the formulation. Alternatively, but less common, we can integrate by parts in the expressions for the matrix entries:

$$\int_0^L{\psi}_i(x){\psi}_j''(x) {\, \mathrm{d}x} = - \int_0^L{\psi}_i'(x){\psi}_j'(x) dx + [{\psi}_i{\psi}_j']_0^L\nonumber$$

$$\begin{equation} = - \int_0^L{\psi}_i'(x){\psi}_j'(x) {\, \mathrm{d}x} + {\psi}_i(L){\psi}_j'(L) - {\psi}_i(0){\psi}_j'(0){\thinspace .} \label{fem:deq:1D:intbyparts0} \tag{46} \end{equation}$$

Integration by parts serves to reduce the order of the derivatives and to make the coefficient matrix symmetric since $({\psi}_i',{\psi}_j') = ({\psi}_j',{\psi}_i')$. The symmetry property depends on the type of terms that enter the differential equation. As will be seen later, integration by parts also provides a method for implementing boundary conditions involving $u'$.

With the choice (31) of basis functions we see that the "boundary terms" ${\psi}_i(L){\psi}_j'(L)$ and ${\psi}_i(0){\psi}_j'(0)$ vanish since ${\psi}_i(0)={\psi}_i(L)=0$.

We therefore end up with the following alternative Galerkin formulation:

$$-(u'',v) = (u', v') = (f,v)\quad \forall v\in V{\thinspace .}$$

Weak form¶

Since the variational formulation after integration by parts make weaker demands on the differentiability of $u$ and the basis functions ${\psi}_i$, the resulting integral formulation is referred to as a weak form of the differential equation problem. The original variational formulation with second-order derivatives, or the differential equation problem with second-order derivative, is then the strong form, with stronger requirements on the differentiability of the functions.

For differential equations with second-order derivatives, expressed as variational formulations and solved by finite element methods, we will always perform integration by parts to arrive at expressions involving only first-order derivatives.

Boundary function¶

So far we have assumed zero Dirichlet boundary conditions, typically $u(0)=u(L)=0$, and we have demanded that ${\psi}_i(0)={\psi}_i(L)=0$ for $i\in{\mathcal{I}_s}$. What about a boundary condition like $u(L)=D\neq0$? This condition immediately faces a problem: $u = \sum_j c_j{\varphi}_j(L) = 0$ since all ${\varphi}_i(L)=0$.

We remark that we faced exactly the same problem where we considered Fourier series approximations of functions that where non-zero at the boundaries. We will use the same trick as we did earlier to get around this problem.

A boundary condition of the form $u(L)=D$ can be implemented by demanding that all ${\psi}_i(L)=0$, but adding a boundary function $B(x)$ with the right boundary value, $B(L)=D$, to the expansion for $u$:

$$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x) {\thinspace .}$$

This $u$ gets the right value at $x=L$:

$$u(L) = B(L) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(L) = B(L) = D{\thinspace .}$$

The idea is that for any boundary where $u$ is known we demand ${\psi}_i$ to vanish and construct a function $B(x)$ to attain the boundary value of $u$. There are no restrictions on how $B(x)$ varies with $x$ in the interior of the domain, so this variation needs to be constructed in some way. Exactly how we decide the variation to be, is not important.

For example, with $u(0)=0$ and $u(L)=D$, we can choose $B(x)=x D/L$, since this form ensures that $B(x)$ fulfills the boundary conditions: $B(0)=0$ and $B(L)=D$. The unknown function is then sought on the form

$$\begin{equation} u(x) = \frac{x}{L}D + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x), \label{fem:deq:1D:essBC:Bfunc:u1} \tag{47} \end{equation}$$

with ${\psi}_i(0)={\psi}_i(L)=0$.

The particular shape of the $B(x)$ function is not important as long as its boundary values are correct. For example, $B(x)=D(x/L)^p$ for any power $p$ will work fine in the above example. Another choice could be $B(x)=D\sin (\pi x/(2L))$.

As a more general example, consider a domain $\Omega = [a,b]$ where the boundary conditions are $u(a)=U_a$ and $u(b)=U_b$. A class of possible $B(x)$ functions is

$$\begin{equation} B(x)=U_a + \frac{U_b-U_a}{(b-a)^p}(x-a)^p,\quad p>0 {\thinspace .} \label{fem:deq:1D:essBC:Bfunc:gen} \tag{48} \end{equation}$$

Real applications will most likely use the simplest version, $p=1$, but here such a $p$ parameter was included to demonstrate that there are many choices of $B(x)$ in a problem. Fortunately, there is a general, unique technique for constructing $B(x)$ when we use finite element basis functions for $V$.

How to deal with nonzero Dirichlet conditions.

The general procedure of incorporating Dirichlet boundary conditions goes as follows. Let $\partial\Omega_E$ be the part(s) of the boundary $\partial\Omega$ of the domain $\Omega$ where $u$ is specified. Set ${\psi}_i=0$ at the points in $\partial\Omega_E$ and seek $u$ as

$$\begin{equation} u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x), \label{fem:deq:1D:essBC:Bfunc:u2} \tag{49} \end{equation}$$

where $B(x)$ equals the boundary conditions on $u$ at $\partial\Omega_E$.

Remark. With the $B(x)$ term, $u$ does not in general lie in $V=\hbox{span}\, \{{\psi}_0,\ldots,{\psi}_N\}$ anymore. Moreover, when a prescribed value of $u$ at the boundary, say $u(a)=U_a$ is different from zero, it does not make sense to say that $u$ lies in a vector space, because this space does not obey the requirements of addition and scalar multiplication. For example, $2u$ does not lie in the space since its boundary value is $2U_a$, which is incorrect. It only makes sense to split $u$ in two parts, as done above, and have the unknown part $\sum_j c_j {\psi}_j$ in a proper function space.

Computing with global polynomials¶

The next example uses global polynomials and shows that if our solution, modulo boundary conditions, lies in the space spanned by these polynomials, then the Galerkin method recovers the exact solution.

Computing with Dirichlet and Neumann conditions¶

Let us perform the necessary calculations to solve

$$-u''(x)=2,\quad x\in \Omega=[0,1],\quad u'(0)=C,\ u(1)=D,$$

using a global polynomial basis ${\psi}_i\sim x^i$. The requirements on ${\psi}_i$ is that ${\psi}_i(1)=0$, because $u$ is specified at $x=1$, so a proper set of polynomial basis functions can be

$${\psi}_i(x)=(1-x)^{i+1}, \quad i\in{\mathcal{I}_s}{\thinspace .}$$

A suitable $B(x)$ function to handle the boundary condition $u(1)=D$ is $B(x)=Dx$. The variational formulation becomes

$$(u',v') = (2,v) - Cv(0)\quad\forall v\in V{\thinspace .}$$

From inserting $u=B + \sum_{j}c_j{\psi}_j$ and choosing $v={\psi}_i$ we get

$$\sum_{j\in{\mathcal{I}_s}} ({\psi}_j',{\psi}_i')c_j = (2,{\psi}_i) - (B',{\psi}_i') - C{\psi}_i(0),\quad i\in{\mathcal{I}_s}{\thinspace .}$$

The entries in the linear system are then

$$\begin{align*} A_{i,j} &= ({\psi}_j',{\psi}_i') = \int_{0}^1 {\psi}_i'(x){\psi}_j'(x){\, \mathrm{d}x} = \int_0^1 (i+1)(j+1)(1-x)^{i+j}{\, \mathrm{d}x}\\ &= \frac{(i+1)(j+1)}{i + j + 1},\\ b_i &= (2,{\psi}_i) - (D,{\psi}_i') -C{\psi}_i(0)\\ &= \int_0^1\left( 2{\psi}_i(x) - D{\psi}_i'(x)\right){\, \mathrm{d}x} -C{\psi}_i(0)\\ &= \int_0^1 \left( 2(1-x)^{i+1} + D(i+1)(1-x)^i\right){\, \mathrm{d}x} -C\\ &= \frac{(D-C)(i+2) + 2}{i+2} = D - C + \frac{2}{i+2} {\thinspace .} \end{align*}$$

Relevant sympy commands to help calculate these expressions are

The output becomes

    A_ij: (i + 1)*(j + 1)/(i + j + 1)
    b_i: ((-C + D)*(i + 2) + 2)/(i + 2)

We can now choose some $N$ and form the linear system, say for $N=1$:

The system becomes

$$\left(\begin{array}{cc} 1 & 1\\ 1 & 4/3 \end{array}\right) \left(\begin{array}{c} c_0\\ c_1 \end{array}\right) = \left(\begin{array}{c} 1-C+D\\ 2/3 -C + D \end{array}\right)$$

The solution (c = A.LUsolve(b)) becomes $c_0=2 -C+D$ and $c_1=-1$, resulting in

$$\begin{equation} u(x) = 1 -x^2 + D + C(x-1), \label{_auto25} \tag{50} \end{equation}$$

We can form this $u$ in sympy and check that the differential equation and the boundary conditions are satisfied:

The output becomes

    u: C*x - C + D - x**2 + 1
    u'': -2
    BC x=0: C
    BC x=1: D

The complete sympy code is found in u_xx_2_CD.py.

The exact solution is found by integrating twice and applying the boundary conditions, either by hand or using sympy as shown in the section Simple model problems and their solutions. It appears that the numerical solution coincides with the exact one.

Abstract notation for variational formulations¶

We have seen that variational formulations end up with a formula involving $u$ and $v$, such as $(u',v')$ and a formula involving $v$ and known functions, such as $(f,v)$. A widely used notation is to introduce an abstract variational statement written as

$$a(u,v)=L(v)\quad\forall v\in V,$$

where $a(u,v)$ is a so-called bilinear form involving all the terms that contain both the test and trial function, while $L(v)$ is a linear form containing all the terms without the trial function. For example, the statement

$$\int_{\Omega} u' v' {\, \mathrm{d}x} = \int_{\Omega} fv{\, \mathrm{d}x}\quad\hbox{or}\quad (u',v') = (f,v) \quad\forall v\in V$$

can be written in abstract form: find $u$ such that

$$a(u,v) = L(v)\quad \forall v\in V,$$

where we have the definitions

$$a(u,v) = (u',v'),\quad L(v) = (f,v){\thinspace .}$$

The term linear means that

$$L(\alpha_1 v_1 + \alpha_2 v_2) =\alpha_1 L(v_1) + \alpha_2 L(v_2)$$

for two test functions $v_1$ and $v_2$, and scalar parameters $\alpha_1$ and $\alpha_2$. Similarly, the term bilinear means that $a(u,v)$ is linear in both its arguments:

$$\begin{align*} a(\alpha_1 u_1 + \alpha_2 u_2, v) &= \alpha_1 a(u_1,v) + \alpha_2 a(u_2, v), \\ a(u, \alpha_1 v_1 + \alpha_2 v_2) &= \alpha_1 a(u,v_1) + \alpha_2 a(u, v_2) {\thinspace .} \end{align*}$$

In nonlinear problems these linearity properties do not hold in general and the abstract notation is then

$$F(u;v)=0\quad\forall v\in V{\thinspace .}$$

The matrix system associated with $a(u,v)=L(v)$ can also be written in an abstract form by inserting $v={\psi}_i$ and $u=\sum_j c_j{\psi}_j$ in $a(u,v)=L(v)$. Using the linear properties, we get

$$\sum_{j\in{\mathcal{I}_s}} a({\psi}_j,{\psi}_i) c_j = L({\psi}_i),\quad i\in{\mathcal{I}_s},$$

which is a linear system

$$\sum_{j\in{\mathcal{I}_s}}A_{i,j}c_j = b_i,\quad i\in{\mathcal{I}_s},$$

where

$$A_{i,j} =a({\psi}_j,{\psi}_i), \quad b_i = L({\psi}_i){\thinspace .}$$

In many problems, $a(u,v)$ is symmetric such that $a({\psi}_j,{\psi}_i) = a({\psi}_i,{\psi}_j)$. In those cases the coefficient matrix becomes symmetric, $A_{i,j}=A_{j,i}$, a property that can simplify solution algorithms for linear systems and make them more stable. The property also reduces memory requirements and the computational work.

The abstract notation $a(u,v)=L(v)$ for linear differential equation problems is much used in the literature and in description of finite element software (in particular the FEniCS documentation). We shall frequently summarize variational forms using this notation.

Examples on variational formulations¶

The following sections derive variational formulations for some prototype differential equations in 1D, and demonstrate how we with ease can handle variable coefficients, mixed Dirichlet and Neumann boundary conditions, first-order derivatives, and nonlinearities.

Variable coefficient¶

Consider the problem

$$\begin{equation} -\frac{d}{dx}\left( {\alpha}(x)\frac{du}{dx}\right) = f(x),\quad x\in\Omega =[0,L],\ u(0)=C,\ u(L)=D{\thinspace .} \label{_auto26} \tag{51} \end{equation}$$

There are two new features of this problem compared with previous examples: a variable coefficient ${\alpha} (x)$ and nonzero Dirichlet conditions at both boundary points.

Let us first deal with the boundary conditions. We seek

$$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_i(x){\thinspace .}$$

Since the Dirichlet conditions demand

$${\psi}_i(0)={\psi}_i(L)=0,\quad i\in{\mathcal{I}_s},$$

the function $B(x)$ must fulfill $B(0)=C$ and $B(L)=D$. The we are guaranteed that $u(0)=C$ and $u(L)=D$. How $B$ varies in between $x=0$ and $x=L$ is not of importance. One possible choice is

$$B(x) = C + \frac{1}{L}(D-C)x,$$

which follows from (48) with $p=1$.

We seek $(u-B)\in V$. As usual,

$$V = \hbox{span}\{{\psi}_0,\ldots,{\psi}_N\}{\thinspace .}$$

Note that any $v\in V$ has the property $v(0)=v(L)=0$.

The residual arises by inserting our $u$ in the differential equation:

$$R = -\frac{d}{dx}\left( {\alpha}\frac{du}{dx}\right) -f{\thinspace .}$$

Galerkin's method is

$$(R, v) = 0,\quad \forall v\in V,$$

or written with explicit integrals,

$$\int_{\Omega} \left(-\frac{d}{dx}\left( {\alpha}\frac{du}{dx}\right) -f\right)v {\, \mathrm{d}x} = 0,\quad \forall v\in V {\thinspace .}$$

We proceed with integration by parts to lower the derivative from second to first order:

$$-\int_{\Omega} \frac{d}{dx}\left( {\alpha}(x)\frac{du}{dx}\right) v {\, \mathrm{d}x} = \int_{\Omega} {\alpha}(x)\frac{du}{dx}\frac{dv}{dx}{\, \mathrm{d}x} - \left[{\alpha}\frac{du}{dx}v\right]_0^L {\thinspace .}$$

The boundary term vanishes since $v(0)=v(L)=0$. The variational formulation is then

$$\int_{\Omega} {\alpha}(x)\frac{du}{dx}\frac{dv}{dx}{\, \mathrm{d}x} = \int_{\Omega} f(x)v{\, \mathrm{d}x},\quad \forall v\in V{\thinspace .}$$

The variational formulation can alternatively be written in a more compact form:

$$({\alpha} u',v') = (f,v),\quad \forall v\in V {\thinspace .}$$

The corresponding abstract notation reads

$$a(u,v)=L(v)\quad\forall v\in V,$$

with

$$a(u,v)= ({\alpha} u',v'),\quad L(v)=(f,v) {\thinspace .}$$

We may insert $u=B + \sum_jc_j{\psi}_j$ and $v={\psi}_i$ to derive the linear system:

$$({\alpha} B' + {\alpha} \sum_{j\in{\mathcal{I}_s}} c_j {\psi}_j', {\psi}_i') = (f,{\psi}_i), \quad i\in{\mathcal{I}_s} {\thinspace .}$$

Isolating everything with the $c_j$ coefficients on the left-hand side and all known terms on the right-hand side gives

$$\sum_{j\in{\mathcal{I}_s}} ({\alpha}{\psi}_j', {\psi}_i')c_j = (f,{\psi}_i) + (\alpha (D-C)L^{-1}, {\psi}_i'), \quad i\in{\mathcal{I}_s} {\thinspace .}$$

This is nothing but a linear system $\sum_j A_{i,j}c_j=b_i$ with

$$\begin{align*} A_{i,j} &= (\alpha {\psi}_j', {\psi}_i') = \int_{\Omega} {\alpha}(x){\psi}_j'(x), {\psi}_i'(x){\, \mathrm{d}x},\\ b_i &= (f,{\psi}_i) + (\alpha (D-C)L^{-1},{\psi}_i')= \int_{\Omega} \left(f(x){\psi}_i(x) + {\alpha}(x)\frac{D-C}{L}{\psi}_i'(x)\right) {\, \mathrm{d}x} {\thinspace .} \end{align*}$$

First-order derivative in the equation and boundary condition¶

The next problem to formulate in terms of a variational form reads

$$\begin{equation} -u''(x) + bu'(x) = f(x),\quad x\in\Omega =[0,L],\ u(0)=C,\ u'(L)=E{\thinspace .} \label{_auto27} \tag{52} \end{equation}$$

The new features are a first-order derivative $u'$ in the equation and the boundary condition involving the derivative: $u'(L)=E$. Since we have a Dirichlet condition at $x=0$, we must force ${\psi}_i(0)=0$ and use a boundary function to take care of the condition $u(0)=C$. Because there is no Dirichlet condition on $x=L$ we do not make any requirements to ${\psi}_i(L)$. The simplest possible choice of $B(x)$ is $B(x)=C$.

The expansion for $u$ becomes

$$u = C + \sum_{j\in{\mathcal{I}_s}} c_j {\psi}_i(x) {\thinspace .}$$

The variational formulation arises from multiplying the equation by a test function $v\in V$ and integrating over $\Omega$:

$$(-u'' + bu' - f, v) = 0,\quad\forall v\in V$$

We apply integration by parts to the $u''v$ term only. Although we could also integrate $u' v$ by parts, this is not common. The result becomes

$$(u',v') + (bu',v) = (f,v) + [u' v]_0^L, \quad\forall v\in V {\thinspace .}$$

Now, $v(0)=0$ so

$$[u' v]_0^L = u'(L)v(L) = E v(L),$$

because $u'(L)=E$. Thus, integration by parts allows us to take care of the Neumann condition in the boundary term.

Natural and essential boundary conditions.

A common mistake is to forget a boundary term like $[u'v]_0^L$ in the integration by parts. Such a mistake implies that we actually impose the condition $u'=0$ unless there is a Dirichlet condition (i.e., $v=0$) at that point! This fact has great practical consequences, because it is easy to forget the boundary term, and that implicitly set a boundary condition!

Since homogeneous Neumann conditions can be incorporated without "doing anything" (i.e., omitting the boundary term), and non-homogeneous Neumann conditions can just be inserted in the boundary term, such conditions are known as natural boundary conditions. Dirichlet conditions require more essential steps in the mathematical formulation, such as forcing all ${\varphi}_i=0$ on the boundary and constructing a $B(x)$, and are therefore known as essential boundary conditions.

The final variational form reads

$$(u',v') + (bu',v) = (f,v) + E v(L), \quad\forall v\in V {\thinspace .}$$

In the abstract notation we have

$$a(u,v)=L(v)\quad\forall v\in V,$$

with the particular formulas

$$a(u,v)=(u',v') + (bu',v),\quad L(v)= (f,v) + E v(L){\thinspace .}$$

The associated linear system is derived by inserting $u=B+\sum_jc_j{\psi}_j$ and replacing $v$ by ${\psi}_i$ for $i\in{\mathcal{I}_s}$. Some algebra results in

$$\sum_{j\in{\mathcal{I}_s}} \underbrace{(({\psi}_j',{\psi}_i') + (b{\psi}_j',{\psi}_i))}_{A_{i,j}} c_j = \underbrace{(f,{\psi}_i) + E {\psi}_i(L)}_{b_i} {\thinspace .}$$

Observe that in this problem, the coefficient matrix is not symmetric, because of the term

$$(b{\psi}_j',{\psi}_i)=\int_{\Omega} b{\psi}_j'{\psi}_i {\, \mathrm{d}x} \neq \int_{\Omega} b {\psi}_i' {\psi}_j {\, \mathrm{d}x} = ({\psi}_i',b{\psi}_j) {\thinspace .}$$

Nonlinear coefficient¶

Finally, we show that the techniques used above to derive variational forms apply to nonlinear differential equation problems as well. Here is a model problem with a nonlinear coefficient $\alpha(u)$ and a nonlinear right-hand side $f(u)$:

$$\begin{equation} -({\alpha}(u)u')' = f(u),\quad x\in [0,L],\ u(0)=0,\ u'(L)=E {\thinspace .} \label{_auto28} \tag{53} \end{equation}$$

Our space $V$ has basis $\left\{ {{\psi}}_i \right\}_{i\in{\mathcal{I}_s}}$, and because of the condition $u(0)=0$, we must require ${\psi}_i(0)=0$, $i\in{\mathcal{I}_s}$.

Galerkin's method is about inserting the approximate $u$, multiplying the differential equation by $v\in V$, and integrate,

$$-\int_0^L \frac{d}{dx}\left({\alpha}(u)\frac{du}{dx}\right)v {\, \mathrm{d}x} = \int_0^L f(u)v {\, \mathrm{d}x}\quad\forall v\in V {\thinspace .}$$

The integration by parts does not differ from the case where we have ${\alpha}(x)$ instead of ${\alpha}(u)$:

$$\int_0^L {\alpha}(u)\frac{du}{dx}\frac{dv}{dx}{\, \mathrm{d}x} = \int_0^L f(u)v{\, \mathrm{d}x} + [{\alpha}(u)vu']_0^L\quad\forall v\in V {\thinspace .}$$

The term ${\alpha}(u(0))v(0)u'(0)=0$ since $v(0)$. The other term, ${\alpha}(u(L))v(L)u'(L)$, is used to impose the other boundary condition $u'(L)=E$, resulting in

$$\int_0^L {\alpha}(u)\frac{du}{dx}\frac{dv}{dx}{\, \mathrm{d}x} = \int_0^L f(u)v{\, \mathrm{d}x} + {\alpha}(u(L))v(L)E\quad\forall v\in V,$$

or alternatively written more compactly as

$$({\alpha}(u)u', v') = (f(u),v) + {\alpha}(u(L))v(L)E\quad\forall v\in V {\thinspace .}$$

Since the problem is nonlinear, we cannot identify a bilinear form $a(u,v)$ and a linear form $L(v)$. An abstract formulation is typically find $u$ such that

$$F(u;v) = 0\quad\forall v\in V,$$

with

$$F(u;v) = (a(u)u', v') - (f(u),v) - a(L)v(L)E {\thinspace .}$$

By inserting $u=\sum_j c_j{\psi}_j$ and $v={\psi}_i$ in $F(u;v)$, we get a nonlinear system of algebraic equations for the unknowns $c_i$, $i\in{\mathcal{I}_s}$. Such systems must be solved by constructing a sequence of linear systems whose solutions hopefully converge to the solution of the nonlinear system. Frequently applied methods are Picard iteration and Newton's method.

Implementation of the algorithms¶

Our hand calculations can benefit greatly by symbolic computing, as shown earlier, so it is natural to extend our approximation programs based on sympy to the problem domain of variational formulations.

Extensions of the code for approximation¶

The user must prepare a function integrand_lhs(psi, i, j) for returning the integrand of the integral that contributes to matrix entry $(i,j)$ on the left-hand side. The psi variable is a Python dictionary holding the basis functions and their derivatives in symbolic form. More precisely, psi[q] is a list of

$$\{\frac{d^q{\psi}_0}{dx^q},\ldots,\frac{d^q{\psi}_{N_n-1}}{dx^q}\} {\thinspace .}$$

Similarly, integrand_rhs(psi, i) returns the integrand for entry number $i$ in the right-hand side vector.

Since we also have contributions to the right-hand side vector (and potentially also the matrix) from boundary terms without any integral, we introduce two additional functions, boundary_lhs(psi, i, j) and boundary_rhs(psi, i) for returning terms in the variational formulation that are not to be integrated over the domain $\Omega$. Examples, to be shown later, will explain in more detail how these user-supplied functions may look like.

The linear system can be computed and solved symbolically by the following function:

Fallback to numerical methods¶

Not surprisingly, symbolic solution of differential equations, discretized by a Galerkin or least squares method with global basis functions, is of limited interest beyond the simplest problems, because symbolic integration might be very time consuming or impossible, not only in sympy but also in WolframAlpha (which applies the perhaps most powerful symbolic integration software available today: Mathematica). Numerical integration as an option is therefore desirable.

The extended solver function below tries to combine symbolic and numerical integration. The latter can be enforced by the user, or it can be invoked after a non-successful symbolic integration (being detected by an Integral object as the result of the integration in sympy).

Note that for a numerical integration, symbolic expressions must be converted to Python functions (using lambdify), and the expressions cannot contain other symbols than x. The real solver routine in the varform1D.py file has error checking and meaningful error messages in such cases. The solver code below is a condensed version of the real one, with the purpose of showing how to automate the Galerkin or least squares method for solving differential equations in 1D with global basis functions:

Example with constant right-hand side¶

To demonstrate the code above, we address

$$-u''(x)=b,\quad x\in\Omega=[0,1],\quad u(0)=1,\ u(1)=0,$$

with $b$ as a (symbolic) constant. A possible basis for the space $V$ is ${\psi}_i(x) = x^{i+1}(1-x)$, $i\in{\mathcal{I}_s}$. Note that ${\psi}_i(0)={\psi}_i(1)=0$ as required by the Dirichlet conditions. We need a $B(x)$ function to take care of the known boundary values of $u$. Any function $B(x)=1-x^p$, $p\in\mathbb{R}$, is a candidate, and one arbitrary choice from this family is $B(x)=1-x^3$. The unknown function is then written as

$$u(x) = B(x) + \sum_{j\in{\mathcal{I}_s}} c_j{\psi}_j(x){\thinspace .}$$

Let us use the Galerkin method to derive the variational formulation. Multiplying the differential equation by $v$ and integrating by parts yield

$$\int_0^1 u'v' {\, \mathrm{d}x} = \int_0^1 fv {\, \mathrm{d}x}\quad\forall v\in V,$$

and with $u=B + \sum_jc_j{\psi}_j$ we get the linear system

$$\begin{equation} \sum_{j\in{\mathcal{I}_s}}\left(\int_0^1{\psi}_i'{\psi}_j' {\, \mathrm{d}x}\right)c_j = \int_0^1(f{\psi}_i-B'{\psi}_i') {\, \mathrm{d}x}, \quad i\in{\mathcal{I}_s}{\thinspace .} \label{_auto29} \tag{54} \end{equation}$$

The application can be coded as follows with sympy:

The printout of u reads -b*x**2/2 + b*x/2 - x + 1. Note that expanding u, before simplifying, is necessary in the present case to get a compact, final expression with sympy. Doing expand before simplify is a common strategy for simplifying expressions in sympy. However, a non-expanded u might be preferable in other cases - this depends on the problem in question.

The exact solution ${u_{\small\mbox{e}}}(x)$ can be derived by some sympy code that closely follows the examples in the section Simple model problems and their solutions. The idea is to integrate $-u''=b$ twice and determine the integration constants from the boundary conditions:

The last line prints 0, which is not surprising when ${u_{\small\mbox{e}}}(x)$ is a parabola and our approximate $u$ contains polynomials up to degree 4. It suffices to have $N=1$, i.e., polynomials of degree 2, to recover the exact solution.

We can play around with the code and test that with $f=Kx^p$, for some constants $K$ and $p$, the solution is a polynomial of degree $p+2$, and $N=p+1$ guarantees that the approximate solution is exact.

Although the symbolic code is capable of integrating many choices of $f(x)$, the symbolic expressions for $u$ quickly become lengthy and non-informative, so numerical integration in the code, and hence numerical answers, have the greatest application potential.