5
$\begingroup$

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!

$\endgroup$
4
$\begingroup$

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:

zero-one

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:

one-two

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]]
$\endgroup$
0
2
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ 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? $\endgroup$ – Jolle Jan 26 '15 at 17:48
  • $\begingroup$ 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. $\endgroup$ – sebastian_g Jan 26 '15 at 17:54
  • 1
    $\begingroup$ You can't solve eigenvalue problems with inhomogenous boundary conditions. This is simply not well defined. $\endgroup$ – Wolfgang Bangerth Jan 27 '15 at 4:10
  • $\begingroup$ @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. $\endgroup$ – sebastian_g Jan 27 '15 at 8:46
  • 1
    $\begingroup$ @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. $\endgroup$ – Wolfgang Bangerth Jan 28 '15 at 3:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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