# Eigenfunctions of Laplacian and Helmholtz equation¶

## Wave equation with time-harmonic forcing¶

Let’s have wave equation with special right-hand side

(1)\begin{align}\begin{aligned}w_{tt} - \Delta w &= f\, e^{i\omega t} &&\quad\text{ in }\Omega\times(0,T),\\w &= 0 &&\quad\text{ on }\partial\Omega\times(0,T)\end{aligned}\end{align}

with $$f \in L^2(\Omega)$$. Assuming ansatz

$w(t, x) = u(x) e^{i\omega t}$

we observe that $$u$$ has to fulfill

(2)\begin{align}\begin{aligned}-\Delta u - \omega^2 u &= f &&\quad\text{ in }\Omega,\\u &= 0 &&\quad\text{ on }\partial\Omega.\end{aligned}\end{align}

Try solving (2) in FEniCS with data

(3)\begin{align}\begin{aligned}\Omega &= (0,1)\times(0,1),\\\omega &= \sqrt{5}\pi,\\f &= x + y\end{aligned}\end{align}

on series of refined meshes. Observe behavior of solution energy $$\|\nabla u\|_2$$ with refinement. Is there a convergence or not?

Define eigenspace of Laplacian (with zero BC) corresponding to $$\omega^2$$ as

$E_{\omega^2} := \biggl\{ u\in H_0^1(\Omega): -\Delta u = \omega^2 u \biggr\}.$

$$E_{\omega^2}\neq\{0\}$$ if and only if $$\omega^2$$ is an eigenvalue. Note that $$E_{\omega^2}$$ is finite-dimensional. Now define $$P_{\omega^2}$$ as $$L^2$$-orthogonal projection onto $$E_{\omega^2}$$. It is not difficult to check that the function

(4)$w(t, x) = \frac{t e^{i\omega t}}{2i\omega} (P_{\omega^2} f)(x) + e^{i\omega t} u(x)$

solves (1) provided $$u$$ fulfills

(5)\begin{align}\begin{aligned}-\Delta u - \omega^2 u &= (1-P_{\omega^2}) f &&\quad\text{ in }\Omega,\\u &= 0 &&\quad\text{ on }\partial\Omega.\end{aligned}\end{align}

Note that problem (5) has a solution which is uniquely determined up to arbitrary function from $$E_{\omega^2}$$.

Construct basis of $$E_{\omega^2}$$ by numerically solving the corresponding eigenproblem with data (3).

Hint

Having forms a, m and boundary condition bc representing eigenvalue problem

\begin{align}\begin{aligned}-\Delta u &= \lambda u &&\quad\text{ in }\Omega,\\u &= 0 &&\quad\text{ on }\partial\Omega.\end{aligned}\end{align}

assemble matrices A, B using function assemble_system

A = assemble_system(a, zero_form, bc)
B = assemble(m)


Then 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)
plt.show()


Implement projection $$P_{\omega^2}$$. Use it to solve problem (5) with data (3).

Construct the solution $$w(t, x)$$ of the wave equations (1) using formula (4). Plot temporal evolution of its real and imaginary part.

## Mesh generation by Gmsh¶

Modify a Gmsh demo to mesh a half ball

$\{(x,y,z), x^2 + y^2 + z^2 < 1, y>0\}$

using the following code:

wget https://gitlab.onelab.info/gmsh/gmsh/blob/ad0ab3d5c310e7048ffa6e032ccd4e8f0108aa12/demos/api/boolean.py
source /LOCAL/opt/gmsh-4.0.0/gmsh.conf
python3 boolean.py
meshio-convert -p -o xdmf-xml boolean.msh boolean.xdmf
paraview boolean.xdmf &

<edit> boolean.py

python3 boolean.py
meshio-convert -p -o xdmf-xml boolean.msh boolean.xdmf


If in a need peek into

>>> import gmsh


Find $$E_{\omega^2}$$ with $$\omega^2 \approx 70$$ on the half ball. Plot the eigenfunctions in Paraview.

Hint

Use Glyph filter, Sphere glyph type, decrease the scale factor to ca. 0.025.

Use Clip filter. Drag the clip surface by mouse, hit Alt+A to refresh.

## Reference solution¶

Show/Hide Code

Download Code

from dolfin import *
import matplotlib.pyplot as plt

def solve_helmholtz(V, lambd, f):
"""Solve Helmholtz problem

-\Delta u - lambd u = f  in \Omega
u = 0  on \partial\Omega

and return u.
"""

bc = DirichletBC(V, 0, lambda x, on_boundary: on_boundary)
u, v = TrialFunction(V), TestFunction(V)
L = f*v*dx
u = Function(V)
solve(a == L, u, bc)
return u

def build_laplacian_eigenspace(V, lambd, maxdim, tol):
"""For given space V finds eigenspace of Laplacian
(with zero Dirichlet BC) corresponding to eigenvalues
close to lambd by given tolerance tol. Return list
with basis functions of the space.
"""

# Assemble Laplacian A and mass matrix B
bc = DirichletBC(V, 0, lambda x, on_boundary: on_boundary)
u, v = TrialFunction(V), TestFunction(V)
L_dummy = Constant(0)*v*dx
m = u*v*dx
A, _ = assemble_system(a, L_dummy, bc)
B = assemble(m)

# Prepare eigensolver for
#
#    A x = lambda B x
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'] = float(lambd)
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['tolerance'] = 1e-6
#eigensolver.parameters['verbose'] = True  # for debugging

# Solve for given number of eigenpairs
eigensolver.solve(maxdim)

# Iterate over converged eigenpairs
space = []
for j in range(eigensolver.get_number_converged()):

# Get eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(j)

# Check that eigenvalue is real
assert near(c/r, 0, 1e-6)

# Consider found eigenvalues close to the target eigenvalue
if near(r, lambd, tol*lambd):
print('Found eigenfunction with eigenvalue {} close to target {} '
'within tolerance {}'.format(r, lambd, tol))

# Store the eigenfunction
eig = Function(V)
eig.vector()[:] = rx
space.append(eig)

# Check that we got whole eigenspace, i.e., last eigenvalue is different one
assert not near(r, lambd, tol), "Possibly don't have whole eigenspace!"

# Report
print('Eigenspace for {} has dimension {}'.format(lambd, len(space)))

return space

def orthogonalize(A):
"""L^2-orthogonalize a list of Functions living on the same
function space. Modify the functions in-place.
Use classical Gramm-Schmidt algorithm for brevity.
For numerical stability modified Gramm-Schmidt would be better.
"""

# Set of single function is orthogonal
if len(A) <= 1:
return

# Orthogonalize overything but the last function
orthogonalize(A[:-1])

# Orthogonalize the last function to the previous ones
f = A[-1]
for v in A[:-1]:
r = assemble(inner(f, v)*dx) / assemble(inner(v, v)*dx)
assert f.function_space() == v.function_space()
f.vector().axpy(-r, v.vector())

# Problem data
f = Expression('x + x', degree=1)
omega2 = 5*pi**2

# Iterate over refined meshes
ndofs, energies = [], []
for n in (2**i for i in range(2, 7)):

mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "Lagrange", 1)
u = solve_helmholtz(V, omega2, f)

# Store energy to check convergence
ndofs.append(u.function_space().dim())
energies.append(norm(u, norm_type='H10'))

# Plot energies against number dofs
plt.plot(ndofs, energies, 'o-')
plt.xlabel('dimension')
plt.ylabel('energy')
plt.show()

# Problem data
f = Expression('x + x', degree=1)
omega2 = 5*pi**2

# Iterate over refined meshes
ndofs, energies = [], []
for n in (2**i for i in range(2, 7)):

mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "Lagrange", 1)

# Build eingenspace of omega2
eigenspace = build_laplacian_eigenspace(V, omega2, 10, 0.1)

# Orthogonalize f to the eigenspace
f_perp = project(f, V)
orthogonalize(eigenspace+[f_perp])

# Find particular solution with orthogonalized rhs
u = solve_helmholtz(V, omega2, f_perp)

# Store energy to check convergence
ndofs.append(u.function_space().dim())
energies.append(norm(u, norm_type='H10'))

# Plot energies against number dofs
plt.plot(ndofs, energies, 'o-')
plt.xlabel('dimension')
plt.ylabel('energy')
plt.show()

# Create and save w(t, x) for plotting in Paraview
omega = omega2**0.5
Pf = project(f - f_perp, V)
T_per = 2*pi/omega
create_and_save_w(omega, Pf, u, 20*T_per, 0.1*T_per)

def create_and_save_w(omega, Pf, u, T, dt):
"""Create and save w(t, x) on (0, T) with time
resolution dt
"""

# Extract common function space
V = u.function_space()
assert V == Pf.function_space()

w = Function(V)

Pf_vec = Pf.vector()
u_vec = u.vector()
w_vec = w.vector()

f = XDMFFile('w.xdmf')
f.parameters['rewrite_function_mesh'] = False
f.parameters['functions_share_mesh'] = True

def c1(t):
return t*sin(omega*t)/(2*omega), -t*cos(omega*t)/(2*omega)

def c2(t):
return cos(omega*t), sin(omega*t)

t = 0
while t < T:
c1_real, c1_imag = c1(t)
c2_real, c2_imag = c2(t)

# Store real part
w.vector().zero()
w.vector().axpy(c1_real, Pf_vec)
w.vector().axpy(c2_real, u_vec)
w.rename('w_real', 'w_real')
f.write(w, t)

# Store imaginary part
w.vector().zero()
w.vector().axpy(c1_imag, Pf_vec)
w.vector().axpy(c2_imag, u_vec)
w.rename('w_imag', 'w_imag')
f.write(w, t)

t += dt

f.close()

if __name__ == '__main__':