Hi

I’m trying to solve a 2d Helmholtz scattering equation for its resonance values for n=2 and other values:

[tex]\Delta u_1+k^2(1+n)u_1=0 \in D[/tex]

[tex]\Delta u_2+k^2u_2=0 \in \mathcal{R}^2 \D[/tex]

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
mesh.SetPML(pml.Radial(rad=1.5,alpha=1j,origin=(0,0)), "PML")
#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)
k += grad(u)*grad(v)*dx
#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?