# Bessel's integral matlab implementation [closed]

I'm trying to implement this integral representation of Bessel function of the first kind of order n.

$$J_n(x)=\frac{1}{2\pi} \int_{-\pi}^{\pi}e^{i(n\tau-x\sin(\tau))}d\tau$$

here is what I tried:

t = -pi:0.1:pi;
n = 1;
x = 0:5:20;
A = @(t) exp(sqrt(-1)*(n*t-x*sin(t)));
B = integral(A(t),-pi,pi);
plot(A(t),x)


the plot i'm trying to get is as shown in the wikipedia page. it said:

Error using * Inner matrix dimensions must agree.

Error in besselfn (line 8) A(t) = exp(sqrt(-1)*(n*t-x*sin(t)));

so i tried putting x=5;

and the output was:

Subscript indices must either be real positive integers or logicals.

Error in besselfn (line 8) A(t) = exp(sqrt(-1)*(n*t-x*sin(t)));

How to get this correct? what am I missing?

• Why don't you just use besselj? – Matt L. Sep 10 '16 at 8:10
• @MattL. I want to know how they implemented it! – iamgr007 Sep 10 '16 at 9:29

Like Matt has said, you have to evaluate the integral for each value of x separately, which can be done with a for loop.

And the other remark is that you can have to pass a function as the first argument to the function integral. So instead of A(t), which will evaluate your function at the values you declared for t (which always returns the same vector), just use A (or @(t) A(t)).

One way you could implement this is:

n = 0 : 2;
x = linspace(0, 20, 200);

A = @(t, x, n) exp(1i * (n * t - x * sin(t)));
B = zeros(length(n), length(x));

for m = 1 : length(n)
for k = 1 : length(x)
B(m,k) = integral(@(t) A(t, x(k), n(m)), -pi, pi);
end
end

plot(x, B(1,:), 'r', x, B(2,:), 'g--', x, B(3,:), 'b-.')
legend('J_0(x)', 'J_1(x)', 'J_2(x)', 'Location', 'NorthEast') • Why the downvote? – fibonatic Sep 10 '16 at 13:02
• @AIaye do not that MATLAB will complain about B being complex, if you do not want that you could use the cosine instead of exp (like in Matt's answer). Also the graphs in your question are normalized by a factor 1/(2pi). – fibonatic Sep 10 '16 at 16:16

You indeed have to use just a single value for the variables x and n that you pass to your integration function. If you want to evaluate the integral for different values of x or n you need to loop over these values.

Below is an example using octave (with quad instead of integral, but I think it should also work in Matlab using integral). If you don't understand the use of the function handle below, check out Parameterizing Using Anonymous Functions.

function J = mybesselj(n,x)

J = quad(@(t) cos(n*t-x*sin(t)), -pi, pi) / (2*pi);


[Note that I used the function $f(t)=\cos(nt-x\sin(t))$ instead of the complex exponential, because the imaginary part of the integrand $\exp(j(nt-\sin(t)))$ is odd, and, consequently, the integral over the imaginary part vanishes (due to the symmetrical integration limits).]

Now you have to use a loop to compute the Bessel function for a range of values of x and n:

x = linspace(0,20,100);
J0 = zeros(100,1);
J1 = zeros(100,1);
J2 = zeros(100,1);
for i=1:100,
xi = x(i);
J0(i) = mybesselj(0,xi);
J1(i) = mybesselj(1,xi);
J2(i) = mybesselj(2,xi);
end
plot(x,J0,x,J1,x,J2); • thanks, but it is not working fine with matlab! – iamgr007 Sep 10 '16 at 15:34
• @AIaye: Can only be a very minor thing, but anyway, you also got another solution. – Matt L. Sep 10 '16 at 15:36