0
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Just for note, the definition for v looks strange, it equals to $v = \sqrt{\frac{2|w|}{\pi}}$, which doesn't look right. And dividing only by $\sqrt{|w|}$. $\endgroup$
    – m0nhawk
    Nov 17, 2018 at 22:55
  • $\begingroup$ Thanks. I didn't know about it, just noted that the equation and Mathematica code are different. $\endgroup$
    – m0nhawk
    Nov 17, 2018 at 23:09
  • $\begingroup$ @m0nhawk, oh wait, you are right, they are different. I probably messed something up in the arguments for the equation. And the coefficient in front is not important, so I left it out. Thank you for catching that $\endgroup$
    – Yuriy S
    Nov 17, 2018 at 23:11
  • $\begingroup$ @m0nhawk, I fixed the arguments in the formula $\endgroup$
    – Yuriy S
    Nov 17, 2018 at 23:12
  • $\begingroup$ This seems to be one of those cases where doing a straight numerical integration yields more accurate results than trying to derive a numerically unstable closed form. (Acton devotes a few paragraphs on this subject in "Numerical Methods That Work".) $\endgroup$ Jan 5, 2019 at 15:35

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ "It would be helpful if you showed all code used." - I did show all the important code. It works on its own and doesn't require anything undefined. //// Thank you for the answer, I'm not sure I understand what you did specifically, but using Integrate is not convenient for me, I need this to work as fast as possbile. $\endgroup$
    – Yuriy S
    Nov 17, 2018 at 23:09
  • $\begingroup$ (1) The question is then unclear: you either do or do not want the integral, and you certainly did not show a computation thereof. If you want speed, for given values of parameters {a,b,c}, the NItegrate should be used. If the goal is to derive a symbolic expression and later do substitutions, then `Integrate can be used (and I did just that). $\endgroup$ Nov 17, 2018 at 23:20
  • $\begingroup$ (2) Possibly the unexplained issue is from using yet another variant: indefii2 = Integrate[ TrigReduce[Cos[a*x - Pi/4]*Cos[b*x - Pi/4]*Cos[c*x - Pi/4]]/ Sqrt[x], {x, 1, Infinity}, Assumptions -> 1 <= a <= b <= c] That might be a problematic antiderivative, I'm not sure. $\endgroup$ Nov 17, 2018 at 23:21
  • $\begingroup$ Please understand, the integral is just the background I decided to show. I was interested in dealing with the loss of significance in the final exact expression with Fresnel integrals. I have since solved the problem in the way I posted here. Thank you again for the answer, my final code now works fine $\endgroup$
    – Yuriy S
    Nov 17, 2018 at 23:26
  • $\begingroup$ Here is why the code matters. I also (finally) obtained an expression that suffered from bad loss of precision. But I have reason to suspect the expression itself was erroneous, and the precision loss was cancellation error insofar as it would give zero. That is to say, possibly you were running up against a buggy integral. But I cannot be certain, because I don't know what exactly you had or how to produce it myself. $\endgroup$ Nov 18, 2018 at 14:23
0
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ "Mathematica using different argument conventions than the commonly accepted ones is really inconvenient." - at the very least, Mathematica's convention is consistent with the one used by the DLMF and Abramowitz and Stegun; perhaps you are merely accustomed to the older (unnormalized) definition for the Fresnel integrals. $\endgroup$ Nov 13, 2019 at 15:23
  • $\begingroup$ @J.M., fair enough. $\endgroup$
    – Yuriy S
    Nov 14, 2019 at 5:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.