1
$\begingroup$

I have an integral of the form

\begin{align} F(\omega) = \int_0^{\infty} f(s,\omega) \mathrm{d}s \end{align}

which I would like to numerically evaluate and plot for a range of $\omega \in [-30,30]$ with an arbitrary step size. Here is my work so far:

integrand[s_, omega_] := 
     Module[{sigma, xs, x1, x2, zeta, Vc, theta, t1, t2}, 
      sigma = Pi/(Pi + 2);
      xs = Exp[-Pi*s/(2*sigma)];
      x1 = -2.0*
        sigma/Pi*(Log[xs/(1 + Sqrt[1 - xs^\.1d2])] + Sqrt[1 - xs^\.1d2]);
      x2 = 1.0 - 2*sigma/Pi*(1 - xs);
      zeta = x2 + x1*I;
      Vc = 1.0/(2*sigma);
      theta = -1*ArcSin[Exp[-Pi/(2*sigma)*s]];
      t1 = 1.0/Sqrt[1 + Tan (theta)^2];
      t2 = -1.0/Sqrt[1 + 1/Tan[theta]^2];
      Re[(t1 - I*t2)/(Sqrt[zeta^2 - 1])]*Exp[I*omega*s/Vc]]


omega = Array[# &, 101, {-30, 30}]
Map[N[Integrate[integrand[s, omega], {s, 0, Infinity}]], omega]

I get no luck with this (the mathematica prompt just goes to the next line when I type shift+enter) Ideally I would like to plot $F(\omega)$ as my result. Furthermore if anyone has any particular suggestions on how to handle this particular integrand (it is singular at the lower bound and oscillatory) that would be great. I would appreciate any help, I am new to mathematica.

$\endgroup$
  • 2
    $\begingroup$ The issue here is not that your integrand is complex, it is the definition of the integrand itself. You might want to try numerical evaluation. For that you could use the pattern _?NumericQ for your integrand parameters. $\endgroup$ – Wizard Oct 9 '14 at 12:18
  • 1
    $\begingroup$ Syntax error too... Map[N[Integrate[integrand[s, omega], {s, 0, Infinity}]], omega] I assume should be Map[N[Integrate[integrand[s, #], {s, 0, Infinity}]]&, omega]. I take it you want to map over each value of the list omega in turn? $\endgroup$ – Ymareth Oct 9 '14 at 12:46
4
$\begingroup$

There are a couple syntax errors in your integrand definition. Here is a corrected version:

integrand[s_, omega_] := 
 Module[{sigma, xs, x1, x2, zeta, Vc, theta, t1, t2}, 
  sigma = Pi/(Pi + 2);
  xs = Exp[-Pi*s/(2*sigma)];
  x1 = -2.0*sigma/Pi*(Log[xs/(1 + Sqrt[1 - xs^2])] + Sqrt[1 - xs^2]);
  x2 = 1.0 - 2*sigma/Pi*(1 - xs);
  zeta = x2 + x1*I;
  Vc = 1.0/(2*sigma);
  theta = -1*ArcSin[Exp[-Pi/(2*sigma)*s]];
  t1 = 1.0/Sqrt[1 + Tan[theta]^2];
  t2 = -1.0/Sqrt[1 + 1/Tan[theta]^2];
  Re[(t1 - I*t2)/(Sqrt[zeta^2 - 1])]*Exp[I*omega*s/Vc]]

In particular, you need to use Tan[theta], and not Tan(theta); in Mathematica, functions use square brackets [], rather than round brackets ().

Also, and somewhat harder to find, was the use of 1 - xs^ 2 instead of 1 - xs^2, which causes weird (and wrong!) things to happen. In general, manually entering spaces is not needed in Mathematica, unless you are using them to separate variable names, and in some cases (like this) causes problems.

The third error is the use of N[Integrate[...]], rather than NIntegrate[...]; the two commands are different (see the Documentation), and you want to use the latter.

The fourth problem is the integrand itself, which is highly oscillatory and, on my computer, caused my laptop to run out of RAM. I'd suggest searching for previously-asked questions on how to handle such integrands.

The easiest approach might be to note that your integral appears to be the Fourier transform of a real-valued function; assuming that you don't care about the relative error of the small-magnitude values and only are interested in a small $L^2$ error, you can simply discretize the integrand (without the $\theta$-dependent exponential) and use the Fourier command to compute the Fast Fourier transform, which will be significantly faster than using NIntegrate. Since the integrand is real, smooth, non-oscillatory, and decaying, it should be trivial to handle using the FFT, with the exception being the singularity near the origin.

$\endgroup$
  • $\begingroup$ Thanks for your answer, it was very helpful. Does it make a difference that my exponential is slightly different from the standard fourier transform in that I have a 1/Vc factor in there? If I discretise, how far in s should I truncate my integrand? In using mathematica I was hoping to get some automatic clever handling integration so that I don't have to implement a numerical routine myself. $\endgroup$ – Dipole Oct 9 '14 at 14:58
  • $\begingroup$ @Jack: My knowledge of NIntegrate is honestly pretty lacking as I don't use it much; there probably is a way to handle the oscillations without causing RAM overload, but someone else will have to answer that. The factor of 1/Vc is constant, so it's not a problem, but you'll have to think carefully about what your axes are really showing. The documentation page for Fourier describes the exact definition used. $\endgroup$ – DumpsterDoofus Oct 9 '14 at 15:13

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.