How to deal with the loss of significant digits in this expression with Fresnel integrals?

Question

I needed to solve this integral:

$$\int_1^\infty \frac{dx}{\sqrt{x}} \cos \left(a x-\frac{\pi}{4} \right) \cos \left(b x-\frac{\pi}{4} \right) \cos \left(c x-\frac{\pi}{4} \right)$$

Coming from approximating the corresponding integral with Bessel functions.

The only way I saw is to expand the product into a sum of trig functions and use Fresnel integrals. With a little help from Mathematica, I've got:

$$\cos \left(a x-\frac{\pi}{4} \right) \cos \left(b x-\frac{\pi}{4} \right) \cos \left(c x-\frac{\pi}{4} \right)= \\ = \frac{1}{4} \left(\cos \left(w_1 x+\frac{\pi}{4} \right)+\cos \left(w_2 x+\frac{\pi}{4} \right)+\cos \left(w_3 x+\frac{\pi}{4} \right)-\cos \left(w_4 x+\frac{\pi}{4} \right) \right)$$

Where: $$w_1=a-b-c \\ w_2=b-a-c \\ w_3=c-a-b \\ w_4=a+b+c$$

Here I use the Mathematica definition for Fresnel integrals which is different from the usual one.

Taking the integrals (again, with the help of Mathematica) I get:

$$\int_1^\infty \frac{dx}{\sqrt{x}} \cos \left(w x+\frac{\pi}{4} \right)= \\ =\frac{\sqrt{\pi}}{2 \sqrt{|w|}} \left(1-\operatorname{sign} w+2\operatorname{sign} w~ S \left(\frac{ \sqrt{2|w|}}{\sqrt{\pi}} \right)-2 C \left(\frac{ \sqrt{2|w|}}{\sqrt{\pi}} \right) \right)$$

This indeed agrees with the numerical integral. However, only when I set the precision of the arguments $>100$.

Can I transform the exact expression in such a way that the cancellation of the digits and the loss of precision doesn't occur? Or at least is minimal?

Here's the code (the coefficient in front of the result is slightly different than in the formula above), and the resul showing the loss of digits for WorkingPrecision->50:

a = RandomReal[{1, 100}, WorkingPrecision -> 50];
b = RandomReal[{1, 100}, WorkingPrecision -> 50];
c = RandomReal[{1, 100}, WorkingPrecision -> 50];
w1 = a - b - c;
w2 = b - a - c;
w3 = c - a - b;
w4 = a + b + c;
F[w_] := Module[{v}, v = N[Sqrt[2/\[Pi]] Sqrt[Abs[w]], 50];
   N[(1 - Sign[w] - 2 FresnelC[v] + 2 Sign[w] FresnelS[v])/Sqrt[
     Abs[w]], 50]];
N[{a, b, c}, 10]
F[w1] + F[w2] + F[w3] - F[w4]

Output:

{33.34824876, 95.28687965, 41.03797714}

-0.011752

For smaller precision I often don't get any significant digits.

I would really like to use this integral in an algorithm where keeping this huge precision is not very convenient, so if you have any ideas, I would be grateful.

Answer 1

It would be helpful if you showed all code used.

Here is what I see from the integral. Without loss of generality (due to the symmetry) I arbitrarily declare a<=b<=c. Not sure if it helped.

ii = 
 Integrate[
  1/Sqrt[x]*Cos[a*x - Pi/4]*Cos[b*x - Pi/4]*Cos[c*x - Pi/4], {x, 1, 
   Infinity}, Assumptions -> 1 <= a <= b <= c]

Let's do a sanity check.

SeedRandom[1234]
replacements = Thread[{a, b, c} -> Sort[RandomReal[{1, 100}, 3]]];
ii /. replacements
NIntegrate[
 1/Sqrt[x]*Cos[a*x - Pi/4]*Cos[b*x - Pi/4]*Cos[c*x - Pi/4] /. 
  replacements, {x, 1, Infinity}]

(* Out[119]= -0.00656369548086 - 3.84339663156*10^-18 I

During evaluation of In[117]:= NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {4.8103669517*10^7}. NIntegrate obtained -0.00651653024521 and 0.0001552783984406559` for the integral and error estimates.

Out[120]= -0.00651653024521 *)

One can get this in terms of Fresnels using the indefinite integral on 1TrigExpanded input (maybe also the definite integral, I didn't have the patience for it to finish).

indefii = 
 Integrate[
  1/Sqrt[x]*TrigExpand[Cos[a*x - Pi/4]*Cos[b*x - Pi/4]]*
   Cos[c*x - Pi/4], x]
hi = Limit[indefii /. replacements, x -> Infinity];
lo = indefii /. replacements /. x -> 1;
hi - lo

(* Out[130]= 1/4 Sqrt[\[Pi]] (FresnelC[
    Sqrt[a - b - c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a - b - c] + 
   FresnelC[Sqrt[a + b - c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a + b - c] + 
   FresnelC[Sqrt[a - b + c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a - b + c] - 
   FresnelC[Sqrt[a + b + c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a + b + c] - 
   FresnelS[Sqrt[a - b - c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a - b - c] + 
   FresnelS[Sqrt[a + b - c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a + b - c] + 
   FresnelS[Sqrt[a - b + c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a - b + c] + 
   FresnelS[Sqrt[a + b + c] Sqrt[2/\[Pi]] Sqrt[x]]/Sqrt[a + b + c])

Out[133]= -0.00656369548086 + 0. I *)

So we are getting good agreement with quadrature, and in all cases using machine arithmetic.

Answer 2

I have solved this problem, it seems, by using the connection of Fresnel integrals to the error function.

Have to note, that Mathematica using different argument convenctions than the commonly accepted ones is really inconvenient.

Here's the code with seemingly no loss of significance:

F[w_] := Module[{v, Ret}, v = N[ Sqrt[Abs[w]]/Sqrt[2], 20];
   If[w > 0, 
    Ret = N[I (Erf[(1 + I) v] - Erf[(1 - I) v])/Sqrt[Abs[w]], 20], 
    Ret = N[(2 - (Erf[(1 + I) v] + Erf[(1 - I) v]))/Sqrt[Abs[w]], 20]];
   Re[Ret]];

If anyone has any further advice, I would be grateful.