## Nested NIntegrate with FindRoot

I am trying to numerically integrate a function using nested NIntegrate:

$$F(N,x,s)=\int_{-\infty}^s \int_{-\infty}^{+\infty} K(N,z’,x,x’) g_{x’,s’} dx’ds’$$

where the kernel of the integration, $$K(N,z’x,x’)$$, is a messy expression defined in the mathematica code below, and $$g_{x’,s’}$$ is a bi-variate gaussian defined by:

$$g_{x,s’}=\frac{n}{2\pi\sigma_{x’}\sigma_{s’}}\exp\left({ -\frac{x’^2}{2\sigma_{x’}^2} }\right)\exp\left({ -\frac{s’^2}{2\sigma_{s’}^2} }\right).$$

The tricky part(s) is that:

1. $$z’$$ in the $$K(N,z’,x,x’)$$ needs to be solved for numerically using FindRoot and will have a $$s’$$ dependence.
2. The integration upper limit over $$ds’$$ is a variable $$s$$.
3. I suspect the kernel is oscillatory with $$N$$ (denoted "Kernel" in the code below) so maybe an averaging of the kernel over $$N$$ can be done to simplify the kernel and eliminate $$N$$ if the integrations prove to be too time consuming.

At the end, I would like a function, F(N,x,s), that would be able to plot across $$s$$ for a given $$(N,x)$$ values i.e. Plot[F[a,b,s,{s,-1e-5,1e-5}].

(*Constants*) e = -1.60217733*10^-19; m = 9.109389699999999*10^-31; epsilon = 8.854187817620391*10^-12; re = 2.81794092*10^-15; c = 2.99792458*10^8; n = -10^-10/e; KK = 1; lw = 0.026; kw = (2 Pi)/lw; gamma = 4000/0.511; beta = Sqrt[1 - 1/gamma^2]; sigmaS = 10^-5; sigmaX = 30*10^-6; coeff = n/(2 Pi*sigmaS*sigmaX) Exp[-(xprime^2/(2 sigmaX^2))]* Exp[-(sprime^2/(2 sigmaS^2))];  (*Preliminary Equations*) rs2 = {zprime, xprime + KK/(gamma*kw) Sin[kw*zprime], 0}; ro2 = {(NN + 10000)*lw, x + KK/(gamma*kw) Sin[kw*(NN + 10000)*lw], 0};  betas = {beta - KK^2/(2 gamma^2) Cos[kw*zprime]^2,KK/gamma Sin[kw*zprime], 0}; betao = {beta - KK^2/(2 gamma^2) Cos[kw*(NN + 10000)*lw]^2,KK/gamma Sin[kw*(NN + 10000)*lw], 0};  betaDot = {(c*KK^2*kw)/(2 gamma^2)Sin[2 kw*zprime], -((KK*c*kw)/gamma) Sin[kw*zprime], 0};  deltar2 = ro2 - rs2; Rgam2 = Sqrt[deltar2[[1]]^2 + deltar2[[2]]^2];  Ec2 = (e/(4 Pi*epsilon)) (deltar2/Rgam2 - betas)/(gamma^2 Rgam2^2 (1 - (deltar2/Rgam2).betas)^3); Erad2 = (e/(4 Pi*epsilon)) Cross[deltar2/Rgam2,Cross[deltar2/Rgam2 - betas, betaDot]]/(c*Rgam2*(1 - (deltar2/Rgam2).betas)^3);  sumElong = (Ec2[[1]] + Erad2[[1]]); sumEtran = (Ec2[[2]] + Erad2[[2]]);  (*Numerical Functions*)  ZPRIME[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ] := zprime /.FindRoot[s - sprime == (Sqrt[gamma^2 + KK^2] (EllipticE[kw*(NN + 10000)*lw,KK^2/(gamma^2 + KK^2)] - EllipticE[kw zprime, KK^2/(gamma^2 + KK^2)]))/(gamma kw) -beta Sqrt[((NN + 10000)*lw - zprime)^2 + (x - xprime + (KK Sin[kw *(NN + 10000)*lw])/(gamma kw) - (KK Sin[kw zprime])/(gamma kw))^2], {zprime, 0}]   Kernel = coeff re/gamma (sumElong*betao[[1]] + sumEtran*betao[[2]])/.{zprime -> ZPRIME[NN, x, xprime, s, sprime]};  FNxprimesprime[NN_?NumericQ, x_?NumericQ, xprime_?NumericQ, s_?NumericQ, sprime_?NumericQ]:= Kernel  FNsprime[NN_?NumericQ, x_?NumericQ, s_?NumericQ, sprime_?NumericQ] :=NIntegrate[FNxprimesprime[NN, x, xprime, s, sprime], {xprime, -300/10^6, 300/10^6}]  FN[NN_?NumericQ,x_?NumericQ, s_?NumericQ] := NIntegrate[FNsprime[NN,x, s, sprime], {sprime,-10^-4, s}]  lst1 = Table[{ss, FN[0,0, ss], PrecisionGoal -> 5] // Quiet}, {ss, -10^-5, 10^-5, 10^-6}] ListPlot[lst1] 

## Speeding up NIntegrate and DensityPlot

I’m making a density plot of a spectral function Sfull[kx, ky, kz, \[CapitalOmega]] along a specified path in {kx,ky,kz}-space. The plot looks as it should, but it takes a long time to run (12+ hours). Sfull contains numerical integrals which I believe are the source of the problem. The evaluation of Sfull along the {kx,ky,kz}-space path could perhaps be made more efficient as well.

==============

Definition of Sfull[kx_, ky_, kz_, \[CapitalOmega]_]:

\[Theta] = 1.119; Qx = 0; Qy = 0; Qz = 0; \[Alpha] = 0.2; S = 0.5; J = 1.2; \[Gamma][kx_, ky_, kz_] := (Cos[kx] + Cos[ky] + \[Alpha] Cos[kz])/(2 + \[Alpha]); \[Omega]k[kx_, ky_, kz_] := Sqrt[(1 + \[Gamma][kx, ky, kz]) (1-Cos[2\[Theta]] \[Gamma][kx, ky, kz])]; \[Epsilon][kx_, ky_, kz_] := 2 J S (2 + \[Alpha]) \[Omega]k[kx, ky,kz]; Ak[kx_, ky_, kz_] := 2 J S (2 + \[Alpha]) (1 + Sin[\[Theta]]^2\[Gamma][kx, ky, kz]); u[kx_, ky_, kz_] := Sqrt[(Ak[kx, ky, kz] + \[Epsilon][kx, ky, kz])/(2 \[Epsilon][kx, ky,kz])]; v[kx_, ky_, kz_] := Sqrt[(Ak[kx, ky, kz] - \[Epsilon][kx, ky, kz])/(2 \[Epsilon][kx, ky,kz])]; \[CapitalPhi]1[kx_, ky_, kz_, qx_,qy_,qz_] := (\[Gamma] [kx, ky, kz] (u[kx, ky, kz] + v[kx, ky,kz]) (u[qx, qy, qz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[qx, qy, qz] u[kx - qx +Qx, ky - qy + Qy,kz - qz + Qz]) + \[Gamma][qx, qy, qz] (u[qx, qy, qz] + v[qx, qy, qz]) (u[kx, ky, kz] u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx, ky, kz] v[kx - qx + Qx, ky - qy + Qy,kz - qz + Qz]) + \[Gamma][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] (u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) (u[kx, ky, kz] u[qx, qy, qz] + v[kx, ky, kz] v[qx, qy, qz])); \[CapitalPhi]2[kx_, ky_, kz_, qx_,qy_,qz_] := \[Gamma][kx, ky, kz] (u[kx, ky, kz] + v[kx, ky, kz])(u[qx, qy, qz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[qx, qy, qz] u[kx - qx + Qx,ky - qy + Qy,kz - qz + Qz]) + \[Gamma][qx, qy, qz] (u[qx, qy, qz] + v[qx, qy, qz]) (u[kx, ky, kz] v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx, ky, kz] u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) + \[Gamma][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] (u[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + v[kx - qx + Qx, ky - qy + Qy, kz - qz + Qz]) (u[kx, ky, kz] v[qx, qy, qz] + v[kx, ky, kz] u[qx, qy, qz]); \[CapitalSigma][kx_, ky_, kz_, \[CapitalOmega]_]:=(1/2)NIntegrate[(-Abs[\[CapitalPhi]1[kx, ky, kz, qx, qy, qz]]^2/(\[CapitalOmega]-\[Epsilon][qx, qy, qz] - \[Epsilon][kx - qx + Qx, ky - qy + Qy, kz - qz + Qz] + I 0.001) + Abs[\[CapitalPhi]2[kx, ky, kz, qx, qy, qz]]^2/(\[CapitalOmega] + \[Epsilon][qx, qy, qz] + \[Epsilon][kx + qx - Qx, ky + qy - Qy, kz + qz - Qz] - I 0.001)),{qy, 0, 2 \[Pi]}, {qz, 0, 2 \[Pi]}, Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 4, MaxRecursion -> 10^4]; G[kx_, ky_, kz_, \[CapitalOmega]_] := (1/(\[CapitalOmega] - \[Epsilon][kx, ky, kz] + \[CapitalSigma][kx, ky, kz, \[CapitalOmega]])); AA[kx_, ky_, kz_, \[CapitalOmega]_] := ((-1/\[Pi]) Im[G[kx, ky, kz,\[CapitalOmega]]]); Sxx[kx_, ky_, kz_, \[CapitalOmega]_] := (\[Pi] S (u[kx, ky, kz] +v[kx, ky, kz])^2 AA[kx, ky, kz, \[CapitalOmega]]); Syy[kx_, ky_, kz_, \[CapitalOmega]_] := (\[Pi] S (u[kx, ky, kz] -v[kx, ky, kz])^2 AA[kx, ky, kz, \[CapitalOmega]]); Szz[kx_, ky_, kz_, \[CapitalOmega]_] := \[Pi]NIntegrate[((u[qx, qy, qz] v[kx - qx, ky - qy, kz - qz] + v[qx, qy, qz] u[kx - qx, ky -qy,kz - qz])^2 (1/Sqrt[2 \[Pi]] (0.1 J)) Exp[-(\[CapitalOmega] - \[Epsilon][qx,qy, qz] -\[Epsilon][kx - qx, ky - qy, qz - kz])^2/(2 (0.1 J)^2)]), {qx, 0, 2 \[Pi]}, {qy, 0, 2 \[Pi]},{qz,0, 2 \[Pi]}, Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 4, MaxRecursion -> 10^4]; Sxx0[kx_, ky_, kz_, \[CapitalOmega]_] := (Sin[\[Theta]]^2 Sxx[kx,ky,kz, \[CapitalOmega]] + Cos[\[Theta]]^2 Szz[kx - Qx, ky - Qy, kz - Qz, \[CapitalOmega]]); Szz0[kx_, ky_, kz_, \[CapitalOmega]_] := (Cos[\[Theta]]^2 Sxx[kx - Qx, ky - Qy, kz - Qz, \[CapitalOmega]] + Sin[\[Theta]]^2 Szz[kx, ky, kz, \[CapitalOmega]]); Syy0[kx_, ky_, kz_, \[CapitalOmega]_] := (Syy[kx, ky,kz,\[CapitalOmega]]); Sfull[kx_, ky_, kz_, \[CapitalOmega]_] := Sxx0[kx, ky, kz, \[CapitalOmega]] + Syy0[kx, ky, kz, \[CapitalOmega]] + Szz0[kx, ky, kz, \[CapitalOmega]]; 

Definition of {kx,ky,kz}-space path:

kpath = {{0, {0, 0, 0, \[CapitalOmega]}}, {1, {\[Pi], \[Pi], 0, \[CapitalOmega]}}, {2, {\[Pi], 0, 0, \[CapitalOmega]}}, {3, {0, \[Pi], 0, \[CapitalOmega]}}, {4, {0, 0, 0, \[CapitalOmega]}}}; if = Interpolation[kpath, InterpolationOrder -> 1];  

Density plot:

DensityPlot[{Sfull@@if[i]}, {i, 0, 4}, {\[CapitalOmega], 0, 10},ColorFunction -> "SunsetColors", PlotLegends -> Automatic] 

This produces the plot

which takes over 12 hours to run. \[CapitalSigma][kx, ky, kz, \[CapitalOmega]] and Szz[kx, ky, kz, \[CapitalOmega]] contain the numerical integrals (which I would like to perform according to the adaptive quasi Monte Carlo scheme with AccuracyGoal -> 4, MaxRecursion -> 10^4) and I’d like to speed these up. I’m also unsure that using Sfull@@if[i] is the most efficient way to produce the density plot. I would eventually like to make the plot smoother by increasing PlotPoints, but this would of course increase the runtime further.

I have tried SymbolicPreprocessing->0 and defining functions with ?NumericQ variables but this does not speed up the integrals noticably.

## Discrepancy between the results of NIntegrate with different methods and options

I am trying to perform a numerical integration on a function defined through a sum of exponential terms. The summation is given by:

sum[z_, z0_, t_, nmax_] :=   1/Sqrt[4 t]*Sum[ Exp[-(z - z0 - 2 n)^2/(4 t)] + Exp[-(z + z0 - 2 n)^2/(4 t) ], {n, -nmax, nmax}]; 

and we define

f[z_, zp_, y_, yp_, z0_, t_] =    ( Exp[-(y - yp)^2/(8t)]/Sqrt[8t] )*D[sum[z, z0, t, 20]*sum[zp, z0, t, 20], t, z, zp]; 

where I have chosen nmax=20.

I wish to perform numerical integration on f. I define

int1[a_] :=   NIntegrate[ f[zp, zpp, 0, ypp, z0, t]*( ypp (zp - zpp) )/( a (zp - zpp)^2 + ypp^2 )^(3/2),  {z0, 0., 1.}, {t, 0., 10.}, {zp, 0., 1.}, {zpp, 0., 1.}, {ypp, 0., Infinity}] 

Based on which I get the following table (no errors generated)

tab1 = Table[{a, int[a]}, {a, 1., 5., .5}]  (* {{1., 0.00135643}, {1.5, 0.000734155}, {2., -0.000611633},  {2.5, 0.0000596739}, {3.,0.0359735}, {3.5, 0.0292143},  {4., 0.01122}, {4.5, 0.00889722}, {5., 0.00649666}} *) 

To check, I tried AccuracyGoal-> 30 and PrecisionGoal -> 30 and got the same results. However, as soon as I include WorkingPrecision, the integrated gives 0 all the time. For example:

intwrk[a_] :=   NIntegrate[ f[zp, zpp, 0, ypp, z0, t]*(ypp (zp - zpp))/(a (zp - zpp)^2 + ypp^2)^(3/2),  {z0, 0, 1}, {t, 0, 10}, {zp, 0, 1}, {zpp, 0, 1}, {ypp, 0, Infinity},  AccuracyGoal -> 30, PrecisionGoal -> 30, WorkingPrecision -> 100] 

Gives all zeros – which I can’t seem to understand. Am I doing something wrong here?

Then I also performed the calculation using LocalAdaptive method:

int2[a_] :=   NIntegrate[ f[zp, zpp, 0, ypp, z0, t]*(ypp (zp - zpp))/(a (zp - zpp)^2 + ypp^2)^(3/2), {z0, 0., 1.}, {t,0., 10.}, {zp, 0, 1}, {zpp, 0., 1.}, {ypp, 0., Infinity }, Method -> "LocalAdaptive"] 

which gives

tab2 = Table[{a, int2[a]}, {a, 1., 5., .5}]  (*  {{1., 1.75934*10^-29}, {1.5, -8.79671*10^-30}, {2., 8.7967*10^-30},  {2.5, -8.7967*10^-30}, {3., 8.7967*10^-30}, {3.5, -2.80635*10^-51},  {4., -8.79671*10^-30}, {4.5, -9.78474*10^-37}, {5., -9.78474*10^-37}} *) 

which is a very different result. This one also runs much faster than the int1. I guess these numbers here are not reliable, but how can I check them?

PS: based on the answer to my previous post, I tried Method -> "GaussKronrodRule" and Method -> {"MultidimensionalRule", "Generators" -> 9}, but they run forever and I couldn’t get any outcome.

## Forcing NIntegrate not to store full solution, in order to save memory

When using NIntegrate to integrate a time-varying dynamical system, I often only care about the final value of the dependent variables at the end of integration. However, NIntegrate returns interpolation functions, implying that it is storing all intermediate values of the dynamical variables in memory.

Is there a way to save memory by forcing NIntegrate to only return the final values of the dependent variables?

## NIntegrate does not evaluate this finite integral composed of divergent parts

I would like to numerically evaluate the following integral:

$$I = \int_{-\infty}^\infty d\tau_3 \int_{-\infty}^\infty d\tau_4 \frac{1}{1+\tau_3^2} \left\lbrace \frac{2}{1+\tau_4^2} \log (\tau_3 – \tau_4)^2 + \left(\frac{1}{1+\tau_3^2} + \frac{1}{1+\tau_4^2} \right) \phi(\tau_3,\tau_4) \right\rbrace \tag{1}$$

with $$\phi(r,s)$$ a complicated function as defined in the code below. Note that the first term with the log is divergent, but that this divergence is canceled by another divergence present in the 2nd term with the $$\phi$$-function. When I try to evaluate the integral, NIntegrate stays unevaluated. Why is that, and what is the numerical value of this integral?

Here is the code I used so far:

S[\[Tau]3_, \[Tau]4_] := (\[Tau]3 - \[Tau]4)^2/(1 + \[Tau]3^2); a[\[Tau]3_, \[Tau]4_] := 1/4 Sqrt[4*R[\[Tau]3, \[Tau]4]*S[\[Tau]3, \[Tau]4] - (1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4])^2] ;  F[\[Tau]3_, \[Tau]4_] := I Sqrt[-((1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] - 4 I*a[\[Tau]3, \[Tau]4])/(1 - R[\[Tau]3, \[Tau]4] - S[\[Tau]3, \[Tau]4] + 4 I*a[\[Tau]3, \[Tau]4]))]; phi[\[Tau]3_, \[Tau]4_] := 1/a[\[Tau]3, \[Tau]4] Im[PolyLog[2, F[\[Tau]3, \[Tau]4]*Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]] + Log[Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]] Log[1 - F[\[Tau]3, \[Tau]4]*Sqrt[R[\[Tau]3, \[Tau]4]/S[\[Tau]3, \[Tau]4]]]]; NIntegrate[1/(1^2 + \[Tau]3^2) (2/(1^2 + \[Tau]4^2)Log[(\[Tau]3 - \[Tau]4)^2] + (1/(1^2 + \[Tau]3^2) + 1/(1^2 + \[Tau]4^2)) phi[\[Tau]3, \[Tau]4]), {\[Tau]3, -\[Infinity], \[Infinity]}, {\[Tau]4, -\[Infinity], \[Infinity]}] 

## Understanding why NIntegrate requires explicit substitution of variables in argument

The problem that leads me here begins with a quantity I have previously defined, let’s call it test, that has many other quantities in its definition. When evaluated, test is an expression that includes two variables, let’s call them k and x. A MWE would be

test = kx; 

I wish to create a function of k that includes a NIntegration of test over x with limits that involve k. An example would be

testint[k_] := NIntegrate[test, {x, k, 2k}]; 

Evaluating this for some k , say k = .1, returns an error:

>>testint[.1]  NIntegrate::inumr: The integrand k x has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.1,0.2}}.  NIntegrate[k x, {x, 0.1, 0.2}] 

However, if I define testint using a temporary variable and perform a replacement in the argument of NIntegrate , then it computes fine:

>>testint[k1_] := NIntegrate[test/.k->k1, {x, k1, 2k1}]; >>testint[.1] 0.0015 

I found this answer, which led me to try the explicit substitution: Replace variable with value prior to evaluating NIntegrate
Another answer addresses the order of NIntegrate with the help of the ?NumericQ pattern check: How do I prevent NIntegrate::inumr errors within other functions?

My question is Why does NIntegrate require this explicit substitution in order to compute?

As a test, I even tried removing NIntegrates HoldAll attribute thinking that would force the evaluation of test before the integration. It did, but not soon enough to help.

>>test = k x; >>ClearAttributes[NIntegrate, HoldAll] >>testint[k_] := NIntegrate[testin, {x, k, 2 k}]; >>testint[.1]//Trace  NIntegrate::inumr: The integrand k x has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.1,0.2}}.  {testint[0.1], NIntegrate[test, {x, 0.1, 2 0.1}], {test, k x}, {{2 0.1, 0.2}, {x, 0.1, 0.2}}, NIntegrate[k x, {x, 0.1, 0.2}], {{x} =., {x =.}, {x =., Null}, {Null}}, {x =., Null}, ... 

## NIntegrate user defined function

I have the the function

f[a_] := Module[{solution, ans, x}, solution = NSolve[x + a == 4]; 

I want to integrate it by writing

NIntegrate[f[x],{x,1,5} 

but I get this error message:

NIntegrate::inumr: The integrand x\$  2682 has evaluated to non-numerical values  for all sampling points in the region with boundaries {{1,5}}. 

In my problem the function is a bit more complicated but the principle is the same. Is there a nice way to integrate user defined equation-solving functions using for instance NIntegrate? I am aware that you can use f[x] to generate a table, interpolate the table, and finally integrate. However, is there perhaps a nicer way of doing this that avoids creating a interpolating function?

## FindRoot with NIntegrate: either recursion depth exceeded or numerical integration converging too slowly

Trying to numerically solve an equation involving an integral, either the error Recursion depth exceeded appears (when AccuracyGoal->10 or smaller) or the error Numerical integration converging too slowly appears. In the MWE below, the functions involved are smooth, bounded and the solution for the equation is 1 (as confirmed by Plot[{hdif[z], 0}, {z, smin, 2}]), so there should not be difficulty numerically integrating or solving.

How to find a solution to similar equations which cannot be solved by hand analytically?

MWE

Clear[smin, smax, h, hdif] smin = -1; h[z_] := Exp[-z/2] hdif[zmax_?NumericQ] :=   NIntegrate[(Exp[z] - 1)*h[z], {z, smin, zmax}, AccuracyGoal -> 5] smax = FindRoot[{hdif[smax] == 0}, {smax, 0, smin, Max[-smin, 1000]}] 

Other solutions recommend increasing AccuracyGoal or WorkingPrecision (both of which I tried), cancelling periodic oscillations, integrating piecewise around points where the derivative is discontinuous (these do not apply to the function in question). Replacing NIntegrate with Integrate works for h[z_] := Exp[-z/2], but not for functions without a closed form integral.

## NIntegrate problems

Sorry, I am boxed in by Mathematica again, every exit being blocked and can’t find any answers on the web. This usually happens once the code is lengthy so it is hard to debug and to show with the questions. Sorry

1) why can’t I define a function that calls Nintegrate? It seems to want to immediately integrate with values it has at hand. Having to put Nintegrate directly inside my table generation makes the code very hard to read and to change. Neither F[k_] or F[k:_] works. I don’t understand the difference.

2) when I tried to swap in a new NIntegrate based on a manual integration by parts, it is more stable than the old in single integrations, but when swapped in the table entries for it have { curly brackets around them} and don’t plot.

## Fast Plot3D, failing NIntegrate, and reckless surgery

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!