# Boundary conditions or other problems

Hi

I’m trying to solve a 2d Helmholtz scattering equation for its resonance values for n=2 and other values:
$$\Delta u_1+k^2(1+n)u_1=0 \in D$$
$$\Delta u_2+k^2u_2=0 \in \mathcal{R}^2 \D$$
with Sommerfeld radiation condition and equality of the u’s on the boundary of D and their partials with respect to the normal.

This turns into a quadratic eigenvalue problem and I use the following code to try and solve for k.

My code is:

#Import needed packages and programs
import netgen.gui
from netgen.geom2d import SplineGeometry
from ngsolve import *
import matplotlib.pyplot as plt
import netgen.geom2d as geom2d
from netgen.geom2d import unit_square
import numpy as np
import scipy.sparse as scs
from scipy.sparse.linalg import eigs, eigsh

#Create geometries, first circle is the scatterer, second begins the PML third truncates the infinite domain
#third circle (largest) gets a boundary condition due to the sommerfeld radiation condition.
geo = SplineGeometry()
geo.AddCircle( (0.0, 0.0), r=1.0, leftdomain=1, rightdomain=2, bc = "scatterer")
geo.AddCircle ( (0.0, 0.0), r=1.5, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.0, 0.0), r=2.0, leftdomain=3, rightdomain=0,bc = "outerbnd")

#Set material 1 to be inner do that only those nodes inside the scatterer recieve the $n*u*v*dx("inner")$ term
geo.SetMaterial(1, "inner")
geo.SetMaterial(3, "PML")

#Generate the mesh
ngmesh = geo.GenerateMesh(maxh=0.04)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle

#Set the extra value of the scatterer
nn=2;

#I think this tells the program what level of accuracy and where our test/trial functions live
fes = H1(mesh, order=4, complex=True, dirichlet="dir")

#create our trial function u and our test function v
u = fes.TrialFunction()
v = fes.TestFunction()

#create matrix K that is the non lambda term (sign for quad-eig)
k = BilinearForm (fes, symmetric=True)

#create matrix C that is our lambda term (sign for quad-eig)
c = BilinearForm (fes, symmetric=True)
c += 1.j*u*v*ds('outerbnd')

#create matrix M which is our squared term and has the nn element (sign for quad-eig)
m = BilinearForm (fes, symmetric=True)
m += u*v*dx
m += nn*u*v*dx('inner')

#convert bilinearform into a matrix form
k.Assemble()
c.Assemble()
m.Assemble()

#tells the program to get ready
u = GridFunction(fes, multidim=20, name='resonances')

#moves the matrices from ngsolve into scipy.sparse because ngsolve hates stacking things
mcsr = scs.csr_matrix(m.mat.CSR(), shape=(m.mat.height, m.mat.width),dtype='complex')
ccsr = scs.csr_matrix(c.mat.CSR(), shape=(c.mat.height, c.mat.width),dtype='complex')
kcsr = scs.csr_matrix(k.mat.CSR(), shape=(k.mat.height, k.mat.width),dtype='complex')

#finds the size of M and creates the needed identity matrix of that size
leng = m.mat.height
I = scs.csr_matrix(scs.eye(leng), dtype='complex')

#Builds our first block matrix in csc form
B1 = scs.csr_matrix(scs.bmat([[mcsr, None], [None, I]], dtype='complex'))

#Builds our second block matrix in csc form
B2 = scs.csr_matrix(scs.bmat([[ccsr, kcsr],[I, None]], dtype='complex'))

omegaval, omegavec = eigs(B2, k=11, M=B1, sigma=2.2-0.2*1.j)

for ome in omegaval:
print(ome)


This results in the following output results:
(2.1782707705367237-0.5074461983082241j)
(2.2227185408797476-0.5094990190123052j)
(2.2675293712400597-0.4965898981642067j)
(2.306187573060749-0.49327372936216596j)
(2.3375387477427036-0.4800143486380862j)
(2.3939839608351714-0.44250305302787263j)
(2.250056254487997+0.10549392553008125j)
(2.2934973319036587+0.09638866757241127j)
(2.3307875643942353+0.08609849468627778j)
(2.386139969059608+0.060068995693344174j)
(2.372968202422313+0.07669589505376706j)

The problem is I can solve for the resonances using Hankel and Bessel functions and a Newton solver to get that a resonance values should be at 2.284220485273998-0.3840888660574083j. I have verified this will a friend and a Muller’s method solver.

First question is does the bilinear form account for the equality along the boundary (in both u and its derivative wrt the normal) or is there something special I need to do in the code

Second, does anyone have an idea as to how to fix the discrepancy?