No but when I do, I get a result:

out=Integrate[BesselJ[2n+1/2,L]BesselJ[2*m+1/2,L]/Sqrt[L^2+K],{L,0,\[Infinity]},Assumptions->n\[Element]Integers&&n>=0&&m\[Element]Integers&&m>=0&&K\[Element]Reals&&K>0]
Series[out,{K,\[Infinity],4}]

outputs a very large expression. Note that you can not substitute values directly into the output, first you have to apply Normal to it. See the detailed documentation of Series.

Of course you can do series expansions around infinity:

Series[function,{x,Infinity,order}]

But getting any kind of precision might be tricky, you might need a lot of terms!

That last plot just uses machine precision, to correctly plot it use:

Plot[f[7, 0, K], {K, 0, 1000}, WorkingPrecision -> 32]

No need for N[ ], as K is numerical value anyhow...

Machine Precision is not enough precision for large K (probably due to summation):

p1 = Plot[f[7, 0, K], {K, 0, 1000}, WorkingPrecision -> 32, PlotStyle -> Red];
p2 = Plot[f[7, 0, K], {K, 0, 1000}];
Show[{p1, p2, p1}]

You see that for K > 250 the error grows considerably.

So for small K, the results are roughly the same, but the larger the K, the larger the error.

I have used N[] in Manipulate[] to decrease the time for graph manipulation!

Hi,

Yes, the huge time is one big problem, apart from the uncertainty how many terms to sum. But I have BAD NEWS: Your trick produces wrong results!

I have tried the calculations again, and here are some values of the integral obtained by your trick:

for m=n=1, K=10:
0.15432194794717644427531871148044089860

for m=2,n=1, K=10: 0.175921073744037741410030108109402757600

The first value is correct. But the second value is identical to what one would obtain from the fully symbolic formula after omitting the term multiplied by Sin[pi*(m-n)]. This is understandable, since as a result of your trick the infinity is reduced to a finite value, which is then multiplied by Sin[pi*(m-n)]=0, giving zero.

The correct value should be: -0.009153189573949895158506406248625

So, it seems that we are back in the starting point. Any ideas?

Leslaw

Dear Leslaw,

If the number of problematic points is not large, then you can use Piecewise[] function to define the result of integration. So, for the problematic points you define specific values, but for the rest domain you use the exact formula obtained by Integrate[].

Hi,

I don't see my previous message posted about two hours ago, but I have another observation. I think the problem with indeterminate results arises, because there is a coefficient 1+m-n (or 1-m+n) in the term

HypergeometricPFQ[{1/2, 1/2, 1}, {1/2 - m - n, 1 + m - n, 1 - m + n, 3/2 + m + n}, K]

The coefficient contributes to the Pochhammer symbol in the denominator of the hypergeometric series expansion. When m != n there is a possibility that the Pochhammer symbol becomes zero. Hence, the whole term becomes infinite. Later it is multiplied by Sin[pi*(m-n)], which is zero, but of course infinity * 0 is indeterminate.

I therefore think that the trick with replacing integers by reals is a bit a cheating, with a likelyhood of a numerical instability. It would be much better to obtain (if this is possible) an alternative expansion formula that avoids the problem. It is strange that MATHEMATICA cannot detect such a problem automatically, since it obtains an information that m>=0 and n>=0. Isn't this some sort of a bug?

My correct formula for the integral is:

(1/2)*Sum[(-1)^(n+m+mu) * Gamma[mu/2+1/2]^2 * (Sqrt[z])^mu/(Gamma[mu/2+n-m+1] * Gamma[mu/2-n+ m+1] * Gamma[mu/2-n-m+1/2] * Gamma[mu/2+n+m+3/2]), {mu,0,?}]

An attempt to evaluate this formula symbolically produces a yet another representation using hypergeometric series, which again turns useless because of indeterminate results. But perhaps there is some alternative method of representing the formula in such a way so that its evaluation would be fast (without the need to actually sum many terms, which is time-consuming for large K) and simultaneously accurate?

Leslaw

So, do I understand properly that there is no way to obtain the integral values with 32 digits precision (or more)?

Sure, we can by using N[2,32]:

In[1]:= Integrate[
 BesselJ[N[2, 32]*n + 1/2, L]*
  BesselJ[N[2, 32]*m + 1/2, L]/Sqrt[L^2 + K], {L, 0, \[Infinity]}, 
 Assumptions -> 
  n \[Element] Integers && n >= 0 && m \[Element] Integers && m >= 0 &&
    K \[Element] Reals && K > 0]


Out[1]= (Gamma[
     0.50000000000000000000000000000000 + 
      1.0000000000000000000000000000000 m + 
      1.0000000000000000000000000000000 n] HypergeometricPFQ[{1, 1/2, 
      1/2}, {0.50000000000000000000000000000000 - 
       1.0000000000000000000000000000000 m - 
       1.0000000000000000000000000000000 n, 
      1.50000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n, 
      1.00000000000000000000000000000000 - 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n, 
      1.00000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m - 
       1.0000000000000000000000000000000 n}, K])/(2 Gamma[
     1.00000000000000000000000000000000 + 
      1.0000000000000000000000000000000 m - 
      1.0000000000000000000000000000000 n] Gamma[
     1.00000000000000000000000000000000 - 
      1.0000000000000000000000000000000 m + 
      1.0000000000000000000000000000000 n] Gamma[
     1.50000000000000000000000000000000 + 
      1.0000000000000000000000000000000 m + 
      1.0000000000000000000000000000000 n]) + (K^(
    0.50000000000000000000000000000000 + 
     1.0000000000000000000000000000000 m + 
     1.0000000000000000000000000000000 n)
     Gamma[-0.50000000000000000000000000000000 - 
      1.0000000000000000000000000000000 m - 
      1.0000000000000000000000000000000 n] Gamma[
     1.00000000000000000000000000000000 + 
      1.0000000000000000000000000000000 m + 
      1.0000000000000000000000000000000 n]^2 Gamma[
     1.50000000000000000000000000000000 + 
      1.0000000000000000000000000000000 m + 
      1.0000000000000000000000000000000 n] \
HypergeometricPFQ[{1.50000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n, 
      1.00000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n, 
      1.00000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n}, \
{1.50000000000000000000000000000000 + 
       1.0000000000000000000000000000000 m + 
       1.0000000000000000000000000000000 n, 
      2.0000000000000000000000000000000 + 
       2.0000000000000000000000000000000 m + 
       2.0000000000000000000000000000000 n, 
      1.5000000000000000000000000000000 + 0.*10^-32 m + 
       2.0000000000000000000000000000000 n, 
      1.5000000000000000000000000000000 + 
       2.0000000000000000000000000000000 m + 0.*10^-32 n}, 
     K])/(2 \[Pi] Gamma[
     1.5000000000000000000000000000000 + 
      2.0000000000000000000000000000000 m + 0.*10^-32 n] Gamma[
     1.5000000000000000000000000000000 + 0.*10^-32 m + 
      2.0000000000000000000000000000000 n] Gamma[
     2.0000000000000000000000000000000 + 
      2.0000000000000000000000000000000 m + 
      2.0000000000000000000000000000000 n])

You use the composition Evaluate@Integrate. Will the code work the same if I first generate the formula using Integrate (and real coefficients), and next call Evaluate?

The meaning of Evaluate[] is to avoid repeated Integrate[] as you firstly Integrate[] and after Integration[] you assign the result to the function.

I still want to understand what is the reason of the error that occurs when evaluating the fully symbolic formula. Is there no method of finding out where the error actually arises? Some sort of tracing analysis?

I don't know exactly the properties of the function you investigate, but I think that it has discontinuity points which must be investigated by applying the Limit[]. When we use the approximate computations, we have the values near the discontinuity points. It's only my supposition! You must be aware of this!

Hello,

First my response to Sander: we have n>=0 and m>=0 here. This was assumed.

Now to Valeriu: well, thanks again, this helps a bit, but not entirely. I am confused and have more questions:

1) I can see that your instruction:

f[m, n, K_] := Evaluate@Integrate[ BesselJ[2.n + 1/2, L]BesselJ[2.*m + 1/2, L]/Sqrt[L^2 + K], {L, 0, [Infinity]}, Assumptions -> n [Element] Integers && n >= 0 && m [Element] Integers && m >= 0 && K [Element] Reals && K > 0]

allows me to generate correct values, but only with an accuracy of several digits, whereas I need to obtain a prescribed number of accurate digits, typically 32 or more. How should I modify the instruction to obtain such results?

2) I completely don't understand what happens here. Is the actual formula returned by Integrate[] correct or not? Why there are factors 0*infinity detected when the formula is used in symbolically derived form, but not detected in your approach?

3) Actually I'd be more happy to first have a fully correct formula that does not cause errors, and next just evaluate its value. Is there no way to do this?

4) Isn't it perhaps that Evaluate@Integrate[] uses NIntegrate rather than symbolic formula evaluation?

5) Isn't the multiplication of integers n, m by real coefficients a sort of violation of the initial assumptions? Why this helps?

Leslaw

I changed a little your formulas (2 by 2.) and obtained the following result:

In[37]:= Integrate[
 BesselJ[2.*n + 1/2, L]*BesselJ[2.*m + 1/2, L]/Sqrt[L^2 + K], {L, 
  0, \[Infinity]}, 
 Assumptions -> 
  n \[Element] Integers && n >= 0 && m \[Element] Integers && m >= 0 &&
    K \[Element] Reals && K > 0]

Out[37]= 1/2 ((
   Gamma[1/2 + 1. m + 1. n] HypergeometricPFQ[{1/2, 1, 1/
      2}, {1/2 - 1. m - 1. n, 3/2 + 1. m + 1. n, 1 - 1. m + 1. n, 
      1 + 1. m - 1. n}, K])/(
   Gamma[1 + 1. m - 1. n] Gamma[1 - 1. m + 1. n] Gamma[
     3/2 + 1. m + 1. n]) + (K^(1/2 + 1. m + 1. n)
       Gamma[-(1/2) - 1. m - 1. n] Gamma[1 + 1. m + 1. n]^2 Gamma[
       3/2 + 1. m + 1. n] HypergeometricPFQ[{1 + 1. m + 1. n, 
        3/2 + 1. m + 1. n, 1 + 1. m + 1. n}, {3/2 + 1. m + 1. n, 
        2 + 2. m + 2. n, 1.5 + 2. n, 1.5 + 2. m}, K])/(\[Pi] Gamma[
       1.5 + 2. m] Gamma[1.5 + 2. n] Gamma[2 + 2. m + 2. n]))

Now, I can animate the result:

f[m_, n_, K_] := 
 Evaluate@Integrate[
   BesselJ[2.*n + 1/2, L]*BesselJ[2.*m + 1/2, L]/Sqrt[L^2 + K], {L, 
    0, \[Infinity]}, 
   Assumptions -> 
    n \[Element] Integers && n >= 0 && m \[Element] Integers && 
     m >= 0 && K \[Element] Reals && K > 0]

t = Flatten@
   Table[Plot[f[m, n, K], {K, 0, 10}], {m, 0, 20, 1}, {n, 0, 20, 1}];

Hope, this will help you!

The general solution (so for every m and n) might not cover all cases for m and n.

Example:

Say I integrate:

Integrate[x^n, x]

Giving:

x^(1 + n)/(1 + n)

However this does not cover the n = -1 case:

Integrate[x^-1, x]

Gives:

Log[x]