# Expression involving delta function

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)$$

and

$$Y_n(w)=b_{2n+2}(w)-a_0(w)$$

where

$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 = IntegrateImproperDumpzero is not a valid limit of integration.


Any help in overcoming these issues would be greatly appreciated.

• Please format your code in Mathematica format not $\TeX$ Aug 19 '16 at 3:29
• This is not the code. This is the statement of the problem I'm trying to solve. Aug 19 '16 at 8:42
• I've added the code now. Aug 19 '16 at 10:31
• 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. Aug 19 '16 at 16:05
• I would appreciate specific advice, pertinent to the problem posted above. 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.

• Also, welcome to the world of Mathematica! Aug 20 '16 at 8:52
• 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! Aug 20 '16 at 9:02
• 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}}. Aug 20 '16 at 9:15