Expression involving delta function

Question

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 = Integrate`ImproperDump`zero is not a valid limit of integration.

Any help in overcoming these issues would be greatly appreciated.

Answer 1

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.