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):
(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�pm, Clay <c...@claysturner.com> wrote:
> On Jul 27, 3:34�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
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 pm, Clay <c...@claysturner.com> wrote:
>> On Jul 27, 2:48 pm, Clay <c...@claysturner.com> wrote:
>>
>>
>>
>>
>>
>> > On Jul 27, 10:00 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 (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 pm, Clay <c...@claysturner.com> wrote:
>> On Jul 27, 10:00 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 (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�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