I've been reading multiple papers and related posts for a while now, but I can't seem to find a specific answer to the issues I'm having so I hope someone can clarify things here. I'll provide some background first before diving into the issues at hand.


Basically, I'm developing a 1D finite element code which calculates eigenvalues and eigenvectors after performing a Fourier analysis on a set of equations. There are quite a few equations/variables, but for simplicity say there are 3 of each and name the variables $u_1, u_2$ and $u_3$. I'm using either quadratic or cubic basis functions depending on the variable (see the figure below, left is quadratic, right is cubic)

Hence it basically comes down to this for a single variable: due to the nature of the basis functions we have a 4x4 matrix, subdivided into four 2x2 blocks that are combining the basis functions, where for example $h3$ stands for $C3$ or $Q3$ depending on the variable.

Here $\alpha$ is a combination of the finite element integrals evaluated using Gaussian quadrature etc. Hence, for 3 variables, this yields a 12x12 block for every point in the grid. Matrix assembly is quite straightforward using the above scheme, and eventually this results in a matrix eigenvalue problem of the form $ A\mathbf{x} = \omega B \mathbf{x}$, where $A$ and $B$ are block-tridiagonal matrices. $B$ is real and symmetric, $A$ is complex and is, in general, not symmetric. Say for simplicity that the domain is $x \in [0, 1]$.

Boundary conditions

Partial integration of the finite element integrals introduces natural boundary conditions, which are already added to their respective elements in the $A$ matrix during assembly. Applying for example Dirichlet boundary conditions on both sides of the domain is easy for 'standard' finite element problems. However, here we have a general eigenvalue problem. This post contains a good overview of the possible approaches. I should note that I use method 2 described there since my problem can vary, which implies that reducing the system would mean handling the different cases separately and I would like to keep it general. Also, the matrix $B$ is not allowed to be singular, since it boils down to solving the problem $B^{-1}A\mathbf{x} = \omega \mathbf{x}$, as for example described here.

Now, assume that $u_1$ and $u_2$ are quadratic, and $u_3$ is cubic. The Dirichlet conditions could state that $u_2(0) = u_2(1) = u_3(0) = u_3(1) = 0$, and I'm handling that by omitting the basis functions that are not zero on the boundaries. For the cubic ones those are $C1$ and $C2$, for the quadratic ones $Q3$ and $Q4$. Looking at the grid above this means setting the odd rows/columns to 0 for the cubic variable $u_3$ and the even rows/columns to 0 for the quadratic variable $u_2$. Nothing is done for $u_1$. Afterwards a 1 is introduced on $B$'s diagonal and a large number on $A$'s diagonal to yield $u_i = (\omega / s)u_i$ where $s$ is the large number, resulting in $u_1 = 0$ for $\omega \neq s$. This works out quite nicely, since the eigenfunctions of $u_2$ and $u_3$ are zero on the boundaries as expected and look fine.

The question

Say I want to replace one Dirichlet condition on the right side by a Neumann condition, for example $\partial_x u_2(1) = 0$. This is again not hard for a 'standard' problem, but I don't see how to do this for an eigenvalue problem. How would one approach implementing this? I'm having two thoughts about this:

  1. Zero out rows/columns containing the derivative of the basis functions for that specific variable.
  2. Perhaps these enter through the natural boundary conditions, which would have to be handled explicitly in that case instead of just adding them to the A-matrix. See also this post.

Is one of these the correct way to handle Neumann conditions?

I'm using LAPACK routines and Fortran to assemble and solve the system, by the way. I'm also linking these posts, which are related (and quite informative) but didn't help:

  1. Generalised eigenvalue problem
  2. Quadratic eigenvalue problem
  3. Boundary conditions and Galerkin method

Edit: As per @nicoguaro's suggestion I split the OP and kept the most relevant question here, the others will be linked later on.

  • $\begingroup$ I would suggest that you split your question into three different posts. That way it would be easier to answer them and more useful as well. $\endgroup$
    – nicoguaro
    Commented May 6, 2020 at 16:33
  • $\begingroup$ @nicoguaro that's a valid point, the idea was to keep them together since they are closely related, but I'm starting to think it's too much information for a single post. I'll edit it and focus on the third part here, that one is bugging me the most. $\endgroup$
    – n-claes
    Commented May 6, 2020 at 17:00
  • $\begingroup$ I think that you can post three questions and mention that they are related and add the hyperlinks to the other ones. $\endgroup$
    – nicoguaro
    Commented May 6, 2020 at 17:02
  • $\begingroup$ The thing is that I might have the answer to one of the questions, but I feel discouraged to not answer because it does not (completely) answer it. $\endgroup$
    – nicoguaro
    Commented May 6, 2020 at 17:02
  • $\begingroup$ I'm no sure I understand the question. Is it: how to apply Neumann boundary conditions for FEM in eigenvalue problems? $\endgroup$
    – nicoguaro
    Commented May 7, 2020 at 0:11

1 Answer 1


The boundary conditions don't depend on the choice of your basis but on the formulation you have for your problem. If you have a "standard" finite element formulation, you don't need to do anything to apply (homogeneous) Neumann boundary conditions, they are already satisfied by your system.

In the most common formulation, Neumann boundary conditions are natural boundary conditions. That means that they appear in your variational formulation and are satisfied implicitly by your problem. For the Helmholtz equation, the weak form is something like the following

$$\int\nabla u\nabla w\, \mathrm{d}V + \int w\frac{\partial u}{\partial \hat{\mathbf{n}}}\,\mathrm{d}S = \omega^2 \int u w\, \mathrm{d}V\, ,$$

but in the case of homogeneous boundary conditions, that you should have in an eigenvalue problem, the second term in the left-hand-side is zero because $\partial u/\partial \hat{\mathrm{n}}=0$.

  • 1
    $\begingroup$ Right, I know that the choice of basis functions does not influence the boundary conditions. I am a bit puzzled though by your answer "you don't need to do anything, they are already satisfied by your system". Dirichlet boundary conditions have to be explicitly implemented, so why shouldn't this be the case for Neumann BC's as well? $\endgroup$
    – n-claes
    Commented May 13, 2020 at 15:30
  • $\begingroup$ Aah I see, this makes sense. In the case of homogeneous Neumann conditions this would indeed automatically be satisfied due to the weak form. And this also solves a follow-up question I was having about applying nonhomogeneous boundary conditions, in which case I just have to take care of the additional terms. Thanks! $\endgroup$
    – n-claes
    Commented May 13, 2020 at 20:25

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.