1-D eigenvalue problem with Neumann conditions


168
views
0
4 months ago by
I am a complete novice (a pure mathematician).
Just to get going with Fenics, I am trying to compute solutions to the Neumann eigenvalue problem
for the operator   $\partial^2_x$2x  on the interval [0,1]. Namely, I want to find  $\lambda$λ  and u so that

$-u''\left(x\right)=\lambda\cdot u\left(x\right)$u''(x)=λ·u(x) with  u'(0)=1 and u'(1)=0.

The solutions ( $\lambda$λ,u) to this problem consist of  $\lambda$λ = (k \pi)^2 and u which
are scalar multiples of cos(k \pi  x) where k is an integer.

The code below produces a "solution" u that "interpolates" between the solutions -cos(\pi  x) and cos(\pi x).
I would have expected either cos(\pi x) or -cos(\pi  x), but not both superimposed in this odd way! 
I'd be grateful for an explanation.  
from __future__ import print_function
from dolfin import *

# Define mesh, function space
n = 250
mesh = UnitIntervalMesh(n)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define basis and bilinear forms
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
b = u*v*dx

# Assemble stiffness forms
A = PETScMatrix()
B = PETScMatrix()
assemble(a,tensor=A)
assemble(b,tensor=B)

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)

# Compute all eigenvalues of A x = \lambda x
print("Computing eigenvalues. This can take a minute.")
eigensolver.solve()

# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(1)
print("eigenvalue = ", r)

# Initialize function and assign eigenvector
u = Function(V)
u.vector()[:] = rx

# Plot eigenfunction
plot(u)
interactive()​

Community: FEniCS Project

3 Answers


3
4 months ago by
By default, the SLEPc eigensolver solves for the eigenvalues in descending order.  What  you're plotting is the second-highest frequency mode in the discrete system (the highest being that obtained from get_eigenpair(0)).  The high-frequency modes are poorly approximated, since the mesh cannot capture the rapid spatial variation of the exact solution.  (This is why some time integrators for semidiscrete PDEs, like the generalized-\(\alpha\) method, deliberately add artificial damping at high frequencies.)  To see the modes that you're probably expecting, you can set the parameter
# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)
eigensolver.parameters["spectrum"]="smallest magnitude"

to get the lowest frequencies first.  The \(0^{\text{th}}\) mode will be flat to within machine precision (so the plot may look like noise after it automatically re-scales the \(y\)-axis), then the \(1^{\text{st}}\), \(2^{\text{nd}}\), \(3^{\text{rd}}\), etc. modes will be (nearly) your textbook cosine waves, until the mesh can't approximate them anymore.  

@ David Kemeny: Thanks very much!  Perfect cosine now appears.   
written 4 months ago by Chris Judge  
3
4 months ago by
What is confusing you is that you are seeing the eigenfunction corresponding to the largest eigenvalue, not the smallest eigenvalue.  If you add the line

eigensolver.parameters["spectrum"] = "smallest magnitude"
​

just before the eigensolver.solve() command, you will see the function you expect.

To solve eigenvalue problems effectively, you need to familiarize yourself with several parameters to eigensolver.  Not only "spectrum", but also "solver", "problem_type", "spectral_transform", and "spectral_shift".  These are well documented in the documentation for SLEPc.
@ Doug Arnold: Thanks! I guess that the echo of the shape of the k^th eigenfunction in the (n-k)^th eigenfunction is due to some symmetry...  
written 4 months ago by Chris Judge  
0
4 months ago by
I think you are getting the largest eigenvalue. Try this small change

# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)
eigensolver.parameters['spectrum'] = 'smallest magnitude'​
@ Praveen C: Thanks!  If David Kamensky is correct, then it would actually be the next to highest. I am looking for next to lowest not the lowest.
written 4 months ago by Chris Judge  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »