### 1-D eigenvalue problem with Neumann conditions

168

views

0

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$∂

$-u''\left(x\right)=\lambda\cdot u\left(x\right)$−

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.

Just to get going with Fenics, I am trying to compute solutions to the Neumann eigenvalue problem

for the operator $\partial^2_x$∂

^{2}_{x}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

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

What is confusing you is that you are seeing the eigenfunction corresponding to the

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

*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.

written
4 months ago by
Chris Judge

0

I think you are getting the largest eigenvalue. Try this small change

@ 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.

```
# Create eigensolver
eigensolver = SLEPcEigenSolver(A, B)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
```

written
4 months ago by
Chris Judge

Please login to add an answer/comment or follow this question.