NDSolve value precision

Question

I have a function of a numerical integral that works perfectly. Because I found out that I can get the same result using NDSolve in a much faster way (around 0.15s using the integral and 0.02s using NDSolve, per point), I translated this integral (and, in the future, other problematic integrals) into an ODE. The thing is that from some values it returns a lot of errors.

The integral is

FAu[q_?NumericQ] := 1/92.57237188637842`20 NIntegrate[
If[q == 0 \[Or] r == 0, r^2/(1 + E^((r - 6.38`20)/0.535`20)), 
Sin[q r ]/q  r/(1 + E^((r - 6.38`20)/0.535`20))], {r, 0, 100}, 
MinRecursion -> 3, MaxRecursion -> 100, WorkingPrecision -> 15, 
PrecisionGoal -> 7, AccuracyGoal -> \[Infinity]];

and the corresponding NDSolve version is

FAu1[q_] := If[q == 0, 1, NDSolveValue[{y'[r] == 
1/92.57237188637842 Sin[q r]/q r/(1 + E^((r - 6.38)/0.535)), 
y[0] == 0}, y[100], {r, 0, 100}]]

Here is a plot to compare the functions:

T = ParallelTable[{q, Abs@FAu[q/0.197]}, {q, 0, 2, 1/1000}];
A = ListLogPlot[T, PlotStyle -> Red]
B = LogPlot[Abs@FAu1[q/0.197], {q, 0, 2}]
Show[A, B]

The ODE versions is so much faster that Plot is able to create the graph in a small amount of time, while the numerical integral version need to be plotted as a list. The blue line corresponds to the ODE, and, as you can see, for values q>1 it starts giving a lot of errors. I would like to find out why and how to solve this problem. The things I tried to do is to add WorkingPrecision->20 (and adding SetPrecision[...,20] around the equations) all that does is pushing the errors a little bit farther. I also tried adding AccuracyGoal->Infinity & PrecisionGoal->15 at the same time but for some reason I get a 'ComplexInfinity' error. Ideally, I would like to be able to get a good precision for values of q at least in [0,2000].

Just to be clear, I do not wish to rewrite the integral version in a more efficient way, I want to focus on the ODE.

Thanks

Answer 1

You're on the right track with AccuracyGoal, which by default is about 10^-8. Not coincidentally, as y[r] approaches 10^-8, noise starts to creep in. You should think of AccuracyGoal as setting a tolerance such that within that tolerance, numbers should be treated as roughly equivalent to zero. But because of rounding error, it should not be so great (such as Infinity) that it is impossible to achieve at the chose working precision. So AccuracyGoal needs to be great enough that the accuracy throughout the interval of integration is great enough to get enough digits of precision for a satisfactory result. At machine precision, setting it to ~16 digits more than -Log10 of the "average" value of y[r] over a "substantial" interval (e.g., ignore zero crossings). So

AccuracyGoal -> 24

seems a reasonable choice. Over a lengthy interval like this, error can still accumulate, and some slight fuzziness can be seen at the end. I've had success sometimes increasing the PrecisionGoal to a little more than half working precision, (Half working precision is the standard default and is a somewhat conservative heuristic.) For machine precision (~15.95), PrecisionGoal -> 10 is often achievable, but usually not one higher than 12.

So with that in mind, this seems to get the expected and sought-after result:

FAu1[q_] := If[q == 0,
   1,
   NDSolveValue[
    {y'[r] == 1/92.57237188637842 Sin[q r]/q r/(1 + E^((r - 6.38)/0.535)), 
      y[0] == 0}, y[100], {r, 0, 100}, 
    AccuracyGoal -> 24, PrecisionGoal -> 10]];

Show[A, B]