# Stability issues a epsilon goes to zero for resonances

Hello,

I’m running some resonance experiments with a base problem and a perturbed problem. I calculate the resonance in the base case (first block of code at the bottom) and then run the perturbed case in block 2 and 4 changing the choices of epsilon. My issue is I know from a 1d example that I need epsilon in the range of 1/20-1/40 to obtain the numerical results that coincide with the analytical results that I already have. The idea is that the epsilon cases are perturbed resonance values that converge to the base resonance as epsilon goes to 0 along the path delta.

The issue is I somehow need to strengthen the code to handle the smaller epsilons but I am not sure how. I have tried expanding the PML final layer further out but that does not work. I have also tried to shrink the layers and shrink the max mesh size but have issue obtaining the base resonance (which is not perfect but close to the true values, verified through alternate means). Any help would be appreciated. I tried attaching a jupyter notebook but it would not attach. As such, I have the code posted below.

Block 1

[code]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

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)
geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)

#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.1)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle
#Set the extra value of the scatterer
nn=10;

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

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

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

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

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

u = GridFunction(fes, multidim=1, name=‘resonances’)

lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
u.vecs, shift=0.1502419606698783-0.19383557950824307*1.j)
Draw(u)

print ("lam: ", lam)
tplam=lam[0]
print(tplam)[/code]

Block 2:

[code]siz=5;
tpepsvec = np.zeros(siz);
tplamvec = np.zeros(siz, dtype = complex);
for i in range(siz):
eps = 1/(15+i+0.05);
tpepsvec[i]=eps;

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)
geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)
#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.1)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle
#Set the extra value of the scatterer
#Set the extra value of the scatterer
f = 10+2np.pisin(2np.pix/eps)+2np.picos(2np.piy/eps)

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

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

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

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

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

u = GridFunction(fes, multidim=5, name='resonances')

lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
u.vecs, shift=tplam)
for l in lam:
if abs(l-(tplam))<0.05:
tplamvec[i]=lam[0]


print(tplamvec) [/i][/i][/code]

[i][i]Block 3:

[code]N=5-1;

sq = np.zeros(N+1);
for j in range(N+1):
sq[j]=np.abs(tplamvec[j]-tplam);
print(sq)

for i in range(N+1):
ax = plt.gca()
ax.scatter(tpepsvec[i] ,sq[i] , c=‘blue’);
ax.set_yscale(‘log’);
ax.set_xscale(‘log’);

m, b = np.polyfit(np.log(tpepsvec), np.log(sq), 1);

plt.plot(tpepsvec, np.exp(b)*tpepsvec**m)
plt.title(‘\delta = 0.05, on loglog scale’)
plt.xlabel(‘Epsilon’)
plt.ylabel('Resonance values ')
plt.show()
print(‘The slope is:’, m)[/i][/i][/code]

The slope output is [/i][/i]0.29766692963062147[i][i][i][i]Block 4:

[code]siz=5;
tepsvec = np.zeros(siz);
tlamvec = np.zeros(siz, dtype = complex);
for i in range(siz):
eps = 1/(2+i+0.05);
tepsvec[i]=eps;

#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.AddRectangle( (0.0, 0.0), (1.0, 1.0), leftdomain=1, rightdomain=2, bc = “scatterer”)
geo.AddCircle ( (0.5, 0.5), r=12.0, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.5, 0.5), r=35.0, leftdomain=3, rightdomain=0,bc = “dir”)
#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.1)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle
#Set the extra value of the scatterer
#Set the extra value of the scatterer
f = 10+2np.pisin(2np.pix/eps)+2np.picos(2np.piy/eps)

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

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

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

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

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

u = GridFunction(fes, multidim=5, name='resonances')

lam = ArnoldiSolver(a.mat, m.mat, fes.FreeDofs(),
u.vecs, shift=tplam)
for l in lam:
if abs(l-(tplam))<0.05:
tlamvec[i]=lam[0]


print(tlamvec)

N=5-1;

sq = np.zeros(N+1);
for j in range(N+1):
sq[j]=np.abs(tlamvec[j]-tplam);
print(sq)

for i in range(N+1):
ax = plt.gca()
ax.scatter(tepsvec[i] ,sq[i] , c=‘blue’);
ax.set_yscale(‘log’);
ax.set_xscale(‘log’);

m, b = np.polyfit(np.log(tepsvec), np.log(sq), 1);

plt.plot(tepsvec, np.exp(b)*tepsvec**m)
plt.title(‘\delta = 0.05, on loglog scale’)
plt.xlabel(‘Epsilon’)
plt.ylabel('Resonance values ')
plt.show()
print(‘The slope is:’, m)

This has slope: [/i][/i][/i][/i][/code][/i][/i][/i][/i]0.9212055445929848