# Find root of function defined via NIntegrate

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: 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?

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}}}
*)

• 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. – Davide Venturelli Feb 16 at 12:29