I've been reading some resources on the web about Galerkin methods to solve PDEs, but I'm not clear about something. The following is my own account of what I have understood.

Consider the following boundary value problem (BVP):

$$L[u(x,y)]=0 \quad \text{on}\quad (x,y)\in\Omega, \qquad S[u]=0 \quad \text{on} \quad (x,y)\in\partial\Omega$$

where $L$ is a 2nd order linear differentiation operator, $\Omega\subset\mathbb{R}^2$ is the domain of the BVP, $\partial\Omega$ is the boundary of the domain, and $S$ is a 1st order linear differential operator. Expess $u(x,y)$ as an aproximation of the form:

$$u(x,y)\approx \sum_{i=1}^N a_i g_i(x,y)$$

where the $g_i$ are a set of functions that we will use to approximate $u$. Substituting in the BVP:

$$\sum_i a_i L[g_i(x,y)]=R(a_1,...,a_N,x,y)$$

Since our approximation is not exact, the residual $R$ is not exactly zero. In the Galerkin-Ritz-Raleigh method we minimize $R$ with respect to the set of approximating functions by requiring $\langle R,g_i \rangle = 0$. Hence

$$\langle R,g_i \rangle = \sum_{j=1}^N a_j \langle L[g_j],g_i \rangle = 0$$

Therefore, to find the coefficients $a_i$, we must solve the matrix equation:

$$\left( \begin{array}{ccc} \left\langle L\left[g_1\right],g_1\right\rangle & \ldots & \left\langle L\left[g_N\right],g_1\right\rangle \\ \ldots & \ldots & \ldots \\ \left\langle L\left[g_1\right],g_N\right\rangle & \ldots & \left\langle L\left[g_N\right],g_N\right\rangle \end{array} \right)\left( \begin{array}{c} a_1 \\ \ldots \\ a_N \end{array} \right)=0$$

My question is: How do I incorporate the boundary conditions into this?

EDIT: Originally the question said that $S[u]$ was a 2nd order linear differential operator. I changed it to a 1st order linear differential operator.

  • 1
    $\begingroup$ becko, welcome to scicomp! Our policy on cross-posting follows that of other Stack Exchange sites. It is permissible to cross-post if you tailor the same question (more or less) to different audiences. It is permissible to ask your question to be migrated to another site after some time, if you feel that your question is not getting answered satisfactorily (or at all) on the site where it is initially posted. $\endgroup$ Commented Feb 18, 2012 at 2:27
  • $\begingroup$ However, it is generally considered abusive behavior to cross-post. If you look at the list of beta sites on Area 51, many of them are still in public beta after a year. We're still going to be around for a while (at least on a longer time scale than it takes most questions on this site to be answered). Furthermore, unless the math users that answer your question are also scicomp users, they won't get proper credit or attribution for their answer here on scicomp if you copy-paste it from math, and vice versa. $\endgroup$ Commented Feb 18, 2012 at 2:33
  • 1
    $\begingroup$ Are you sure that $S$ is also a second order operator? In general, that is not a well-posed problem. For example, if $S = L$, then you are asking for $L u = 0$ in $\bar\Omega$ which has horribly non-unique solutions (e.g. the solution of the PDE on any larger domain with any boundary conditions is also a solution). Usually we ask for $S$ to be a (possibly nonlinear) first order operator. $\endgroup$
    – Jed Brown
    Commented Feb 18, 2012 at 6:49
  • 2
    $\begingroup$ Even if $S \ne L$, you're still looking at non-unique solutions. Consider if $L$ is the laplace operator and $S$ is any other second order linear operator. Then any $u$ such that $\nabla u=\mathbf{k}$ for some constant vector $\mathbf{k}$ can be added to a solution to make another solution. $\endgroup$
    – Dan
    Commented Feb 18, 2012 at 18:43
  • 1
    $\begingroup$ @GeoffOxberry For your peace of mind, the duplicate question at math has been deleted. Obviously you were right about keeping the question here. I have received very helpful responses. $\endgroup$
    – a06e
    Commented Feb 19, 2012 at 16:09

4 Answers 4


A quick and general answers without mathematical abstractions. There are several options to impose boundary conditions, e.g.

Strictly speaking the Galerkin method requires that you choose a set of basis functions which satisfy the BC of the problem (e.g. via basis recombination and/or splitting of the approximation $u_h=u_0+u_N$ wit $u_0$ responsible for inhomogenous solutions and $u_N$ a partial sum which relies on basis functions which satisfies the homogenous conditions)

  • Penalty methods/Lagrange multiplies where one essentially add a penalty term which incorporated the boundary condition, e.g. $A + \tau \cdot B = b + \tau\cdot b_p$ where $B$ is a matrix responsible for the discrete boundary condition and $b_p$ is responsible for inhomogenous terms. In the limit $\tau\to\infty$ the conditions is strongly imposed and otherwise it is weakly imposed. Choice of $\tau$ affects conditioning of the system.

  • Tau method where a number of equations are exchanged (modification of rows in Galerkin system) with discrete versions of boundary conditions which is then enforced explicitly. Note: one option is also to make the systems overdetermined with additional boundary conditions.

  • Before discretization (Ritz Method) rewrite Galerkin formulation via Gauss divergence theorem to transform volume integrals to boundary integrals and then incorporate ( exact or approximately) boundary conditions directly in formulation before discretization.

  • Finally, by exploiting connection between nodal/modal expansions it is also possible to derive a nodal Galerkin method where the solution to the system is the coefficients of a Lagrange basis rather than a modal basis.

  • 1
    $\begingroup$ I guess $\tau$ is $\lambda$, isn't it? $\endgroup$
    – shuhalo
    Commented Feb 21, 2012 at 1:34
  • $\begingroup$ Yep. corrected. $\endgroup$ Commented Feb 21, 2012 at 20:50
  • 1
    $\begingroup$ Should it read "Galerkin method requires that you choose a set of basis functions which satisfy the BC of the problem"? $\endgroup$
    – knl
    Commented Apr 19, 2017 at 12:53
  • $\begingroup$ @knl: I think so as well, whereas the other sentence makes no sense. I'll do an edit. $\endgroup$
    – davidhigh
    Commented Mar 21, 2018 at 21:12

One possibility is to assemble the system matrix $A$ and right-hand side vector $b$, with the prescribed degrees of freedom as unknowns, like any other degree of freedom. Then, $A$ and $b$ are modified by zeroing rows and columns associated with the prescribed dofs, and putting a one into the corresponding diagonal entry, and appropriately modifying the rhs vector $b$.

When you zero rows, put one into the diagonal and change the rhs so that you enforce the prescribed value, the system is no longer symmetric. That´s why you zero columns and modify the rhs vector $b$ to account for the prescribed value.

Another possibility is to add a very large number $p$ (usually 1e10) to the diagonal of the prescribed dof and then set the rhs entry to p*$\bar{u}$, where $\bar{u}$ is the prescribed value of that dof.


The general problem of dealing with boundary conditions with the finite element method can get pretty complicated. But if:

  • $S(u)$ is such that the only imposition $S(u)=0$ makes on the form of $u$ is that it is equal to some $f(x,y)$ on $\delta \Omega$.

  • You can finagle your elements so that $\delta \Omega$ is entirely on the boundaries of various elements

it's actually very simple. Your equation:

$$\left( \begin{array}{ccc} \left\langle L\left[g_1\right],g_1\right\rangle & \ldots & \left\langle L\left[g_N\right],g_1\right\rangle \\ \ldots & \ldots & \ldots \\ \left\langle L\left[g_1\right],g_N\right\rangle & \ldots & \left\langle L\left[g_N\right],g_N\right\rangle \end{array} \right)\left( \begin{array}{c} a_1 \\ \ldots \\ a_N \end{array} \right)=0$$ needs to be replaced with $$\left( \begin{array}{ccc} \left\langle L\left[g_1\right],g_1\right\rangle & \ldots & \left\langle L\left[g_N\right],g_1\right\rangle \\ \ldots & \ldots & \ldots \\ \left\langle L\left[g_1\right],g_N\right\rangle & \ldots & \left\langle L\left[g_N\right],g_N\right\rangle \end{array} \right)\left( \begin{array}{c} a_1 \\ \ldots \\ a_N \end{array} \right)=\mathbf{b}$$

where the right hand side vector $\mathbf{b}$ represents the boundary conditions.

To determine $\mathbf{b}$, set the elements of your basis that determine the value of $u$ on $\delta \Omega$ to whatever values they need to be to satisfy the boundary conditions. In $\langle L[g_j],g_i \rangle$, you should exclude them from the $g_j$ but not the $g_i$ (the elements of $\mathbf{a}$ that correspond to these functions have already been determined, so they shouldn't be included in the matrix equation) . Then, set up $$\langle R,g_i \rangle = \sum_{j=1}^N a_j \langle L[g_j],g_i \rangle = 0$$ as a matrix equation, and the values of the elements of $\mathbf{b}$ should pop right out as the inner products of $L$ operating on your your interior basis with elements of your boundary basis.

  • $\begingroup$ Thanks for the answer Dan. I don't understand the paragraph that begins "To determine $\mathbf{b}$ ..." (which is the essential part, I think). Could you make it a little more explicit? $\endgroup$
    – a06e
    Commented Feb 18, 2012 at 4:28
  • $\begingroup$ On the other hand, the problem I'm trying to solve only satisfies the second condition you set: The boundary is a rectangle. As for the first condition, the boundary conditions do not specify the values of the function at the boundary. The boundary conditions specify the values of linear combinations of the second order derivatives of the function (something like $a\partial^2 u/\partial x\partial y +b\partial^2 u/\partial x ^2 = 0)$, where $a,b$ are constants. Also the boundary conditions are homogeneous. $\endgroup$
    – a06e
    Commented Feb 18, 2012 at 4:32
  • $\begingroup$ @becko: You might want to be more explicit about what $L$ and $S$ are in your question. Different types of boundary conditions can be handled in different ways. $\endgroup$
    – Dan
    Commented Feb 18, 2012 at 4:40
  • 1
    $\begingroup$ I edited the question a little to make it clearer, I think. I don't want to post the exact problem I'm trying to solve because I want to keep the question as general as I can. I think I can understand the method better that way. $\endgroup$
    – a06e
    Commented Feb 18, 2012 at 4:49
  • $\begingroup$ @becko: We might want to move this to chat, as it's getting kind of long. $\endgroup$
    – Dan
    Commented Feb 18, 2012 at 5:15

Here is a method known as basis recombination, which has not been mentioned in the present thread. I'm citing from the book of J.P. Boyd, "Chebyshev and Fourier Spectral Methods", 2nd Ed., Chapter 6.5.:

If the problem $$L u = f$$ has inhomogeneous boundary conditions, then it may always be replaced by an equivalent problem with homogeneous boundary conditions, so long as the boundary conditions are linear. The first step is to choose a simple function $B(x)$ which satisfies the inhomogeneous boundary conditions. One may then define a new variable $v(x)$ and new forcing function $g(x)$ via $$u(x) ≡ v(x) + B(x)\\ g(x) ≡ f(x) − L B(x)$$ so that the modified problem is $$L v = g$$ where $v(x)$ satisfies homogeneous boundary conditions. ...

The shift function $B(x)$ is arbitrary except for the constraint that it must satisfy all the inhomogeneous boundary conditions. However, the simplest choice is the best: polynomial interpolation of the lowest order that works.

Next comes my own explanation:

  • "Inhomogeneous boundary condition" means a condition which contains a constant, e.g. $$\partial_x u(x,y) |_{x=x_0} = 1.$$

    According to the above program, by choosing a convenient function $B(x)$, you get that down to $$\partial_x u(x,y) |_{x=x_0} = 0.$$

  • Once you made all boundary conditions homogeneous in this way, you can turn to your basis expansion (which I assume is done in terms of a product basis): $$ u(x,y) \ =\ \sum_{ij} a_{ij} \, \phi_{i}(x) \, \varphi_{j}(y) $$ By applying the corresponding BC operator, one obtains $$ \partial_x u(x,y) \ =\ \sum_{ij} a_{ij} \, \phi^\prime_{i}(x) \, \varphi_{j}(y) $$ and this should be zero for $x=x_0$ according to the above example.

  • Now comes the crucial step: by using a basis $\phi_i(x)$ which already satisfies the BC by itself, i.e. $\phi^\prime_i(x)|_{x=x_0} = 0$ for all $i$, the BC of the (transformed) two-dimensional problem are satisfied automatically! Basis sets of this and similar kind can be found, e.g., by a procedure called "basis recombination" (that is often used in combination with collocation methods).

  • Note that this is the point where homogeneous boundary conditions really matter, because otherwise one would need to impose further constraints. For example, suppose we would be working with the "$=1$" condition above, and, correspondingly, let's try to use a basis with $\phi^\prime_i(x)|_{x=x_0} = 1$. Then $$ \partial_x u(x,y)|_{x=x_0} \ =\ \sum_{ij} a_{ij} \, \varphi_{j}(y)\, $$ and in order to make this expression equal to $1$ for all $y$, one would have to constrain the expansion coeffcients $a_{ij}$ as well. Thus, for inhomogeneous BCs, there is a general way to apply the constraints to the one-dimensional parts but use it for to the full problem.

The nice thing about this whole approach is that it is working on a relatively abstract level. Necessary ingredients are only linearity of the BC operator and an ansatz in terms of product basis functions. As such, it is also applicable to approximate methods.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.