Integration of a function contains bessel functions from 0 to inf
9 views (last 30 days)
Show older comments
Ghs Jahanmir
on 31 Aug 2021
Commented: Ghs Jahanmir
on 5 Sep 2021
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*-')
0 Comments
Accepted Answer
Walter Roberson
on 31 Aug 2021
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)
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))
F = int(myfunc, 0, inf);
val=1+(2/Pi).*F ;
valY = subs(val, x, X.'/R)
valYn = double(valY)
plot(X, valYn, 'g*-')
3 Comments
Walter Roberson
on 31 Aug 2021
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.
More Answers (0)
See Also
Categories
Find more on Calculus in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!