Search code examples
matlabbessel-functions

Plotting and finding roots of bessel functions


I am trying to plot roots of a function that is composed of multiple bessel functions being added and multiplied in Matlab. The equation is Jm(omega)*Ik(omega)+Im(omega)*Jk(omega) where Jm is the bessel function of the first kind of order m (besselj). Im is the modified bessel function of the first kind of order m (besseli). For each mode m=o,1,2,...and n=1,2,3... The frequency omega(mn) is the corresponding root of the listed equation. m=0,1,2 n-1,2,3,4. I need to solve the equation for the 12 roots. I am new to Matlab and this is a little out of my league. So far I have this code but I wasn't sure if I needed the variable omega in the script or not. I have also looked at other people's questions on the suject but didn't see any quite like this. The plots I have seen look nothing like mine which tells me I am probably wrong. Thanks for any help.

m=(0:2); k=(1:3); n=(1:4);
Jm=besselj(m,n');
Ik=besseli(k,n');
Jk=besselj(k,n');
Im=besseli(m,n');
g=Jm.*Ik+Im.*Jk
plot(g)

Solution

  • Plotting

    besselj and besseli take what you call omega as their second parameter, so to plot your function you should try something like

    m=0; k=1; omega=0:0.02:10;
    Jm=besselj(m,omega);
    Ik=besseli(k,omega);
    Jk=besselj(k,omega);
    Im=besseli(m,omega);
    g=Jm.*Ik+Im.*Jk;
    plot(omega,g);
    hold all;
    plot(omega,0,'k');
    axis([min(omega) max(omega) -100 100]);
    

    This shows you that for m=1, k=1 the first zeros are around 3.2, 6.3 and 9.4:

    screenshot matlab figure

    Finding the roots numerically

    You could implement Halley's method for your function g, similar to how the roots of besselj are determined in the MatlabCentral file linked by Cheery.