NIntegrate fails while Integrate works

Question

I have a function $f(t)$ defined as

$f(t)=\int\limits_0^t(t-\xi)^{\alpha-1}\ \cos(\xi)\ d\xi$

where $0<\alpha<1$. I now want to evaluate this integral at various values of time. Therefore, my code for the function definition reads:

f[t_] := NIntegrate[(t-x)^(a-1) Cos[x], {x, 0, t}]

For the sake of this example, let us choose $\alpha=.3$, and let us attempt to evaluate $f(t)$ at $t=0.2$. Calling f[.2] returns the following error:

NIntegrate::zeroregion: Integration region {{0.2,0.200000000000000011102230246246115891736237825076626353161415959706}} cannot be further subdivided at the specified working precision. NIntegrate assumes zero integral there and on any further indivisible regions.

NIntegrate::inumri: The integrand Cos[x]/(0.2 -x)^0.7 has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.200000000000000011102230246246115891736237825076626353161415959706,0.200000000000000011092200115363618871738896620143270802392916628941}}.

What is most puzzling is that there exists an analytical solution this integral in terms of the hyper-geometric function. Defining $f(t)$ as

f[t_] := Integrate[(t-x)^(a-1) cos[x], {x, 0, t}]

and then calling f[.2] returns the correct value of the integral with no errors: $f(.2)=2.02934$. In fact, the analytical solution is

f[t] = ConditionalExpression[3.33333 t^0.3 HypergeometricPFQ[{1}, {1.15, 0.65}, -(t^2/4)], 
    Re[t] > 0 && Im[t] == 0]

Does this error has something to do with WorkingPrecision? I am relatively new to Mathematica, so forgive me if this is trivial.

UPDATE

My ultimate goal is to make this function Listable so that I can pass a list of times as the function argument. Consider again the function defined in terms of Integrate (and not NIntegrate). Here $\alpha-0.1$. I now define an array to form my list of times.

time = Array[#/5 - .2 &, 11];

Passing this to $f(t)$ as f[time] gives me the error

Power::infy: "Infinite expression 1/0^1.9 encountered."

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

I note that it does give me the output, despite the error. However, if I pass the list of times manually I get the desired output:

f[{0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2}]

{0, 8.36644, 8.50035, 8.06263, 7.20356, 6.00083, 4.52232, 2.83781, 1.02105, -0.851216, -2.70136}

Is this also related to WorkingPrecision and how Array defines its entries?

Answer 1

f[t_] := NIntegrate[(t - x)^(0.3 - 1) Cos[x], {x, 0, t}, WorkingPrecision -> 10]

f[0.2`20]

Gives:

2.029340978

For the listable situation, it is due to WorkingPrecision as before. You can do the following:

time = Array[#/5 - .2`20 &, 11];

Then Map f across as follows:

f /@ time

Which gives:

{0, 2.029340978, 2.398210571, 2.524050347, 2.479729745, 2.294755393,
1.989951583, 1.585123960, 1.101227368, 0.5608006799, -0.01230648010}

Answer 2

Introduction

The question exposes an interesting issue in writing code for NIntegrate, and I wish to discuss the more general issues, rather than focusing on the localized problem of how best to deal with OP's particular integral. I've encountered this issue before with singular integrals. While I won't claim it is a common problem, the more general issue underlying the OP's code deserves to have an explanation published somewhere.

Analysis of the problem

If you look closely at the inputs reported in the error message, you will notice that the interval appears to be outside the interval of integration. One can also see that NIntegrate has extended the precision of the input by 50 digits. My guess is that some sort of rounding led to the slight increase in the interval, such as the following:

SetPrecision[0.2, 50 + $MachinePrecision]
(*  0.200000000000000011102230246251565404236316680908203125000000000000  *)

It turns out this is not a problem in itself, but it is a clue that figures in understanding the second message.

The second error arises in the definition of f in that f[0.2] constructs an integrand with machine-precision parameters. These do not appear to be promoted to high precision, when NIntegrate goes to high-precision abscissae. This seems reasonable when one considers what sort of arbitrary expressions can be used as integrands and that the WorkingPrecision is MachinePrecision, so I would not call it a bug. But it is a pitfall likely to be overlooked by most users, and probably deserves to be better known. It leads to the second error because (t - x) rounds to zero at machine precision and the exponent (a - 1) is negative, when t = 0.2 and x is in the interval reported in the error message; here are the values at the end points:

(0.2` - x) /. x ->
 {0.200000000000000011102230246246115891736237825076626353161415959706`65.954589770191,   
  0.200000000000000011092200115363618871738896620143270802392916628941`65.954589770191}
(*  {0., 0.}  *)

Solution

This suggests that the precision of the parameters in the integrand need to be larger to avoid the catastrophic cancellation observed in (t - x). Here's a way:

ff[a_, t_] := With[{tt = SetPrecision[t, Infinity]}, 
  NIntegrate[(tt - x)^(a - 1) Cos[x], {x, 0, tt}]
  ];

ff[0.3, 0.2]
(*  2.02934  *)

It should be noted that it is not simply a question of raising WorkingPrecision. Because of the nature of the singularity in the case a = 0.3 and how NIntegrate attacks it, the precision of t has to be much higher than the working precision. (I picked Infinity above so I wouldn't have to worry if it was high enough, even if I want to adjust the integral later.) For instance, define f with WorkingPrecision -> 20:

f[a_, t_] := NIntegrate[(t - x)^(a - 1) Cos[x], {x, 0, t}, WorkingPrecision -> 20]

Then we get the following results:

f[0.3`20, 0.2`20] (* gives BOTH errors as in the OP *)
f[0.3`30, 0.2`30] (* gives ONE error, NIntegrate::ncvb, "failed to converge" *)
f[0.3`40, 0.2`40] (* gives NO errors, ans = 2.0293415670719667237 *)

To make WorkingPrecision -> 20 work, one needs to keep the PrecisionGoal low, by setting PrecisionGoal -> 8.

Conclusion

In general, to construct effective functions for numerical solvers in Mathematica is a bit of an art, or at least requires some foresight and knowledge of the workings of the solver. Consider the question of whether to use exact or approximate numbers (and for approximate numbers, whether to use machine precision or arbitrary precision). Machine precision numbers are generally calculated the fastest, and exact numbers tend to lead to large expressions, potentially using large amounts of time and memory. Thus on the one hand, the standard advice is to use machine precision numbers as soon as possible, although it silently assumes an accurate value will result. On the other hand, we can see that the assumption is not met in this use case. My own approach is to hold on to exact expressions up until numerical computation begins or other operations that are slow on symbolic input. It makes raising WorkingPrecision easier. It makes symbolic analysis easier when needed.

Postscript

For the particular problem at hand, integrating by parts removes the singularity:

ClearAll[fff];
SetAttributes[fff, Listable];
fff[a_, t_] := 
 First@Differences[-(t - x)^(a)/a Cos[x] /. {{x -> 0.}, {x -> t}}] - 
  NIntegrate[(t - x)^(a)/a Sin[x], {x, 0, t}];

fff[0.3, 0.2]
(*  2.02934  *)

And why not put the simplest approach in an answer, from J.M.'s comment:

hardly[a_, t_] := t^a HypergeometricPFQ[{1}, {(a + 1)/2, a/2 +1}, -t^2/4]/a

Answer 3

Instead of increasing WorkingPrecision you can just adjust the integration contour:

f[t_] := NIntegrate[(t-x)^(.3-1) Cos[x], {x, 0, t(1+I), t}]

For the OP example:

f /@ {0, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2} //Chop

{0, 2.02934, 2.39821, 2.52405, 2.47973, 2.29476, 1.98995, 1.58512, 1.10123, 0.5608, -0.012307}