I have a square matrix, `m`

which depends on `kx`

and `ky`

. It isn’t Hermitian, but it does have real eigenvalues. I would like to obtain the integral of the sum of the eigenvalues of this matrix over a region in the space of those two variables. After several minutes of evaluation I begin to receive notifications that not all eigenvalues are being found.

`Message[Eigenvalues::eival] `

Presumably this is for only a few points in the region which are probably of vanishing measure.

In contrast to the trouble I have with `NIntegrate`

, when I plot the sum of eigenvalues, `Plot3D`

returns a highly-resolved and smooth plot in a very modest amount of time. Here’s my code:

`omega[k_?VectorQ] := Abs[Chop[Eigenvalues[m /. {kx -> k.{1, 0}, ky -> k.{0, 1}}]]] Plot3D[Total[omega[{kx, ky}]], {kx, ky} \[Element] region, AspectRatio -> Automatic, PlotPoints -> 60, Mesh -> All] NIntegrate[Total[omega[{kx, ky}]], {kx, ky} \[Element] region, PrecisionGoal -> 1] `

I thought I might just use the data in the plot to coarsely approximate the integral myself. So I set `MaxRecursion->0`

and took the average of the points from the plot.

`plot = Plot3D[Total[omega[{kx, ky}]], {kx, ky} \[Element] region, AspectRatio -> Automatic, PlotPoints -> 60, Mesh -> All] dataPoints = Flatten[Cases[plot, x_GraphicsComplex :> First@x, Infinity], 1]; zData = dataPoints[[;; , 3]]; Mean[zData]*Area[region] `

Obviously this is the dumbest way to do the integral, but for some matrices, `m`

, it converges to three decimal places with `PlotPoints->90`

(and with far fewer points in other cases). I think that’s evidence that my integrand isn’t too pathological.

Shouldn’t integration of a smooth function take about as much time as plotting a smooth function? Especially if I can get a sensible integral from such reckless surgery on Plot3D.

So what gives, NIntegrate???

Thanks!