1
$\begingroup$

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:

enter image description here

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

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Jason B.
    Feb 9, 2016 at 9:06

1 Answer 1

4
$\begingroup$

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]

enter image description here

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 $\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]

enter image description here

$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Feb 9, 2016 at 10:25
  • $\begingroup$ @Farhad - I did the same thing (can't believe I left out that part from my explanation lol) $\endgroup$
    – Jason B.
    Feb 9, 2016 at 10:29
  • 1
    $\begingroup$ an alternate approach is to make your parameter dependence explicit when you define the symbols, ie define s[d_] = , etc $\endgroup$
    – george2079
    Feb 9, 2016 at 14:58

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.