2

I have the following type of integral

Integrate[ r BesselJ[n, a r] BesselJ[n, b r], {r, 0, Infinity}

(where a and b are real and $ n $ are integers) which Mathematica tells me it diverges.

We know that the result is actually

$$ \int_0^{\infty} J_\mu (a r) J_\mu (b r) \quad r \quad dr = \frac{\delta(a-b)}{a} $$

The problem is the Integrate function is unable to recognize a delta function. I have tried a few options that were suggested as using FourierTransform (which does not work because the expression seems to complicated to Fourier transform it) and also the TagSetDelayed option (Teaching Mathematica more about DiracDelta and KroneckerDelta ).

In any case I did not manage to solve that. Is there any form of doing it?

xzczd
  • 65,995
  • 9
  • 163
  • 468
Marc Borrell
  • 21
  • 1
  • 1
    The integral doesn't depend on n? – Ulrich Neumann Nov 04 '20 at 10:30
  • I corrected the expression, there was a mistake. – Marc Borrell Nov 04 '20 at 10:36
  • @ MarcBorrell Still the righthandside of your equation doesn't depend on \[Mu]. Where did you find this formula? – Ulrich Neumann Nov 04 '20 at 10:41
  • https://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/21/02/02/0006/ It is from here. The result does not depend on mu. – Marc Borrell Nov 04 '20 at 10:49
  • Here you find a "confirmation" of your formula Integrate[ t BesselJ[\[Nu], a t] BesselJ[\[Nu], b t], {t, 0, Infinity}] == (1/a) DiracDelta[a - b] /; Element[a, Reals] && Element[b, Reals] && Element[\[Nu], Reals] – Ulrich Neumann Nov 04 '20 at 11:08
  • Does that come from Mathematica? – Marc Borrell Nov 04 '20 at 11:13
  • Yes, from link you gave! – Ulrich Neumann Nov 04 '20 at 11:20
  • I see. But, is there way of making Mathematica "understand" that integral? I mean if you plug in the integral it says that it does not converge. – Marc Borrell Nov 04 '20 at 11:27
  • Something strange. Integrate[BesselJ[1, 2*r]*BesselJ[1, 3*r]*r, {r, 0, Infinity}] says "Integrate::idiv: Integral of r BesselJ[1,2 r] BesselJ[1,3 r] does not converge on {0,[Infinity]}." and NIntegrate[BesselJ[1, 2*r]*BesselJ[1, 3*r]*r, {r, 0, Infinity}, WorkingPrecision -> 50, AccuracyGoal -> 3, PrecisionGoal -> 3, Method -> "ExtrapolatingOscillatory"] performs a warning "NumericalMathNSequenceLimit::seqlim: The general form of the sequence could not be determined, and the result may be incorrect." and 0.096 and the numeric integral over $(0,2000)$ confirms it. – user64494 Nov 04 '20 at 11:34
  • Gradshtein&Ruezhik includes 6.521.1 with the identical integrand. The result is rather like the Kronecker delta. I have no time to verify that. – user64494 Nov 04 '20 at 12:17
  • @user64494 The formula is correct, and it is not the Kronecker delta. – yarchik Nov 04 '20 at 12:23
  • @yarchik: Can you kindly ground your claim? TIA. – user64494 Nov 04 '20 at 12:53
  • Normal[Series[BesselJ[1, 2*r]^2*r, {r, Infinity, 2}]] results in (1/\[Pi] + 15/(256 \[Pi] r^2)) Cos[\[Pi]/4 + 2 r]^2 - ( 3 Cos[\[Pi]/4 + 2 r] Sin[\[Pi]/4 + 2 r])/(8 \[Pi] r) + ( 9 Sin[\[Pi]/4 + 2 r]^2)/(256 \[Pi] r^2) and this implies the diververgence (if I or Mathematica are not mistaken). – user64494 Nov 04 '20 at 14:37
  • Normal[Series[r*BesselJ[1, 2*r]*BesselJ[1, 3*r], {r, Infinity, 2}]];Integrate[%, {r, 1, Infinity}] produces an error communication "Integrate::idiv: Integral of <<1>>/(1536 Sqrt[6] [Pi] r^2) does not converge on {1,[Infinity]}.". In view of it the formula under consideration does not seem to be true. Likely a typo. – user64494 Nov 04 '20 at 20:11
  • 1
    @user64494 Why do you think the series can be integrated here? What you are doing is not mathematically rigorous. Moreover, you comments discourage other people (me including) to provide answers! – yarchik Nov 05 '20 at 06:46
  • @yarchik: See here (in Russian). – user64494 Dec 04 '20 at 12:20
  • @user64494 I am sorry, I do not see the connection. According to the Fubini theorems, one can exchange series expansion and integration when either left or right hand sides are finite. But your series expansion contains $1/r^2$, therefore MA is correct in saying that integral of series does not converge. But this makes no implication about the convergence or divergence of the original integral. At least I do not see. – yarchik Dec 04 '20 at 13:28
  • @yarchik: Up to to the cited link, the integrals of r*BesselJ[1, 2*r]*BesselJ[1, 3*r] and Normal[Series[r*BesselJ[1, 2*r]*BesselJ[1, 3*r], {r, Infinity, 2}]] converge or diverge simultaneously. – user64494 Dec 04 '20 at 13:41
  • This integral "almost" converges, insofar as it converges if your factor changes from r to r^(1-epsilon). So you can do e.g. regularizedDelta = Integrate[r^eps1*BesselJ[nu, a*r]*BesselJ[nu, b*r], {r, 0, Infinity}, Assumptions -> {a > b > 0, nu > 0, 0 < eps1 < 1}]. A plot will then confirm the delta-like result: Plot[Re[Evaluate[ regularizedDelB /. {nu -> 1, a -> 3, eps1 -> .99}]], {b, 2, 3}, PlotRange -> All]. To be sure it has the correct factor (which I would expect to be 1/(a-b) if a>b by the way), one might numerically integrate against test functions. – Daniel Lichtblau Dec 04 '20 at 23:58
  • Actually what I wrote might be off in terms of regularization: possibly a needs to be less than 0 rather than 1.. It still seems to work... – Daniel Lichtblau Dec 05 '20 at 00:12
  • 1
    @yarchik No need to break integrand into an infinite series and then have to worry about swapping integral with summation. Just split into two (or possibly three) lead terms and all the rest. The "all the rest" part can be shown to converge. The lead term part can be shown to not converge due to oscillatory behavior. Ergo, divergence for the integral in the classical sense. – Daniel Lichtblau Dec 05 '20 at 15:33

2 Answers2

2

I didn't succeed in making Mathematica "understand" your integral.

The wellknown simpler case Integrate[Exp[I \[Omega] t], {t, -Infinity, Infinity}] == 2 Pi DiracDelta[\[Omega]] isn't understood by Mathematica too:

Integrate[Exp[I ω t], {t, -Infinity, Infinity}] 
(*Integrate::idiv: Integral of E^(I t ω) does not converge on {-∞,∞}.*)

But workaround AsymptoticIntegrate evaluates to

AsymptoticIntegrate[Exp[I ω t], {t, -m, m}, {m, Infinity, 1}] 
(*(2 Sin[m ω])/ω*)

which is a limit definition of DiracDelta[\[Omega]

Avrana
  • 297
  • 1
  • 4
  • 14
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • What do you mean by "a limit definition"? – user64494 Nov 04 '20 at 11:44
  • 1
    Something like DiracDelta[x]=Limit[Sin[x/eps]/(Pi x),eps->0](see Wikipedia). – Ulrich Neumann Nov 04 '20 at 11:53
  • Unfortunately, this is not a usual limit, but a limit in the weak topology. See Wiki, Represetations of the delta function to this end. Mathematica does not deal with such limits. – user64494 Nov 04 '20 at 12:12
  • 2
    @user64494 The limit is wellknown in mathematics, physics, signal theory,... Sorry I' m not interested in refreshing this neverending story concerning DiracDelta – Ulrich Neumann Nov 04 '20 at 12:27
  • UlrichNeumann (@ does not work), can you give a reference to ground your claim? – user64494 Dec 04 '20 at 12:32
  • @user64494 the Wikipedia page on Fourier transforms states that the Fourier transform of $1$ is $\delta(x)$ (in the sense of distributions). – Roman Dec 04 '20 at 12:41
  • @Roman: Sorry, but Ulrich Neumann confuses a usual limit with a weak limit. – user64494 Dec 04 '20 at 13:24
  • @user64494 Sorry I didn't confuse anything! Perhaps my answer is based on a weak limit. But since Dirac introduced the concept, it is well accepted and succesfully applied in many applications! – Ulrich Neumann Dec 04 '20 at 13:34
  • @UlrichNeumann:-1. No CAS deals with a weak limit. – user64494 Dec 04 '20 at 13:46
  • 1
    I think this is quite a fine answer (and, apparently, I negated a downvote). It's a nice use of AsymptoticIntegrate to generalized functions. And a reminder (to me) that the cardinal sine gives a way to approximate the Dirac delta as an actual function. (I usually think about that in terms of Gaussians. Maybe because my head is bell-shaped...) – Daniel Lichtblau Dec 04 '20 at 14:59
  • 1
    @DanielLichtblau Thanks for your support. It's difficult for me to get used of "downvoting" as a feedback. By the way a very useful limit for Dirac Delta is gaussian ;-) : dirac[x]~ Limit[ Exp[-(x/eps)^2]/(eps Sqrt[Pi]),eps->0] – Ulrich Neumann Dec 04 '20 at 15:22
  • @DanielLichtblau:I prefer formulas and references over emotional words. – user64494 Dec 04 '20 at 20:09
  • @user64494 No mystery intended. What I alluded to is a well known regularization used to compute/approximate generalized functions. Multiply integrand by Exp[-a*x^2], integrate under assumption that a>0, take limit as a->0. Check `In[5]:= regularizedFT = Integrate[Exp[-a*x^2]Exp[Iwx]Cos[x], {x, -Infinity, Infinity}, Assumptions -> {a > 0, Element[w, Reals]}]

    Out[5]= (E^(-((1 + w)^2/(4 a))) (1 + E^(w/a)) Sqrt[[Pi]])/(2 Sqrt[a])and then do for examplePlot[Abs[regularizedFT /. a -> .01], {w, -3, 3}, PlotRange -> All]` to see it approximate the FT of the cosine function.

     – Daniel Lichtblau Dec 04 '20 at 23:32
2

A workaround and under certain assumptions we have:

func= r*BesselJ[n, a r]*BesselJ[n, b r];

InverseMellinTransform[Integrate[MellinTransform[func, a, s], {r, 0, Infinity}, 
Assumptions -> {s > 1, b > 0, n [Element] Integers, n >= 0, 2 + n > s}], s, a]


(DiracDelta[a - b]/b)

Maple 2020.2 Can deal:

enter image description here

Mariusz Iwaniuk
  • 13,841
  • 1
  • 25
  • 41