NIntegrate and MaxRecursion

I’m using NIntegrate to get the value of a very complicated function f[k,T,M], where T and M are some parameters:

S[T_, M_, MaxRec_]:= NIntegrate[f[k, T, M], {k,0,Infinity}, PrecisionGoal->10, MaxRecursion->MaxRec] 

What puzzles me is the fact that, for the same parameters, I get different values depending on MaxRec:

S[10^12, 10^12, 10] 

enter image description here enter image description here

1.30494*10^30 

versus

S[10^12, 10^12, 50] 

enter image description here

-1.161*10^45 

How should I interpret these results? I could provide the explicit form of f if needed. Thanks a lot!

Calculating Residue with NIntegrate fails

Trying to answer this interesting question Principal value from two different axis I observed a problem using Nintegrate:

The function func[p_] := 1/(Sinh[p/2] Sqrt[Cosh[p]]) has a pol at p==0.

The residue of this point evaluates to

Residue[func[z], {z, 0}] (*2*) Limit[func[z] z, z -> 0] (*2*) 

The result might be confirmed by integrating along a path in the complex plane which contains the pol. For example integrating along a square path

NIntegrate[func[z], {z, 1, I, -1, -I, 1}]/(2 Pi I) (*2*) 

evaluates correct value , whereas integrating along a circle

NIntegrate[func[ Exp[I \[CurlyPhi] ]]/(2 Pi I), {\[CurlyPhi], 0, 2 Pi}]  (*~0*) 

gives a message NIntegrate failed to converge... and a wrong result 0!

What’s wrong with this last integration?

How to modify to get the correct result?

Thanks!

NIntegrate with variable in it

I would like to NIntegrate with a variable in the function. Later I will be series expanding it. Can it be done in Matehematica? I am getting errors for a sample integration as,

NIntegrate[Series[Cosh[x]*Exp[-z*Cosh[x]], {z, 0, 2}], {x, 0, \[Infinity]}]

The error is,

NIntegrate::inumr: The integrand Cosh[x]-Cosh[x]^2 z+1/2 Cosh[x]^3 z^2+O[z]^3 has evaluated to non-numerical values for all sampling points in the region with boundaries {{\[Infinity],0.}}.

Is there a way out? I know that, analytically, I can obtain expressions in terms of modified Bessel functions of second kind. But Mathematica does not pick that up.

Instability in NIntegrate

enter image description here I am solving the following equation with NIntegrate. However, the solution branches off at a certain time instance. When I increase the time step size, the instability occurs at larger time instant. Any solution for getting rid of the instability? NB: I have already changed AccuracyGoal and PrecisionGoal but does not work.

Function providing input and integration limits for NIntegrate

I am trying to define some custom function that evaluates the numeric integral of some complicated function. The problem is that I would like to include as an input the integration limits. A MWE follows

f[x_,y_]:= Exp[-2 x^2 y^2] test[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],IntL] myIntL1= {y,-4 x^2,4x^2}; myIntL2= {y,-4 x^4,4x^4}; 

Then, if I evaluate for instance test[3,myIntL1] I get a problem concerning an invalid limit of integration.

Is there a clever way to fix this without defining several functions including the integration limits, such as

    test1[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],{y,-4 x^2,4x^2}]     test2[x_?NumericQ,IntL_]:= NIntegrate[f[x, y],{y,-4 x^4,4x^4}] 

etc?

Of course, here everything looks simple, but in my case all the functions are rather long; some even purely numeric ones. Since I have multiple choices for the integration limits, it would be more practical to avoid defining test1[x_], test2[x_], ..., testN[x_]

Thanks in advance, Pablo

NSolve and NIntegrate, or a better approach

I need to define and plot the following function

$ $ a(t) := \exp\left(\int Z(t)\; dt\right) $ $

where $ Z(t)$ is the solution to the equation

$ $ 0 = t – 2 \int^Z_1 F(x)\; dx $ $

with $ F = F(x)$ being a known (but complicated and non-integrable) function.

How do I define and plot the function $ a(t)$ in Mathematica?

Here is my attempt with a particular function $ F(x)$ that I need to work with:

A = 0; F[x_] = - ((4*A*x^(9/2) + 64*x^6 - 4*Sqrt[A*x^9*(32*x^(3/2) + A)])^(1/3)/(16*x^4 - 4*x^2*(4*A*x^(9/2) + 64*x^6 - 4*Sqrt[A*x^9*(32*x^(3/2) + A)])^(1/3) - (4*A*x^(9/2) + 64*x^6 - 4*Sqrt[A*x^9*(32*x^(3/2) + A)])^(2/3)));  A = 0; Int[Z_?NumericQ] := NIntegrate[F[x], {x, 1, Z}] S[t_?NumericQ] := NSolve[t - 2*Int[Z] == 0, Z] a[t_] := Exp[Integrate[S[t], t]] 

However, when trying to evaluate for example $ a(2)$ I get the following error:

The error

Help with NIntegrate settings to evaluate integral

I am trying to evaluate this integral: \begin{align*} \alpha_{2}=\int_{-\infty}^{1.645}\left[1-\Phi\left(\frac{\sqrt{25}}{\sqrt{15}} 1.645-\frac{\sqrt{10}}{\sqrt{15}} z_{1}\right)\right] \phi\left(z_{1}\right) d z_{1} \end{align*} where $ \Phi$ is the Normal CDF and $ \phi$ is the normal PDF. I know that the answer should be 0.03325.

I used the following code, but it doesn’t converge to an answer. Any suggestions?

pdf[x_] := PDF[NormalDistribution[0, 1], x]  cdf[x_] := CDF[NormalDistribution[0, 1], x]  NIntegrate[ 1 - cdf[Sqrt[25]/Sqrt[15] 1.645 - Sqrt[10]/Sqrt[15] x]* pdf[x], {x, -Infinity, 1.645}]  

which returns

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.  NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {-8.16907*10^224}. NIntegrate obtained 8.582639123060543`15.954589770191005*^27949 and 8.582639123060543`15.954589770191005*^27949 for the integral and error estimates. 

The following code in R gives me the correct answer:

inside <- function(z1, n1, n2, cv) {   nt <- n1 + n2   (1 - pnorm(sqrt(nt/n2) * cv - sqrt(n1/n2) * z1)) * dnorm(z1) }  additional.error <- function(n1, n2, cv) {   integrate(inside, lower = -Inf, upper = cv, n1 = n1, n2 = n2, cv = 1.645)$  value }  additional.error(n1 = 10, n2 = 15, cv = qnorm(0.95)) ``` 

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

enter image description here

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.