# Boundary conditions in a finite element eigenvalue problem

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.

## Background

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:

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

• 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. May 6, 2020 at 16:33
• @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. May 6, 2020 at 17:00
• I think that you can post three questions and mention that they are related and add the hyperlinks to the other ones. May 6, 2020 at 17:02
• 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. May 6, 2020 at 17:02
• I'm no sure I understand the question. Is it: how to apply Neumann boundary conditions for FEM in eigenvalue problems? May 7, 2020 at 0:11

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