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$
  • $\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 '18 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 '18 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 '18 at 23:11
  • $\begingroup$ @m0nhawk, I fixed the arguments in the formula $\endgroup$ – Yuriy S Nov 17 '18 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$ – J. M. is away Jan 5 at 15:35
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$
  • $\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 '18 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$ – Daniel Lichtblau Nov 17 '18 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$ – Daniel Lichtblau Nov 17 '18 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 '18 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$ – Daniel Lichtblau Nov 18 '18 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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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