# Approximation to Bessel Function in Integrand

Started by July 27, 2009
```Hello Clay,

I verified your values for the integrand using Maple for p of 0 to 10. I
envision p to be defined as p>-1. Is there a fast way to implement
numerical integration in Matlab? I am using the Maple commands in Matlab,
but my implementation runs so much slower that way than simply typing in
the integral for various p values directly into Maple. Unfortunately, I
have to use Matlab and its vector capabilities. Here is my implementation:

p=2;
ftn=['x^' num2str(pp) '*ln(x)*exp(-x^2)*BesselI(0,2*x)'];
[r,s]=maple(['evalf(Int(' ftn ',x=0..infinity))']);
int_val=sym(r);

I have not had any success with the command "quad" in Matlab, and so I
resorted to using the evalf(Int(*)) in Maple but with Matlab. Any
suggestions on speeding up the implementation?

Thank you again,

Marek
```
```On Jul 27, 2:48&#2013266080;pm, Clay <c...@claysturner.com> wrote:
> On Jul 27, 10:00&#2013266080;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 &#2013266080;(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

```
```On Jul 27, 3:23&#2013266080;pm, Clay <c...@claysturner.com> wrote:
> On Jul 27, 2:48&#2013266080;pm, Clay <c...@claysturner.com> wrote:
>
>
>
>
>
> > On Jul 27, 10:00&#2013266080;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 &#2013266080;(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

```
```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
```
```"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
```
```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

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

```
```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
```
```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
```
```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
```
```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