# How to evaluate integral from 0 to inf of besselj(0,kr) * besselj(0,kR) * 1/k * (2 - e^-kz - ek(z-L)) dk

19 views (last 30 days)

Show older comments

Hello,

I am trying to evaluate the integral given in equation 6 in this paper (see also eq 10 for g()):

"Coulomb potential and energy of a uniformly charged cylindrical shell"

I've tried it symbolically with Python Sympy library, and the code never finished running. I tried numerically with scipy.integrate.quad , but the margin of error in the answer was huge, 1/4 of the answer. Finally I tried in MATLAB, and got the same errors that Python was giving.

syms r R z L k

f = @(k) besselj(0, r * k) * besselj(0, R * k) * (1/k) * (2 - exp(-k*z) - exp(k*(z-L)))

f =

function_handle with value:

@(k)besselj(0,r*k)*besselj(0,R*k)*(1/k)*(2-exp(-k*z)-exp(k*(z-L)))

int(f,k,0,inf)

ans =

int(-(besselj(0, R*k)*besselj(0, k*r)*(exp(-k*z) + exp(-k*(L - z)) - 2))/k, k, 0, Inf)

Where:

k = Variable of Integration

R = Radius of Charged Cylinder

L = Length of Cylinder

r = r-coordinate of Point at which the Potential is being measured, (In Cylindrical Coordinates)

z = z-coordinate of Point

Is it possible to get a symbolic solution? If not, how could I get a numerical solution, if I specied values for R, L, r, and z?

Thank you!

##### 2 Comments

cr
on 19 Dec 2022

### Accepted Answer

Fabio Freschi
on 19 Dec 2022

Edited: Fabio Freschi
on 19 Dec 2022

For the numerical integration, you can use integral, that also accepts Inf as integration upper bound

clear variables, close all

% some randoms values for the params

r = 1;

R = 2;

z = 1.5;

L = 2;

% function handle

f = @(k)besselj(0,r*k).*besselj(0,R*k).*(1./k).*(2-exp(-k*z)-exp(k*(z-L)));

% indefinite integral

Vinf = integral(f,0,Inf);

One could worry about the warning message. However, without playing with the tolerances and other optional inputs of integral, it seems that the result is reasonable:

% definite integral with increasing upper bound

N = 100;

V = zeros(N,1);

for i = 1:N

V(i) = integral(f,0,max([r,R,z,L])*i);

end

figure, hold on

plot(1:N,V)

plot([1 N],[Vinf,Vinf],'--')

legend('V','Vinf')

##### 2 Comments

Fabio Freschi
on 22 Dec 2022

@Sean: my pleasure! if this answer solves your original problem, please accept it

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!