Question about Integrate

I apologize if this seems naive but I'm running the following code but every time I rerun the code the value of the integral changes. On occasions the value of the integral has positive real part, on others it has negative real parts, etc. I have the the out resonance of the Arnoldi unchanged through each run but the u still changes. 

Could someone pinpoint what is occurring to cause this lack of consistency? 

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.AddCircle( (0.0, 0.0), r=1.0, leftdomain=1, rightdomain=2, bc = "scatterer")
geo.AddCircle ( (0.0, 0.0), r=2.0, leftdomain=2, rightdomain=3)
geo.AddCircle ( (0.0, 0.0), r=3.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.04)
mesh = Mesh(ngmesh)

#Tell the program that we are using PML at the second circle
mesh.SetPML(pml.Radial(rad=2.0, alpha=1j, origin=(0,0)), "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.111)
Draw(u)

print ("lam: ", lam)

# print(Integrate(u, mesh))
print(Integrate(u, mesh))

The Arnold method starts with a random initial guess, which leads to random scaling of eigenvectors