2
$\begingroup$

I have a function defined as

rpd[r_, OptionsPattern[]] :=
  Module[
   {A = OptionValue[A]},
   NIntegrate[ 
    q/(q^2 + r) E^(-q^2/2)*
     Cos[3*q] BesselJ[1, A*q], {q, 0, \[Infinity]}]];

Options[rpd] = {A -> 1};

and plotting it for r \Elem {0,0.001} I see it crosses zero for some value of r (call it r0) which depends on the parameter a<1, like this:

Plot of the function

I'd like to plot a graph of r0(a), so I'm trying to produce a table

rpdA = Table[{a, FindRoot[rpd[r, A -> a], {r, 0.000001}]}, {a, 0.1, 1,
     0.1}];

but I keep getting instead a lot of error messages, namely NIntegrate::ncvb and NIntegrate::slwcon. Could you help me sort this out?

$\endgroup$
0
0
$\begingroup$

I guess this works:

  • Subdivide the interval at 10 to separate the significant oscillatory part from the superexponential decay part.
  • Increase WorkingPrecision to handle the round-off error from the oscillatory part
  • Use the secant method in FindRoot to prevent bad choices for r in trying to numerically approximate the gradient.
  • ?NumericQ protection for rpd[].

ClearAll[rpd];
rpd[r0_?NumericQ, OptionsPattern[]] := 
  Module[{A = SetPrecision[OptionValue[A], 32], 
    r = SetPrecision[r0, 32]}, 
   NIntegrate[
     q/(q^2 + r) E^(-q^2/2)*Cos[3*q] BesselJ[1, A*q], {q, 0, 
      10, \[Infinity]}, MaxRecursion -> 20, WorkingPrecision -> 32, 
     PrecisionGoal -> 6] /; NumericQ[A]];

Options[rpd] = {A -> 1};

rpdA = Table[{a, 
    FindRoot[rpd[r, A -> a], {r, 0.000001, 0.0000001}]}, {a, 0.1, 1, 
    0.1}];

rpdA
(*
{{0.1, {r -> 0.0000846907}}, {0.2, {r -> 0.0000899893}},
 {0.3, {r -> 0.0000993154}}, {0.4, {r -> 0.000113453}},
 {0.5, {r -> 0.000133586}}, {0.6, {r -> 0.000161389}},
 {0.7, {r -> 0.000199148}}, {0.8, {r -> 0.000249928}},
 {0.9, {r -> 0.000317773}}, {1., {r -> 0.000407975}}}
*)
$\endgroup$
1
  • $\begingroup$ I just add for future readers that pdA = Table[{a, Values@@FindRoot[rpd[r, A -> a], {r, 0.000001, 0.0000001}]}, {a, 0.1, 1, 0.1}]; produces precisely the table I was looking for, which can be plotted with ListPlot. $\endgroup$ – Davide Venturelli Feb 16 at 12:29

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.