Integration of a function contains bessel functions from 0 to inf

Dear All,

I am trying to calculate the following formula using "vpaintegral".

(I used this post as my reference https://www.mathworks.com/matlabcentral/answers/709098-indefinite-integrals-of-bessel-function )

D=1e-11;
R=0.001
x=1:.1:6; %(a vector)
t=[1 5 7 24 30] %(one or several values)

x is argument of the below function.

I wrote the following code to evaluate the integral for just one single value of t. But it dose not give me the val as vector based on x values:

[ 0, vpaintegral(-(exp(-(4760145414732603*x^2)/18889465931478580854784)*(besselj(0, x/1000)*bessely(0, (3*x)/2000) - bessely(0, x/1000)*besselj(0, (3*x)/2000)))/(x*(besselj(0, x/1000)^2 + bessely(0, x/1000)^2)), x, 0, Inf), vpaintegral((exp(-(4760145414732603*x^2)/18889465931478580854784)*(besselj(0, x/500)*bessely(0, x/1000) - bessely(0, x/500)*besselj(0, x/1000)))/(x*(besselj(0, x/1000)^2 + bessely(0, x/1000)^2)), x, 0, Inf), vpaintegral((exp(-(4760145414732603*x^2)/18889465931478580854784)*(besselj(0, x/400)*bessely(0, x/1000) - bessely(0, x/400)*besselj(0, x/1000)))/(x*(besselj(0, x/1000)^2 + bessely(0, x/1000)^2)), x, 0, Inf), vpaintegral(-(exp(-(4760145414732603*x^2)/18889465931478580854784)*(besselj(0, x/1000)*bessely(0, (3*x)/1000) - besselj(0, (3*x)/1000)*bessely(0, x/1000)))/(x*(besselj(0, x/1000)^2 + bessely(0, x/1000)^2)), x, 0, Inf)]

any help is appreciated.

Thanks in advance for your suggestions,

---------------------------------------------

syms x
f1=exp((-D*t)*x.^2)
b_y=bessely(0,x.*R);
b_j=besselj(0,x.*R);
yj=(bessely(0,x.*R).^2)+(besselj(0,x.*R).^2);
myfunc=f1.*((besselj(0,x.*r).*b_y)-(bessely(0,x.*r).*b_j))./(x.*yj);
F = vpaintegral(myfunc, 0, inf);
val=1+(2/pi).*F  ;
plot((x/R),val, 'g*-')

Answer 1

[Walter Roberson]

Q = @(v) sym(v);
D = Q(1e-11);
R = Q(0.001);
X = 1:.1:6; %(a vector)
t = [1 5 7 24 30]; %(one or several values)

syms x r
r = 3;  %maybe ???
Pi = sym(pi);
f1=exp((-D*t)*x.^2)

f1 = 

b_y=bessely(0,x.*R);
b_j=besselj(0,x.*R);
yj=(bessely(0,x.*R).^2)+(besselj(0,x.*R).^2);
myfunc = simplify(f1.*((besselj(0,x.*r).*b_y)-(bessely(0,x.*r).*b_j))./(x.*yj))

myfunc = 

F = int(myfunc, 0, inf);
val=1+(2/Pi).*F  ;
valY = subs(val, x, X.'/R)

valY = 

valYn = double(valY)
plot(X, valYn, 'g*-')

[Ghs Jahanmir]

@Walter Roberson

Thanks a lot for your detailed code. Actually I have a type in the problem definition.

instead of X x=1:.1:6; %(a vector)

it should be

r=1:.1:6; %(a vector).

x is the integral variable and I dont assign any range for that except integration over [0,inf].

[Walter Roberson]

The values increase rapidly near 0. For x = 10^-N then near -140-ish, the value of the function is roughly -10^(N-5) . That makes it difficult to integrate.

The integral from 10^-145 to 10^-100 is about 0.0139 -- so neglecting even remarkably small coefficients is enough to change the value in the second decimal place. That makes integration difficult. Even with lowering the abstol and reltol and increasing the maximum function evaluations, and putting boundaries like 10^-100 to 10^7, a single integral can take about 15 minutes, and you have an array of integrals to do.

My tests found that you still got relatively meaningful function values up to about x = 10^7 but by 10^8 it had dropped off sharply and was not work computing 10^7 to infinitity . Interestingly, it is towards the upper bound that really takes the long integration time -- 1e-100 to 1e4 took about 90 seconds, 1e4 to 1e7 took about 750 seconds. Almost all the large scale wiggles happen by x = 25 or so, but the small scale wiggles over that long of a period add up to make a difference in the overall integral.

[Ghs Jahanmir]

Dear @Walter Roberson

you are right. Thanks for your effort and test. I actually try to solve this to reproduce a graph in a book (for 3 times (t). y-axis is the "val" and x-axis is "r" for three "t" values.

I wonder how they were able to obtain this integration and drew the graph long time ago.

Best Regards,

Ghs