-2
$\begingroup$

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.enter image description here

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?

$\endgroup$

closed as off-topic by Gilles, MBaz, Jason R, Peter K. Sep 12 '16 at 12:07

This question appears to be off-topic. The users who voted to close gave these specific reasons:

  • "General programming questions are off-topic here, but can be asked on Stack Overflow." – MBaz, Peter K.
  • "This question does not appear to be about signal processing within the scope defined in the help center." – Gilles, Jason R
If this question can be reworded to fit the rules in the help center, please edit the question.

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

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')

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Why the downvote? $\endgroup$ – fibonatic Sep 10 '16 at 13:02
  • 1
    $\begingroup$ @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). $\endgroup$ – fibonatic Sep 10 '16 at 16:16
2
$\begingroup$

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);

enter image description here

$\endgroup$
  • $\begingroup$ thanks, but it is not working fine with matlab! $\endgroup$ – iamgr007 Sep 10 '16 at 15:34
  • 1
    $\begingroup$ @AIaye: Can only be a very minor thing, but anyway, you also got another solution. $\endgroup$ – Matt L. Sep 10 '16 at 15:36

Not the answer you're looking for? Browse other questions tagged or ask your own question.