numerical integration errors

I have plotted the following integrals with numerical integration nicely without any error:

Λ = 10^-6;
Δ = 10^-3;
θ = 1/2 ArcTan[Δ/δ];
h = 10^-1;
t = 10^3;
s = -h Sqrt[Δ^2 + δ^2]

a[δ_?NumericQ] := NIntegrate[(ω E^(-ω/Λ) (1-Cos[(ω-s)t]))/(ω-s)^2,
{ω, 0, ∞}, MaxRecursion -> 200]

b[δ_?NumericQ] := NIntegrate[(ω E^(-ω/Λ)Sin[(ω – s) t])/(ω – s)^2,
{ω, 0, ∞}, MaxRecursion -> 200]

e[δ_?NumericQ] := NIntegrate[(E^(-ω/Λ) (1-Cos[ω t]))/(ω+s),
{ω, 0, ∞}, MaxRecursion-> 200]

f[δ_?NumericQ] := NIntegrate[(E^(-ω/Λ) Sin[ω t])/(ω+s),
{ω, 0, ∞}, MaxRecursion -> 200]

m[δ_?NumericQ] := NIntegrate[(E^(-ω/Λ) (1-Cos[(ω-s)t]))/(ω-s),
{ω, 0, ∞}, MaxRecursion -> 200]

n[δ_?NumericQ] := NIntegrate[(E^(-ω/Λ) Sin[(ω-s)t])/(ω-s),
{ω, 0, ∞}, MaxRecursion -> 200]

But when I want to plot the combination of them:

Plot[Cos[θ]^2 (Abs[Sin[θ] (1 – (I t)/(π h) Λ +
Cos[2θ]]^2 ((I t)/(π h) Λ – 1/(π h) Log[1 + I t Λ]) +
Sin[2θ]^2 ((I t)/(π h) (Λ + E^(-(s/Λ))s Gamma[0, -(s/Λ)]) –
1/(π h) (a[δ] + I b[δ]))) – Cos[θ] (1/2 Sin[4θ] ((I t)/(π h)
(Λ + E^(s/Λ)s ExpIntegralEi[-(s/Λ)]) – 1/(π h) (e[δ] + I f[δ])) –
1/2 Sin[4θ] ((I t)/(π h) Λ – 1/(π h) (m[δ] + I n[δ])))])^2,
{δ, 0, 100}]

I encounter these errors:

Are these errors essential to the final result? If yes, how can I get rid of them?

=================

  

 

This is far from a minimal working example. Try to narrow down the problem, the smallest code that will give the error will make it more likely for people to help. In your case, the very first example I tried, a[5] gives a convergence error and does not return a numerical value.
– JasonB
Feb 9 at 9:06

=================

1 Answer
1

=================

You want to hunt down the error? Here is the best piece of advice: don’t plot a function until you know it works.

Okay, that’s out of the way, now let’s go through the process of finding out why your code gives an error.

First we can look at just one integral,

Λ = 10^-6;
Δ = 10^-3;
θ = 1/2 ArcTan[Δ/δ];
h = 10^-1;
t = 10^3;
s = -h Sqrt[Δ^2 + δ^2]

a[δ_?NumericQ] := NIntegrate[(ω E^(-ω/Λ) (1-Cos[(ω-s)t]))/(ω-s)^2,
{ω, 0, ∞}, MaxRecursion -> 200]

I see you want to plot for values of δ\delta between 0 and 100, so let’s test it on a value of 5:

a[5.0]

The integrand evaluated to non-numerical values. This is probably due to the fact that the parameter, δ\delta isn’t being replaced by the numerical value. Try this,

a[δ_?NumericQ] :=
With[{s = -h Sqrt[Δ^2 + δ^2], θ =
1/2 ArcTan[Δ/δ]},
NIntegrate[(ω E^(-ω/Λ) (1 –
Cos[(ω – s) t]))/(ω – s)^2, {ω,
0, ∞}, MaxRecursion -> 200]];

That seems to work,

a[5.0]
(* 7.53157*10^-12 *)

and it works for b, e, m, and n but not f. Looking at your integrands I’m not sure what makes f different from the others. But I can say that all of your integrands decrease with \omega very quickly, so there is no reason to take the integral out to infinity.

Look at what happens when we substitute {δ -> 2, ω -> 1} in all six of the integrands:

{(ω E^(-ω/Λ) (1 –
Cos[(ω – s) t]))/(ω – s)^2,
(ω E^(-ω/Λ) Sin[(ω –
s) t])/(ω – s)^2,
(E^(-ω/Λ) (1 – Cos[ω t]))/(ω +
s),
(E^(-ω/Λ) Sin[ω t])/(ω + s),
(E^(-ω/Λ) (1 –
Cos[(ω – s) t]))/(ω – s),
(E^(-ω/Λ) Sin[(ω – s) t])/(ω –
s)} /. {δ -> 2, ω -> 1}
(* {8.933430357305704*10^-434298, \
-2.020538732271493*10^-434296, 1.803453102301659*10^-434295,
3.407603228754232*10^-434295,
1.072011665210259*10^-434297, -2.424646529239256*10^-434296} *)

Since they decrease so quickly with ω\omega you can set the upper integration limit to 0.01 and you get converged results.

Now you can try to make your plot, but I would still try to evaluate the function before plotting (check me here, I had to remove a ] bracket that I assumed was a typo in your code)

Cos[θ]^2 (Abs[
Sin[θ] (1 – (I t)/(π h) Λ +
Cos[2 θ]^2 ((I t)/(π h) Λ –
1/(π h) Log[1 + I t Λ]) +
Sin[2 θ]^2 ((I t)/(π h) (Λ +
E^(-(s/Λ)) s Gamma[
0, -(s/Λ)]) –
1/(π h) (a[δ] + I b[δ]))) –
Cos[θ] (1/2 Sin[
4 θ] ((I t)/(π h) (Λ +
E^(s/Λ) s ExpIntegralEi[-(s/\
Λ)]) – 1/(π h) (e[δ] + I f[δ])) –
1/2 Sin[4 θ] ((I t)/(π h) Λ –
1/(π h) (m[δ] +
I n[δ])))])^2 /. δ -> .1
(* 0.0000249393 *)

So it will return a number (but it will not return a number for δ=0\delta=0 because it is infinite). Now you could try to plot it, but I would instead suggest creating a table and then plotting that:

list = Table[{δ,
Cos[θ]^2 (Abs[
Sin[θ] (1 – (I t)/(π h) Λ +
Cos[2 θ]^2 ((I t)/(π h) Λ –
1/(π h) Log[1 + I t Λ]) +
Sin[2 θ]^2 ((I t)/(π h) (Λ +
E^(-(s/Λ)) s Gamma[
0, -(s/Λ)]) –
1/(π h) (a[δ] + I b[δ]))) –
Cos[θ] (1/2 Sin[
4 θ] ((I t)/(π h) (Λ +
E^(s/Λ) s ExpIntegralEi[-(s/\
Λ)]) – 1/(π h) (e[δ] + I f[δ])) –
1/2 Sin[
4 θ] ((I t)/(π h) Λ –
1/(π h) (m[δ] +
I n[δ])))])^2}, {δ, 0.1,
100, .1}];~Monitor~δ

ListLinePlot[list]

  

 

Thanks man. This was really helpful. The main point is that integrands decrease with ω very quickly, so I evaluated the integrals over {ω,0 ,1} and I got no error. By the way, I got no non-numerical error for my code.
– Farhad
Feb 9 at 10:25

  

 

@Farhad – I did the same thing (can’t believe I left out that part from my explanation lol)
– JasonB
Feb 9 at 10:29

1

 

an alternate approach is to make your parameter dependence explicit when you define the symbols, ie define s[d_] = , etc
– george2079
Feb 9 at 14:58