Reply by Raymond Toy August 17, 20092009-08-17
>>>>> "Rune" == Rune Allnor <allnor@tele.ntnu.no> writes:
Rune> On 27 Jul, 16:32, Clay <c...@claysturner.com> wrote: >> Of course I would try looking in Watson[1] as he wrote the definitive >> work on Bessel functions. I would look for you, but my copy is at the >> farm. >> [1] Watson,G.N., "A Treatise on the Theory of Bessel Functions", 1922 >> Cambridge University Press. Rune> A genuine CUP edition, not a Dover reprint? Rune> Impressive. It was reprinted (not Dover) recently. I have a copy. Ray
Reply by mbtrawicki July 28, 20092009-07-28
Hello Scott,

I really appreciate your response about the closed-form solution. Here is
my results from Maple using on that integral (without the ln(x) function):


> f:=int(x^p*exp(-x^2)*BesselJ(0,sqrt(-1)*2*x),x=0..infinity);
GAMMA(p/2 + 1/2) (p + 1) LaguerreL(- p/2 + 1/2, 1) f := 1/2 -------------------------------------------------- p - 1 GAMMA(p/2 + 1/2) LaguerreL(- p/2 + 1/2, 1, 1) - --------------------------------------------- p - 1
> combine(convert(f,hypergeom));
(1/4 hypergeom([p/2 - 1/2], [1], 1) p + 1/4 hypergeom([p/2 - 1/2], [1], 1) + 1/4 hypergeom([p/2 - 1/2], [2], 1) p - 3/4 hypergeom([p/2 - 1/2], [2], 1)) GAMMA(p/2 - 1/2) I am sure that my result simplifies even further into your result too. Are there maybe any mathematical papers that possibly deal with integrals involving ln(x)? I just really believe that there has to be a closed-form solution for the integral integral(x^p*ln(x)*exp(-x^2)*besseli(0,2*x),x=0..infinity) even if Maple and/or Mathematica do not provide one. Thank you again, Marek
Reply by mbtrawicki July 28, 20092009-07-28
Hello Clay,

Thank you for the hint with performing the numerical integration on the
integral

integral(x^p*ln(x)*exp(-x^2)*besseli(0,2*x),x=0..infinity)

using Gaussian Quadrature. I know that Matlab has some build in functions
("quad" and variants), but I do not think that they have Gaussian
Quadrature. I am using Maple commands in Matlab to perform the integration,
but I probably have to code it up myself. Is there any way to simplify
maybe parts of the integrand (not necessarily using the Bessel function
approximations for large inputs)? I use the approximation

Io(y) ~ (1/sqrt(2*pi*y))*exp(y) for large y

but still have issues finding the integral in a table. I wonder whether
somehow substituting in another simplification for any combination of
functions in the integral would then allow me to use an entry in a table.
The problems are naturally with the ln(x) and/or Bessel function. I have
even tried integration by parts but ran into trouble with the ln(x) in u*v
and corresponding integral for dv in the relationship

integral(u*v) = u*v - integral(v*du),

where the two integrals and u*v are all evaluated from 0 to infinity. I do
not see really any way to find a closed-form solution but am still working
on it. 

Thank you again for the advice and supporting documents,

Marek
Reply by Clay July 27, 20092009-07-27
On Jul 27, 4:40&#2013266080;pm, Clay <c...@claysturner.com> wrote:
> On Jul 27, 3:34&#2013266080;pm, "mbtrawicki" <mbtrawi...@yahoo.com> wrote: > > > Hello Clay, > > > I verified your results and have very similar numbers. I am using Maple to > > compute those integrals. Is there no closed-form solution? I still have not > > found one in any of my references for the integral. Do you know of a faster > > implementation of Maple in Matlab? I am using Maple functions in Matlab > > since the integration is a major component of a bigger program, but my > > Matlab is running so slowly this way. I am not sure whether there is any > > way around it. I welcome all suggestions. > > > Thank you again, > > > Marek > > Hello Marek, > > An obvious closed form solution does not pop into my head. I'm using > MathCad but not using Maple. > > A chunk of 'c' code to generate the Gaussian Quad coefs is not too > hard to come up with. I.e., the integrand can be sampled at about 100 > to 200 points, and from that you would get a very precise result and > should execute in fractions of a second. If &#2013266080;-1 < p < 20, then you > will only need to integrate from 0 to 8. The product I0(2x)*e^-(x*x) > acts alot like &#2013266080;e^-(x/1.7)^3, so unless p is big, the exponential term > crushes x^p to 0. > > Look up how to calculate the Gaussian Quadrature coefs and code that > up. The integral becomes trivial to do from there. > > IHTH, > Clay
See here for how to quickly integrate this function. There is a graph of the integral for p = -1 up to 5, just so you can see how it looks. This MathCad sheet ran in under a second with 1000 values for p. It uses Gaussian Quadrature. http://www.claysturner.com/dsp/x.pdf IHTH, Clay
Reply by Scott Hemphill July 27, 20092009-07-27
Clay <clay@claysturner.com> writes:

> On Jul 27, 3:23&nbsp;pm, Clay <c...@claysturner.com> wrote: >> On Jul 27, 2:48&nbsp;pm, Clay <c...@claysturner.com> wrote: >> >> >> >> >> >> > On Jul 27, 10:00&nbsp;am, "mbtrawicki" <mbtrawi...@yahoo.com> wrote: >> >> > > Hello, >> >> > > I am trying to find a closed-form solution (or reasonable approximation) >> > > to the following integral: >> >> > > integral(x^p*ln(x)*exp(-x^2)*besseli(0,2*x),x=0..infinity), >> >> > > where p is a constant and besseli(0,2*x) is the zeroth order Bessel >> > > function of the first kind. >> >> > > I have already looked through the traditional resources (Table of >> > > Integrals, Series, and Products by Gradshteyn and Ryzhik and Handbook of >> > > Mathematical Functions by Abramowitz and Stegun) for a closed-form solution >> > > and approximations to the Bessel function and have attempted numerical >> > > integration but have not had any luck determining the final result. Do any >> > > of you have suggestions? >> >> > > I would greatly appreciate any feedback and assistance. >> >> > > Thank you again, >> >> > > Marek >> >> > Here are some rough estimates for integer values of p from 0 to 10: >> >> > integral_0_infty of &nbsp;(x^p)*(e^-x^2)*(ln(x))*J0(2x) >> >> > 1st value is for p==0, 2nd for p==1, and so on >> >> > -0.864 >> > -0.173 >> > -0.079 >> > -0.066 >> > -0.086 >> > -0.142 >> > -0.259 >> > -0.484 >> > -0.949 >> > -1.831 >> > -3.51 >> >> > I hope this helps. After looking at the integrand, it is not a bad one >> > to partition (split the integral into several smaller pieces). What >> > values do you envision "p" to have? This looks like it can be done by >> > Gaussian Quadrature easily enough. I did these estimates by Monte- >> > Carlo. >> >> > IHTH, >> >> > Clay- Hide quoted text - >> >> > - Show quoted text - >> >> With the integrang changed to reflect the modified Bessel function, I >> find these rough values for p==0 up to 10 >> >> -0.764 >> 0.155 >> 0.622 >> 1.41 >> 3.007 >> 6.848 >> 15.919 >> 39.671 >> 100.676 >> 275.312 >> 760.466 >> >> Does this help? >> >> Clay- Hide quoted text - >> >> - Show quoted text - > > Disreragard the earlier values. > > > Now for the integrad I have > > (x^p)(e^-(x*x))*I0(2x) > > This yields for the p==0 up to p==10 > > 1.558 > 1.344 > 1.739 > 2.682 > 4.86 > 9.573 > 20.183 > 46.424 > 112.976 > 280.073 > 756.176 > > Clay
And Mathematica's numbers for these are: 1.5538993500654319234 1.3591409142295226177 1.7423093452556455141 2.7182818284590452354 4.8384531982505785615 9.5139863996066583238 20.272069964427690401 46.210791083803769001 111.66415726192771680 284.06045107397022710 756.64455829311024375 Scott -- Scott Hemphill hemphill@alumni.caltech.edu "This isn't flying. This is falling, with style." -- Buzz Lightyear
Reply by Scott Hemphill July 27, 20092009-07-27
Clay <clay@claysturner.com> writes:

> On Jul 27, 2:48&nbsp;pm, Clay <c...@claysturner.com> wrote: >> On Jul 27, 10:00&nbsp;am, "mbtrawicki" <mbtrawi...@yahoo.com> wrote: >> >> >> >> >> >> > Hello, >> >> > I am trying to find a closed-form solution (or reasonable approximation) >> > to the following integral: >> >> > integral(x^p*ln(x)*exp(-x^2)*besseli(0,2*x),x=0..infinity), >> >> > where p is a constant and besseli(0,2*x) is the zeroth order Bessel >> > function of the first kind. >> >> > I have already looked through the traditional resources (Table of >> > Integrals, Series, and Products by Gradshteyn and Ryzhik and Handbook of >> > Mathematical Functions by Abramowitz and Stegun) for a closed-form solution >> > and approximations to the Bessel function and have attempted numerical >> > integration but have not had any luck determining the final result. Do any >> > of you have suggestions? >> >> > I would greatly appreciate any feedback and assistance. >> >> > Thank you again, >> >> > Marek >> >> Here are some rough estimates for integer values of p from 0 to 10: >> >> integral_0_infty of &nbsp;(x^p)*(e^-x^2)*(ln(x))*J0(2x) >> >> 1st value is for p==0, 2nd for p==1, and so on >> >> -0.864 >> -0.173 >> -0.079 >> -0.066 >> -0.086 >> -0.142 >> -0.259 >> -0.484 >> -0.949 >> -1.831 >> -3.51 >> >> I hope this helps. After looking at the integrand, it is not a bad one >> to partition (split the integral into several smaller pieces). What >> values do you envision "p" to have? This looks like it can be done by >> Gaussian Quadrature easily enough. I did these estimates by Monte- >> Carlo. >> >> IHTH, >> >> Clay- Hide quoted text - >> >> - Show quoted text - > > With the integrang changed to reflect the modified Bessel function, I > find these rough values for p==0 up to 10 > > -0.764 > 0.155 > 0.622 > 1.41 > 3.007 > 6.848 > 15.919 > 39.671 > 100.676 > 275.312 > 760.466
Mathematica's numbers are: -0.76932709667515764264 0.14908684058079851859 0.63074368868670269879 1.4073145953911196549 3.0499225104518270595 6.8393124552132027185 16.055428433081164515 39.484039092403305545 101.50192031004547648 271.98733252729141392 757.55019687159582342 Scott -- Scott Hemphill hemphill@alumni.caltech.edu "This isn't flying. This is falling, with style." -- Buzz Lightyear
Reply by Scott Hemphill July 27, 20092009-07-27
Scott Hemphill <hemphill@hemphills.net> writes:

> "mbtrawicki" <mbtrawicki@yahoo.com> writes: > >> Hello Clay, >> >> I checked my references a few times and do see that >> >> I0(y) ~ (1/sqrt(2*pi*y))*exp(y) for large y. >> >> In my integral, the integrand is actually written as >> >> INTEGRAL(x^p*exp(-x^2)*J0(sqrt(-1)*2*x),x=0..infinity), >> >> which can be written as >> >> INTEGRAL(x^p*exp(-x^2)*I0(2*x),x=0..infinity), >> >> where >> >> In(z)=(sqrt(-1)^-n)*Jn(sqrt(-1)*z) >> >> or >> >> I0(z)=J0(sqrt(-1)*z). >> >> I take it that there is no closed-form solution to this integral? > > It depends what you accept for a closed-form. An old version of > Mathematica gives > > 1 + p 1 + p > Gamma[-----] Hypergeometric1F1[-----, 1, 1] > 2 2 > ------------------------------------------- > 2 > > as long as the real part of p is greater than -1.
However, this integrand doesn't contain the factor "ln(x)", which is what you were really looking for. Scott -- Scott Hemphill hemphill@alumni.caltech.edu "This isn't flying. This is falling, with style." -- Buzz Lightyear
Reply by Clay July 27, 20092009-07-27
On Jul 27, 3:34&#2013266080;pm, "mbtrawicki" <mbtrawi...@yahoo.com> wrote:
> Hello Clay, > > I verified your results and have very similar numbers. I am using Maple to > compute those integrals. Is there no closed-form solution? I still have not > found one in any of my references for the integral. Do you know of a faster > implementation of Maple in Matlab? I am using Maple functions in Matlab > since the integration is a major component of a bigger program, but my > Matlab is running so slowly this way. I am not sure whether there is any > way around it. I welcome all suggestions. > > Thank you again, > > Marek
Hello Marek, An obvious closed form solution does not pop into my head. I'm using MathCad but not using Maple. A chunk of 'c' code to generate the Gaussian Quad coefs is not too hard to come up with. I.e., the integrand can be sampled at about 100 to 200 points, and from that you would get a very precise result and should execute in fractions of a second. If -1 < p < 20, then you will only need to integrate from 0 to 8. The product I0(2x)*e^-(x*x) acts alot like e^-(x/1.7)^3, so unless p is big, the exponential term crushes x^p to 0. Look up how to calculate the Gaussian Quadrature coefs and code that up. The integral becomes trivial to do from there. IHTH, Clay
Reply by Scott Hemphill July 27, 20092009-07-27
"mbtrawicki" <mbtrawicki@yahoo.com> writes:

> Hello Clay, > > I checked my references a few times and do see that > > I0(y) ~ (1/sqrt(2*pi*y))*exp(y) for large y. > > In my integral, the integrand is actually written as > > INTEGRAL(x^p*exp(-x^2)*J0(sqrt(-1)*2*x),x=0..infinity), > > which can be written as > > INTEGRAL(x^p*exp(-x^2)*I0(2*x),x=0..infinity), > > where > > In(z)=(sqrt(-1)^-n)*Jn(sqrt(-1)*z) > > or > > I0(z)=J0(sqrt(-1)*z). > > I take it that there is no closed-form solution to this integral?
It depends what you accept for a closed-form. An old version of Mathematica gives 1 + p 1 + p Gamma[-----] Hypergeometric1F1[-----, 1, 1] 2 2 ------------------------------------------- 2 as long as the real part of p is greater than -1. Scott -- Scott Hemphill hemphill@alumni.caltech.edu "This isn't flying. This is falling, with style." -- Buzz Lightyear
Reply by mbtrawicki July 27, 20092009-07-27
Hello Clay,

I verified your results and have very similar numbers. I am using Maple to
compute those integrals. Is there no closed-form solution? I still have not
found one in any of my references for the integral. Do you know of a faster
implementation of Maple in Matlab? I am using Maple functions in Matlab
since the integration is a major component of a bigger program, but my
Matlab is running so slowly this way. I am not sure whether there is any
way around it. I welcome all suggestions.

Thank you again,

Marek