In the course of some physics research I've been working on, a very annoying integral has appeared that I'm having difficulty evaluating numerically. Any help you could offer would be greatly appreciated. The integral to be evaluated is as follows:

$$ \int_{0}^{\infty}d\omega\int_{-\infty}^{\infty}(dk_{x})k_{x}\frac{e^{-2|k_{\bot}|}(2\operatorname{Im} R_{1})(2\operatorname{Im} R_{2})}{|1-e^{-2|k_{\bot}|}R_{1}R_{2}|^{2}}\Theta(vk_{x}-\omega) $$

Here $v$ is a constant, $\Theta$ is the Heaviside theta function, and

$$k_{\bot} = \sqrt{\frac{\omega^{2}}{c^{2}}-k_{x}^{2}} $$

where $c$ is a constant which we can take to be one. Finally, the function $R_{1}$ is given by

$$ R_{1}[\omega,k_{x}] = \cfrac{\sqrt{\cfrac{\omega^{2}}{c^{2}}-k_{x}^{2}}-\sqrt{\cfrac{\omega^{2}\left(1-\cfrac{\omega_{p}^{2}}{\omega^{2}-\omega_{0}^{2}}\right)}{c^{2}}-k_{x}^{2}}}{\sqrt{\cfrac{\omega^{2}}{c^{2}}-k_{x}^{2}}+\sqrt{\cfrac{\omega^{2}\left(1-\cfrac{\omega_{p}^{2}}{\omega^{2}-\omega_{0}^{2}}\right)}{c^{2}}-k_{x}^{2}}} $$

The function $R_{2}(\omega,k_{x}) = R_{1}(\omega-vk_{x},k_{x})$. The parameters $\omega_{0},\omega_{p}$ which appear in the $R$ functions are parameters of the theory which we are free to choose. For obvious reasons, this is not a set of functions which are easy to visualize. However, careful examination of where the integrand goes to zero can substantially reduce the range of integration which we need to consider: the simplified integral is

$$ \int_{0}^{\omega_{0}}d\omega\int_{0}^{\infty}(dk_{x})k_{x}\frac{e^{-2|k_{\bot}|}(2\operatorname{Im} R_{1})(2\operatorname{Im} R_{2})}{|1-e^{-2|k_{\bot}|}R_{1}R_{2}|^{2}}\Theta(vk_{x}-\omega) $$

The integrand is exponentially decaying with increasing $k_{x}$ for $k_{x}>\omega_{0}$, so it should be possible to truncate the $k_{x}$ integral with little error for $k_{x}\gg\omega_{0}$. My strategy for computing the integral so far has been as follows: I declare the function integrand[$\omega,k_{x}$] in Mathematica to be the aforementioned integrand; I then define

integral[p_?NumberQ]:=NIntegrate[integrand[\omega,kx]/.{v->p, other substitutions},{\omega,0,\omega_0},{kx,0,10\omega_0}]

Where "other substitutions" stands for replacements for constants in the integrand, and $10\omega_{0}$ was chosen as an arbitrary upperbound on the $k_{x}$ integral. At the suggestion of readers, I'll post the substitutions which I have been using, but please keep in mind that for all of the constants, any positive real value is allowed in principle. The parameters I've chosen are

$$ \omega_{p}\rightarrow 10^{5},\omega_{0}\rightarrow 50 $$

Mathematica can then, in principle, compute the integral given a choice of $v$. The problem with this is that there seems to be consistent underestimation of the result. I say this because the integrand is positive definite in the region of interest, so increasing the size of the integration region should never decrease the magnitude of the integral. However, I have found situations where increasing the size of the integration region causes the value of the integral to drop, which I assume is because there is some sort of inverse relationship between the accuracy of integrand sampling and the size of the region of integration.

If anyone can provide any insight into how to improve my approach, I would really appreciate the help. I have played with options like AccuracyGoal, but I'm open to suggestions if there are clever options to take besides just "increase accuracy". Thanks!

  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Jan 25 '15 at 12:27
  • 3
    $\begingroup$ I recommend that you include those "other substitutions', so that readers can run your code. $\endgroup$ – bbgodfrey Jan 25 '15 at 12:40

Your Answer

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

Browse other questions tagged or ask your own question.