Find the smallest eigenvalues and eigenfunctions of the follwing problem
with \(\Omega\) the unit circle for example.
from dolfin import * import mshr n=40 geometry=mshr.Circle(Point(0.0,0.0),1.0) mesh = mshr.generate_mesh(geometry,n) V = FunctionSpace(mesh, 'Lagrange', 1) bc = DirichletBC(V, 0.0, DomainBoundary()) u, v = TrialFunction(V), TestFunction(V) a = inner(grad(u), grad(v))*dx L = Constant(0.0)*v*dx m = u*v*dx A, _ = assemble_system(a, L, bc) B = assemble(m) print "Discrete space size: %d"%V.dim() eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B)) #prm = eigensolver.parameters #info(prm, True) eigensolver.parameters['spectrum'] = 'smallest magnitude' eigensolver.solve(20) eig = Function(V) eig_vec = eig.vector() print 'converged: %d'%eigensolver.get_number_converged() for j in range(eigensolver.get_number_converged()): r, c, rx, cx = eigensolver.get_eigenpair(j) eig_vec[:] = rx print r, c plot(eig, interactive=True)
Let’s have wave equation with special right-hand side
with \(f \in L^2(\Omega)\). Such a problem has a solution (in some proper sense; being unique when enriched by initial conditions), see [Evans], chapter 7.2.
Task 1. Try seeking for a particular solution of this equation while taking advantage of special structure of right-hand side. Assuming ansatz
\[w := u\, e^{i\omega t}, \quad u\in H_0^1(\Omega)\]derive non-homogeneous Helmholtz equation for \(u\) using the Fourier method and try solving it using FEniCS with
- \(\Omega = [0,1]\times[0,1]\),
- \(\omega = \sqrt{5}\pi\),
- \(f = x + y\).
Task 2. Plot solution energies against number of degrees of freedom.
Hint. Having list of number of degrees of freedom ndofs and list of energies energies do
import matplotlib.pyplot as plt plt.plot(ndofs, energies, 'o-') plt.show()What does it mean? Is the problem well-posed?
Define eigenspace of Laplacian (with zero BC) corresponding to \(\omega^2\)
If \(E_{\omega^2}\neq{0}\) then \(\omega^2\) is eigenvalue. Then by testing the non-homogeneous Helmholtz equation (derived in previous section) by non-trivial \(v\in E_{\omega^2}\) one can see that \(f\perp E_{\omega^2}\) is required (check it!), otherwise the problem is ill-posed. Hence the assumed ansatz is generally wrong. In fact, the condition \(f\perp E_{\omega^2}\) is sufficient condition for well-posedness of the problem, see [Evans], chapter 6.2.3.
The resolution is to seek for a particular solution for \(f^\parallel\) and \(f^\perp\) (\(L^2\)-projections of \(f\) to \(E_{\omega^2}\) and \(^\perp E_{\omega^2}\) respectively) separately. As \(E_{\omega^2}\) has finite dimension (due to the Fredholm theory), the former can be obtained by solving forced harmonic oscillator equation which is easily done in hand (once \(E_{\omega^2}\) is known). This is equivalent to picking blowing-up ansatz \(w = u\, t\, e^{i t\omega},\, u\in H_0^1(\Omega)\). So let’s focus to the latter part.
Task 3. Use SLEPc eigensolver to find \(E_{\omega^2}\).
Hint. Having assembled matrices A, B, the eigenvectors solving
\[A x = \lambda B x\]with \(\lambda\) close to target lambd can be found by
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B)) eigensolver.parameters['problem_type'] = 'gen_hermitian' eigensolver.parameters['spectrum'] = 'target real' eigensolver.parameters['spectral_shift'] = lambd eigensolver.parameters['spectral_transform'] = 'shift-and-invert' eigensolver.parameters['tolerance'] = 1e-6 #eigensolver.parameters['verbose'] = True # for debugging eigensolver.solve(number_of_requested_eigenpairs) eig = Function(V) eig_vec = eig.vector() space = [] for j in range(eigensolver.get_number_converged()): r, c, rx, cx = eigensolver.get_eigenpair(j) eig_vec[:] = rx plot(eig, title='Eigenvector to eigenvalue %g'%r) interactive()Task 4. Write function which takes a tuple of functions and \(L^2\)-orthogonalizes them using Gramm-Schmidt algorithm.
Task 5. Compute \(f^\perp\) for \(f\) from Task 1 and solve the Helmholtz equation with \(f^\perp\) on right-hand side. Again, plot energies of solutions against number of degrees of freedom.
Note
Lecturer note. Student must not include eigenvectors corresponding to other eigenvalues. SLEPc returns these after last targeted one. For this case the dimension of \(E_{\omega^2}\) is 2. Let’s denote this bunch of vectors by E.
GS orthogonalization is called to tuple E+[f]. This first orthogonalizes eigenvectors themself (for sure – SLEPc doc is not conclusive about this) and then orthogonalizes f to \(E_{\omega^2}\).
| [Evans] | (1, 2) Lawrence C. Evans. Partial Differential Equations. Second edition. 1998, 2010 AMS, Rhode Island. |
from dolfin import *
def helmholtz_illposed(n):
"""For given mesh division 'n' solves ill-posed problem
(-Laplace - 5*pi^2) u = x + y on [0, 1]*[0, 1],
u = 0 on boundary,
and returns space dimension, energy_error (on discrete subspace) and energy."""
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'Lagrange', 1)
bc = DirichletBC(V, 0.0, lambda x, b: b)
u, v = TrialFunction(V), TestFunction(V)
a = (inner(grad(u), grad(v)) - Constant(5.0*pi*pi)*u*v)*dx
L = Expression('x[0] + x[1]')*v*dx
u = Function(V)
solve(a == L, u, bc)
energy_error = assemble(action(Constant(1.0)*action(a, u) - L, u))
energy = assemble(action(Constant(0.5)*action(a, u) - L, u))
return V.dim(), energy_error, energy
def helmholtz_wellposed(n):
"""For given mesh division 'n' solves well-posed problem
(-Laplace - 5*pi^2) u = f on [0, 1]*[0, 1],
u = 0 on boundary,
with f orthogonal to eigenspace of 5*pi^2.
and returns space dimension, energy_error (on discrete subspace) and energy."""
# Assemble Laplacian
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, 'Lagrange', 1)
bc = DirichletBC(V, 0.0, lambda x, b: b)
u, v = TrialFunction(V), TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = Constant(0.0)*v*dx
m = u*v*dx
A, _ = assemble_system(a, L, bc)
B = assemble(m)
# Search for eigenspace for eigenvalue close to 5*pi*pi
# NOTE: A x = lambda B x is proper FE discretization of the eigenproblem
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['spectral_shift'] = 5.0*pi*pi
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['tolerance'] = 1e-6
#eigensolver.parameters['verbose'] = True
eigensolver.solve(5) # Find 5 eigenpairs close to target
eig = Function(V)
eig_vec = eig.vector()
space = []
for j in range(eigensolver.get_number_converged()):
r, c, rx, cx = eigensolver.get_eigenpair(j)
assert near(c/r, 0.0, 1e-6)
assert near(cx.norm('linf')/rx.norm('linf'), 0.0, 1e-6)
if near(r, 5*pi*pi, 0.5*pi*pi):
eig_vec[:] = rx
space.append(eig.copy(deepcopy=True))
# Check that we got whole eigenspace - last eigenvalue is different one
assert not near(r, 5*pi*pi, 0.5*pi*pi), \
"Possibly don't have whole eigenspace!"
print 'Eigenspace for 5*pi^2 has dimension', len(space)
# Orthogonalize right-hand side to 5*pi^2 eigenspace
f = Expression('x[0] + x[1]')
f = project(f, V)
orthogonalize(space+[f])
# Solve well-posed resonant Helmoltz system
a = (inner(grad(u), grad(v)) - Constant(5.0*pi*pi)*u*v)*dx
L = f*v*dx
u = Function(V)
solve(a == L, u, bc)
energy_error = assemble(action(Constant(1.0)*action(a, u) - L, u))
energy = assemble(action(Constant(0.5)*action(a, u) - L, u))
return V.dim(), energy_error, energy
def orthogonalize(A):
"""L^2-orthogonalizes set of Functions A. Stores the result in-place to A.
Uses classical Gramm-Schmidt algorithm for brevity. For numerical stability
modified Gramm-Schmidt would be better."""
assert all(isinstance(a, Function) for a in A)
if len(A) <= 1:
return
orthogonalize(A[:-1])
f = A[-1]
for v in A[:-1]:
# NOTE: L^2 inner product could be preassembled to reduce computation
# of r to matvecs.
r = assemble(inner(f, v)*dx)/assemble(inner(v, v)*dx)
if f.function_space() == v.function_space():
f.vector().axpy(-r, v.vector())
else:
raise NotImplementedError
if __name__ == '__main__':
import numpy as np
import matplotlib.pyplot as plt
# Demonstrate that energy of ill-posed Helmholtz goes to minus infinity
results = np.array(map(helmholtz_illposed, [2**i for i in range(2, 9)]))
plt.subplot(2, 1, 1)
plt.plot(results[:, 0], results[:, 2], 'o-')
plt.title('Ill-posed Helmholtz')
plt.xlabel('dimension')
plt.ylabel('energy')
plt.show(block=False)
# Demonstrate that energy of well-posed Helmholtz converges
results = np.array(map(helmholtz_wellposed, [2**i for i in range(2, 9)]))
plt.subplot(2, 1, 2)
plt.plot(results[:, 0], results[:, 2], 'o-')
plt.title('Well-posed Helmholtz')
plt.xlabel('dimension')
plt.ylabel('energy')
plt.tight_layout()
plt.show(block=True)