How to impose boundary conditions on eigenfunction problems?

Question

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!

Answer 1

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]]

Answer 2

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.