Stability issues a epsilon goes to zero for resonances

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 mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), "PML") #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) a += grad(u)*grad(v)*dx #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 a.Assemble() m.Assemble() u = GridFunction(fes, multidim=1, name='resonances') with TaskManager():     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)

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     mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), "PML") #Set the extra value of the scatterer #Set the extra value of the scatterer     f = 10+2*np.pi*sin(2*np.pi*x/eps)+2*np.pi*cos(2*np.pi*y/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)     a += grad(u)*grad(v)*dx #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 += f*u*v*dx('inner') #convert bilinearform into a matrix form     a.Assemble()     m.Assemble()     u = GridFunction(fes, multidim=5, name='resonances')     with TaskManager():         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]

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]

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     mesh.SetPML(pml.Radial(rad=12.0, alpha=1j, origin=(0.5,0.5)), "PML") #Set the extra value of the scatterer #Set the extra value of the scatterer     f = 10+2*np.pi*sin(2*np.pi*x/eps)+2*np.pi*cos(2*np.pi*y/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)     a += grad(u)*grad(v)*dx #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 += f*u*v*dx('inner') #convert bilinearform into a matrix form     a.Assemble()     m.Assemble()     u = GridFunction(fes, multidim=5, name='resonances')     with TaskManager():         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]