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
Approximation to Bessel Function in Integrand
Started by ●July 27, 2009
Reply by ●July 27, 20092009-07-27
Reply by ●July 27, 20092009-07-27
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
Reply by ●July 27, 20092009-07-27
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
Reply by ●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
Reply by ●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 ●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, > > MarekHello 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 ●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 ●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.466Mathematica'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 ●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 > > ClayAnd 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 ●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, > ClaySee 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
