# NIntegrate: Catastrophic loss of precision

Below is a simplified code that I am trying to evaluate. I am basically getting the rms value of a Wakefield function. When I finally try to evaluate the rms, I get an error. Any help would be greatly appreciated.

sigmaS = 50*10^-6;

gaussian2[s_?NumericQ] := 1/(Sqrt[2 Pi] sigmaS) Exp[-(s^2/(2 sigmaS^2))]

fo2[s_?NumericQ] :=
NIntegrate[gaussian2'[sprime]/(s - sprime)^(1/3), {sprime, -Infinity, s}]

rms2 = Sqrt[
NIntegrate[
gaussian2[x] (fo2[x])^2, {x, -Infinity, Infinity}] - (NIntegrate[
gaussian2[x] fo2[x], {x, -Infinity, Infinity}])^2]

• Please provide a complete and working (even if that means working to an error) example. As it is, your code is non-functional. – ciao Jun 19 '14 at 7:34
• Whoops, that should do it. I was missing a constant! – user1886681 Jun 19 '14 at 8:13
• How long does it take to spit out the error? – dr.blochwave Jun 19 '14 at 11:57

The integral for fo2 can be computed exactly, which makes the numerical computation of rms straightforward:

Clear[s, gaussian2, fo2];

sigmaS = 50*10^-6;

gaussian2[s_] := 1/(Sqrt[2 Pi] sigmaS) Exp[-(s^2/(2 sigmaS^2))];

fo2[s_] = Integrate[gaussian2'[sprime]/(s - sprime)^(1/3), {sprime, -Infinity, s},
Assumptions -> s < 0];

rms2 = Sqrt[
NIntegrate[gaussian2[x] (fo2[x])^2, {x, -Infinity, Infinity}] -
(NIntegrate[gaussian2[x] fo2[x], {x, -Infinity, Infinity}])^2
]

(*
77003.5
*)


Note: The integral for fo2 evaluates to a ConditionalExpression with Re[s] > 0, if no assumptions are given or if the assumption that s is real is given. However we get two results with either assumption that s > 0 or s < 0, both of which are equivalent. Hence the assumption s < 0 above.

intpos = Integrate[
gaussian2'[sprime]/(s - sprime)^(1/3), {sprime, -Infinity, s},
Assumptions -> s > 0]
(*
180000/7 10^(1/3) ((7 Gamma[2/3] Hypergeometric1F1[2/3, 1/2, -200000000 s^2])/
Gamma[-(5/3)] - (
7500 2^(1/6)
s Gamma[13/3] Hypergeometric1F1[7/6, 3/2, -200000000 s^2])/
Sqrt[π])
*)

intneg = Integrate[
gaussian2'[sprime]/(s - sprime)^(1/3), {sprime, -Infinity, s},
Assumptions -> s < 0]
(*
(1/Sqrt[π])200000 5^(1/3) (Gamma[5/6] Hypergeometric1F1[2/3, 1/2, -200000000 s^2] -
10000 Sqrt
s Gamma[4/3] Hypergeometric1F1[7/6, 3/2, -200000000 s^2])
*)

intpos - intneg // FullSimplify
(*
0
*)