# NIntegrate gives two different answers (errors) for the same integrand

The integrand I want to integrate from u=183/985 to u=5 is:

integrand[u]=(1/u)*Sqrt[u^2-33489/970225]*BesselJ[1, 125*Sqrt[u^2-33489/970225]]*FAu[u^2])


where FAu is an interpolated function of a table of data (more info below). Here is a plot of the integrand from 183/985 to 5: If I integrate this function directly using NIntegrate

NIntegrate[integrand[u],{u,183/985, 5}, WorkingPrecision->10, AccuracyGoal->10]//AbsoluteTiming


I get:

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in u near {u} = {0.533503182102844852657037466919066363800523835958162522877040}. NIntegrate obtained -1.98098626595186890458305013843987581889411200403876177007049*10^-9 and 4.4047371056124174089219217622104219743740461184602648715017360.*^-8 for the integral and error estimates.

{5.87632, -1.980986266*10^-9}

But, if I split the integral in 5: [183/985,1], [1,2], ..., [4,5] and add up the results like this:

Total[NIntegrate[I[u], {u, ##}, WorkingPrecision->10, AccuracyGoal->10]&@@@Partition[Flatten[{183/985,Range[1,5]}],2,1]]//AbsoluteTiming

1. it takes 100x more time to evaluate the integral (565.563s),
2. it gives no error,
3. it gives a different answer (-1.538*10^-11).

Context

To obtain the function FAu[u] I create a table from the integral

TAu[q_]:=(1/104.2591938591825620)NIntegrate[If[q == 0 || r == 0, r^2/(1 + E^((r - 6.64220.)/0.54920.)), (Sin[q*r]/(q*r))*(r^2/(1 + E^((r - 6.64220.)/0.54920.)))], {r, 0, 100}, MinRecursion -> 3, MaxRecursion -> 100, WorkingPrecision -> 15,PrecisionGoal -> 7, AccuracyGoal -> Infinity]


and then interpolate it

FAu=Interpolation[TAu];


Questions

As I have been using Mathematica for numerical calculations for not very long I have a few questions:

1. How can I be sure of the result NIntegrate returns? Since FAu[u] is evaluated using NIntegrate, it itself may contain errors. I've tried several NIntegrate methods and all of them return different results (except GaussKronrodRule, which I think NIntegrate automatically uses it for this case).
2. Why splitting the integral takes much longer to evaluate and why it returns a different answer?
3. I'm still confused with PrecisionGoal and AccuracyGoal options. If I want a result with at least 3 - 4 decimal digits after the period, this means I should be using AccuracyGoal->4 instead of PrecisionGoal?

Thanks

• Copying and pasting your code in Mathematica does not produce valid expressions and some quantities are missing. (E.g. q, TAu.) Please, fix those. Also, you are using WorkingPrecision in a wrong way -- you probably want PrecisionGoal. – Anton Antonov Jun 17 '17 at 1:47
• Recommend reading Mathematica's help on NIntegrate and oscillatory integrands, to include Levin's method. – MikeY Jun 17 '17 at 5:47
• Please upload code that actually runs, including code for generating plots. – MikeY Jun 17 '17 at 13:32
• I'd look into Method -> "InterpolationPointsSubdivision" in the doc center. – Michael E2 Jun 17 '17 at 19:17
• @AntonAntonov I'm very sorry for the mistyping, now it should be correct. In it's definition TAu is a function of 'q', which is then interpolated into FAu[q]. Then, in the integrand, FAu is a function of u^2. I've seen in several places that WorkingPrecision need to be set for PrecisionGoal to work, i.e. WorkingPrecision->wp says that NIntegrate should work internally with precision wp, and returns the answer when PrecisionGoal is satisfied. – Pierre Jun 17 '17 at 19:17

OK, it seems to work fine here. Note the simplification of TAu. The extra options on precision, accuracy, etc., essentially had no effect.

TAu[q_] := (1/104.2591938591825620) NIntegrate[ (Sin[ q*r]/(q*r))*(r^2/(1 + E^((r - 6.642)/0.549))), {r, 0, 100}]

dat = Table[{u, TAu[u]}, {u, (183/985)^2, 5^2+.1, .1}];

FAu = Interpolation[dat];


Define aa=33489/970225 to clean up the integrand.

integrand[u_] :=  (1/u) Sqrt[u^2 - aa] BesselJ[1, 125*Sqrt[u^2 - aa]] FAu[u^2];

NIntegrate[integrand[u], {u, (183/985), 5}]

(*    -5.07583*10^-7    *)

NIntegrate[integrand[u], {u, (183/985), 1}] +
NIntegrate[integrand[u], {u, 1, 2}] +
NIntegrate[integrand[u], {u, 2, 3}] +
NIntegrate[integrand[u], {u, 3, 4}] +
NIntegrate[integrand[u], {u, 4, 5}]

(*    -5.07583*10^-7    *)
`