4
$\begingroup$

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: Plot of the integrand 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.40473710561241740892192176221042197437404611846026487150173`60.*^-8 for the integral and error estimates.

{5.87632, -1.980986266*10^-9}

as my answer.

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.25919385918256`20)NIntegrate[If[q == 0 || r == 0, r^2/(1 + E^((r - 6.642`20.)/0.549`20.)), (Sin[q*r]/(q*r))*(r^2/(1 + E^((r - 6.642`20.)/0.549`20.)))], {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

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Anton Antonov Jun 17 '17 at 1:47
  • $\begingroup$ Recommend reading Mathematica's help on NIntegrate and oscillatory integrands, to include Levin's method. $\endgroup$ – MikeY Jun 17 '17 at 5:47
  • $\begingroup$ Please upload code that actually runs, including code for generating plots. $\endgroup$ – MikeY Jun 17 '17 at 13:32
  • 1
    $\begingroup$ I'd look into Method -> "InterpolationPointsSubdivision" in the doc center. $\endgroup$ – Michael E2 Jun 17 '17 at 19:17
  • $\begingroup$ @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. $\endgroup$ – Pierre Jun 17 '17 at 19:17
1
$\begingroup$

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.25919385918256`20) 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    *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.