I have to evaluate the object

$$ b(q,m,n)=\int_{-\infty}^{\infty} \mathrm{d}k |c_n(w(k,q))|^2J_m(q)\delta(Y_n(w(k,q))-m) $$

Where $q$ is a continuous variable, $m$, $n$ are integers and

$$ w(k,q)=k^2+q^2 $$

The exact definitions of $c_n$ and $Y_n$ are irrelevant, the important point being that they are both known only numerically (no closed form expression).

Nevertheless in my problem they are

$$ c_n(w)= \int_{0}^{2\pi}\mathrm{d}z\, S(b_{2n+2}(w), \frac{1}{w^2}, z) C'(a_0(w), \frac{1}{w^2}, z) $$




$C(a,q,z)$, $S(a,q,z)$ - Mathieu even and odd functions (aka elliptic cosine and sine functions), $a_0(z)$, $b_{2n+2}(z)$ - Mathieu characteristic values for even and odd functions.

Here is how I am trying to solve this:

First I define the variables $Y_n, w, c_n$ above

w[q_, k_] := q^2+k^2
Y[n_, q_,k_] :=  MathieuCharacteristicB[2 n + 2, w[q,k]] - 
MathieuCharacteristicA[0, w[q,k]]
c[n_,q_,k_] := NIntegrate[MathieuS[MathieuCharacteristicB[2n + 2, w[q,k]], 1/w[q,k]^2, x] MathieuCPrime[MathieuCharacteristicA[0, w[q,k]], 1/w[q,k]^2, x], {x, 0, 2 Pi}]

Until this point everything seems to work fine. I can evaluate the above objects. But then I try to evaluate $b$ using

b[q_,m_,n_ ] := Integrate[Abs[c[n,q,k]]^2 BesselJ[m, q] DiracDelta[Y[n,q,k] - m], {k, -Infinity,Infinity}]

But this returns the error messages

NIntegrate::inumri: The integrand Indeterminate has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0,6.8}}. 

NIntegrate::nlim: x = Integrate`ImproperDump`zero is not a valid limit of integration.

Any help in overcoming these issues would be greatly appreciated.

  • $\begingroup$ Please format your code in Mathematica format not $\TeX$ $\endgroup$ Aug 19 '16 at 3:29
  • $\begingroup$ This is not the code. This is the statement of the problem I'm trying to solve. $\endgroup$ Aug 19 '16 at 8:42
  • $\begingroup$ I've added the code now. $\endgroup$ Aug 19 '16 at 10:31
  • $\begingroup$ Your code is garbled. For example DiracDelta[E_n - m*r}] makes no sense. Please repair your code and format it for this forum. I suggest you concentrate on the integral containing the delta function. $\endgroup$
    – John Doty
    Aug 19 '16 at 16:05
  • 2
    $\begingroup$ I would appreciate specific advice, pertinent to the problem posted above. $\endgroup$ Aug 20 '16 at 2:05

An extended comment:

Since you're working numerically, a quick and dirty approximation might work. Simply replace the delta function with an approximation.

For example, you might use (Edit: I replaced DiracDelta with an undefined function diracDelta)

diracDelta[x_]:>Exp[-x^2/a^2]/(a Sqrt[Pi]) /. a-> 10^(-15)

You can replace 10^(-15) with whatever value you need for your precision.

Full code:

b[q_,m_,n_ ] = Integrate[Abs[c[n,q,k]]^2 BesselJ[l, q] diracDelta[Y[n,q,k] - m]/.diracDelta[x_]:>Exp[-x^2/a^2]/(a Sqrt[Pi]) /. a-> 10^(-15), {k, -Infinity,Infinity}]

A few more tips:

  • Another approximation (such as a triangle function) might be better suited for your purposes, if the delta function approximation need not be differentiable.
  • You might want to use NIntegraterather than Integrate, if the arguments to Integrate are numeric.
  • You might want to use := (SetDelayed) rather than = (Set). Your first block of code fails to run unless I make that swap.

There seem to be some problems with your code, making it difficult to run and test this solution. Could you post a few sets (or ranges) of example parameters k,q,n? Is k a continuous variable? Note that I got an error MathieuCharacteristicB::zord: There is no zero-order MathieuCharacteristicB.

  • $\begingroup$ Also, welcome to the world of Mathematica! $\endgroup$
    – jjc385
    Aug 20 '16 at 8:52
  • $\begingroup$ Thanks! That is an interesting suggestion, that I will now look into. There were some remaining typos in the code I had presented, which I have now corrected. Now the part before the delta function works fine. k and q are continuous variables, n and m are integers. Thanks again! $\endgroup$ Aug 20 '16 at 9:02
  • $\begingroup$ I tried your suggestion, NIntegrating over k from -100 to 100. I got The integrand (1000000000000000E^(-1000000000000000000000000000000 (-1+Y[1,1,k])^2)\Abs[c[1,1,k]]^2\BesselJ[1,1])/Sqrt[\[Pi]] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-100,100}}. $\endgroup$ Aug 20 '16 at 9:15

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.