# How to impose boundary conditions on eigenfunction problems?

I am trying to solve for the eigenfunctions of a (1D) differential operator using finite differences:

$$A \, f(x) = \lambda f(x)$$

Here is an example in Python where $A = \partial_x^4$:

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import spdiags
from scipy import linalg

def derivative_matrix(n):
stencil = np.array((1,-4,6,-4,1))
diags = range(-2,3)
bands = np.tile(stencil,(n,1)).T
A = spdiags(bands,diags,n,n).todense()
A = np.array(A,dtype=float)
return A

if __name__=="__main__":
n = 100
A = derivative_matrix(n)
vals,vecs = linalg.eig(A)
plt.plot(np.real(vecs[:,-3:None]))
plt.show()


This gives me really nice eigenfunctions:

My questions is: what do I have to change in the implementation to change the boundary condition of $f$? Currently, the method seems to implicitely assume I want $f(x)=f'(x)=0$. But what if I instead wanted $f'(x)=f''(x)=0$? How would I implement that?

Any help will be much appreciated. Thanks in advance!

Consider what happens when you approximate a function's derivative using finite differences near a boundary. If the boundary is at the point $x_0$, the point $x_1$ is just outside the boundary, then the approximation at the point just before the boundary, $x_{-1}$ is: $$h^{-4}(f_{-3} -4f_{-2}+6f_{-1}-4f_0+f_1),$$ where $f_0$ and $f_1$ must be given by boundary conditions.

The boundary condition $f=f'=0$ can be approximated to order $O(h^2)$ with the conditions $$f_0 = 0, \qquad \frac{f_1 - f_{-1}}{2h} = 0,$$ using the usual finite difference formulas. Therefore, to get the boundary condition $f=f'=0$, you adjust the matrix coefficients at the point $x_{-1}$ by substituting the $f_0$ and $f_1$ that you get from BCs: $$f_{-3}-4f_{-2}+7f_{-1}.$$ Note that this is different from your code, where you use $6$ instead of $7$. If you work it out backwards, that's a bit like using $2f'+hf''=0$ instead of $f'=0$. Here is what I get:

For the boundary condition $f'=f''=0$, the procedure is the same. The approximate boundary conditions are: $$\frac{f_1-f_{-1}}{2h} = 0, \qquad \frac{f_1-2f_0+f_{-1}}{h^2} = 0,$$ so that $f_1 = f_0 = f_{-1}$, and therefore the stencil coefficients are $$x_{-2}:\ [1, -4, 6, -3], \qquad x_{-1}:\ [1, -4,3].$$ Here is the outcome:

Matrix corners for case 1:

[[ 7 -4  1  0  0  0  0]    [[ 6 -4  1  0  0  0  0]
[-4  6 -4  1  0  0  0]     [-4  6 -4  1  0  0  0]
[ 1 -4  6 -4  1  0  0]     [ 1 -4  6 -4  1  0  0]
[ 0  1 -4  6 -4  1  0]     [ 0  1 -4  6 -4  1  0]
[ 0  0  1 -4  6 -4  1]     [ 0  0  1 -4  6 -4  1]
[ 0  0  0  1 -4  6 -4]     [ 0  0  0  1 -4  6 -4]
[ 0  0  0  0  1 -4  6]]    [ 0  0  0  0  1 -4  7]]


Matrix corners for case 2:

[[ 3 -4  1  0  0  0  0]    [[ 6 -4  1  0  0  0  0]
[-3  6 -4  1  0  0  0]     [-4  6 -4  1  0  0  0]
[ 1 -4  6 -4  1  0  0]     [ 1 -4  6 -4  1  0  0]
[ 0  1 -4  6 -4  1  0]     [ 0  1 -4  6 -4  1  0]
[ 0  0  1 -4  6 -4  1]     [ 0  0  1 -4  6 -4  1]
[ 0  0  0  1 -4  6 -4]     [ 0  0  0  1 -4  6 -3]
[ 0  0  0  0  1 -4  6]]    [ 0  0  0  0  1 -4  3]]


I solved a similar problem a few days ago. I am going to present the implementation for a simpler case. Consider the eigenvalue problem equation $$-y''=\lambda y\,,\quad x\in\left(a,b\right)\\ y\left(a\right)=0\,\quad y'\left(b\right)=y\left(b\right)$$ The boundary condition at $x=a$ is a homogenous Dirichlet boundary condition. In finite differences we don't have to worry too much about this one.

The Boundary condition at $x=b$ is called Robin boundary condition. The value $y\left(b\right)$ is unknown. Let $y_i$ be the discrete solution at the $i$-th grid point i.e. $y_i=y\left(x_i\right)\,,\, x_i=ih\,,\,i>0\,.$ Using symmetric 2nd order finite differences, we may write the derivative at $x=b$, i.e. $x=x_N$, as $$y'_N=\frac{y_{N+1}-y_{N-1}}{2h}\,.$$ The value $y_{N+1}$ is outside of the domain. Using the boundary Robin condition we obtain $$y'_N=y_N\Leftrightarrow y_{N+1}=y_{N-1}+2h y_{N}\,.$$ Substitute in formulas for 2nd derivative yields $$y''_{N}=\frac{y_{N+1}-2y_N+y_{N-1}}{h^2} =\frac{-2\left(1-h\right)y_{N}+2y_{N-1}}{h^2}\,.$$ In conclusion just modify the last line of your finite difference matrix for this case.

I hoped I could show how to do it. I think the scheme is straight forward for symmetric 4th order finite differences.

• This definitely makes sense, but I don't see how to generalize this to, say, y'(b)=k, where k is some constant. Your method only works because it is a Robin BC, don't you agree? – Jolle Jan 26 '15 at 17:48
• That is correct. To admit, I asked myself the same question. But didn't find a satisfying answer. Actually in physics/ engineering I never had to solve an eigenvalue problem with inhomogeneous boundary conditions. Luckily, boundary conditions you gave above are homogeneous. – sebastian_g Jan 26 '15 at 17:54
• You can't solve eigenvalue problems with inhomogenous boundary conditions. This is simply not well defined. – Wolfgang Bangerth Jan 27 '15 at 4:10
• @WolfgangBangerth That is good remark. Can you recommend literature about eigenvalue problems? This topic is neihter covered intensively in standard engineering math books nor in books on numerics of odes/pdes for beginners. I just remember a book by Lothar Collatz Eigenvalueproblems with Technical Applications, german title Eigenwertaufgaben mit technischen Anwendungen . But this is an ancient book from 1949. – sebastian_g Jan 27 '15 at 8:46
• @sebastian_g: I don't really have a good suggestion. If you think about where you want to solve eigenvalue problems, it's almost always because you want to expand something in a series (i.e., write one quantity as a sum of multiples of others) and in these cases it doesn't make sense if you have a nonzero right hand side defining the "other quantities" because you're going to take (as yet unspecified) multiples of them. – Wolfgang Bangerth Jan 28 '15 at 3:32