Funciones de Bessel

La ecuación diferencial de segundo orden

x 2 d 2 y d x 2 +x dy dx +( x 2 ν 2 )y=0

se conoce como ecuación de Bessel. La solución de esta ecuación diferencial se escribe

y=A J ν (x)+B Y ν (x)

donde A y B son constantes que se determinan a partir de las condiciones iniciales. El índice ν es un número real, aunque en los libros se proporcionan tablas de las funciones Jν(x) y Yν(x) para valores enteros ν=0,1,2,3....

J ν (x)= k=0 ( 1 ) k ( x/2 ) ν+2k k!( ν+k )! Y ν (x)= J ν (x)cos( νx ) J ν (x) sin( νx )

Para ν entero se cumple

J ν (x)= ( 1 ) ν J ν (x) Y ν (x)= ( 1 ) ν Y ν (x)

Funciones de Bessel de primera especie

Representamos J0(x), J1(x), J2(x) y J3(x) llamando a la función besselj(n,x) de MATLAB, dando valores a su primer parámetro n=0,1,2,3,

hold on
for n=0:3
    f=@(x) besselj(n,x);
    fplot(f,[0,20]);
end
xlabel('x')
ylabel('J(x)')
title('Funciones J_\nu(x) de Bessel')
grid on
hold off

Cuando x se hace grande la función Jν(x), tiende hacia

J 0 (x) 2 πx cos( x π 4 ) J ν (x) 2 πx cos( x νπ 2 π 4 )

En la siguiente tabla, se proporcionan los primeros ceros de las funciones de Bessel J0(x), J1(x) y J2(x).

function r = raices(f, x)
    y=f(x);
    indices=find(y(1:end-1).*y(2:end)<0);
    r=zeros(1,length(indices));
    for k=1:length(indices)
        r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
    end
end
x=linspace(0,40,40);
f=@(x) besselj(0,x);
k=raices(f,x); 
disp('J_0(x)')
disp(k)

f=@(x) besselj(1,x);
k=raices(f,x); 
disp('J_1(x)')
disp(k)

f=@(x) besselj(2,x);
k=raices(f,x); 
disp('J_2(x)')
disp(k)
J0(x) J1(x) J2(x)
2.4048 3.8317 5.1356
5.5201 7.0156 8.4172
8.6537 10.1735 11.6198
11.7915 13.3237 14.7960
14.9309 16.4706 17.9598
18.0711 19.6159 21.1170
21.2116 22.7601 24.2701
24.3525 25.9037 27.4206
27.4935 29.0468 30.5692
30.6346 32.1897 33.7165

Relaciones de ortogonalidad

Comprobamos que la integral

0 1 x J 0 ( k m x) J 0 ( k n x)dx=0mn

donde km y kn son raíces distintas de J0(k)=0

x=linspace(0,40,40);
f=@(x) besselj(0,x);
k=raices(f,x); 

f=@(x) x.*(besselj(0,k(1)*x).*besselj(0,k(2)*x));
integral(f,0,1)
ans =   8.9311e-18

Cuando m=n

0 1 x J 0 2 (kx) dx= 1 2 ( J 0 2 (k)+ J 1 2 (k) )

>> syms x k;
>> int('x*besselj(0,k*x)^2',x,0,1)
ans= besselj(0, k)^2/2 + besselj(1, k)^2/2

Otras relaciones

0 1 x J 0 (kx)dx= 1 k J 1 (k)

>> syms k x;
>> int('x*besselj(0,k*x)',x,0,1)
ans =besselj(1, k)/k

Derivadas

d J 0 (kx) dx =k J 1 (kx)

>> syms x k;
>> diff(besselj(0,k*x))
ans =-k*besselj(1, k*x)

Desarrollo en serie Fourier-Bessel

En este apartado explicamos el caso más simple, el de una función f(r) tal que f(R)=0 entonces se puede expresar en términos de un desarrollo en serie de Fourier-Bessel

f(r)= n=1 a n J 0 ( k n r R )

donde kn son los ceros de la función J0(x) o las raíces de la ecuación transcendente J0(k)=0

Para determinar los coeficientes an se utiliza las relaciones de ortogonalidad. Multiplicamos ambos miembros por J0(knr/R) e integramos entre 0 y R

0 R r·f(r) J 0 ( k m r R ) dr= n=1 a n 0 R r· J 0 ( k m r R ) J 0 ( k n r R )dr

Cuando m≠n la integral se hace cero, resultando

0 R r·f(r) J 0 ( k m r R ) dr= a m 0 R r· J 0 2 ( k m r R )dr a m = 2 R 2 J 1 2 ( k m ) 0 R r·f(r) J 0 ( k m r R ) dr

Ejemplo 1

Consideremos la función de la figura, f(r)=1-|r-1| 0≤rR=2

Calculamos los coeficientes an

a n = 2 R 2 J 1 2 ( k n ) 0 R r( 1| r1 | ) J 0 ( k n r R ) dr

Integramos numéricamente mediante la función integral de MATLAB

x=linspace(0,40,40);
f=@(x) besselj(0,x);
k=raices(f,x); 
 
R=2;
a=zeros(1,length(k));
for i=1:length(k)
    f=@(r) (r.*(1-abs(r-1))).*besselj(0,k(i)*r/R);
    a(i)=2*integral(f,0,R)/(R*besselj(1,k(i)))^2;
end
r=linspace(0,R,200);
y=zeros(1,length(r));
for i=1:10
    y=y+a(i)*besselj(0,k(i)*r/R);
end
hold on
fplot('1-abs(r-1)',[0,R]);
plot(r,y);
hold off
grid on
xlabel('r')
ylabel('f(r)')
title('Serie Fourier-Bessel')

Nos aproximamos con los 10 primeros términos del desarrollo en serie a la función f(r). Los valores de los primeros cinco coeficientes son los siguientes:

>> a(1:5)
ans =    1.0149   -0.6448   -0.3029    0.0365    0.0865
Ejemplo 2

Consideremos la función

f(r)={ 10r<b 0brR

Los coeficientes se calculan del siguiente modo

a n = 2 R 2 J 1 2 ( k n ) 0 R r·f(r) J 0 ( k n r R ) dr= 2 R 2 J 1 2 ( k n ) 0 b r· J 0 ( k n r R ) dr

x=linspace(0,40,40);
f=@(x) besselj(0,x);
k=raices(f,x); 
 
R=2;
b=1;
a=zeros(1,length(k));
for i=1:length(k)
    f=@(r) r.*besselj(0,k(i)*r/R);
    a(i)=2*integral(f,0,b)/(R*besselj(1,k(i)))^2;
end
r=linspace(0,R,200);
y=zeros(1,length(r));
for i=1:10
    y=y+a(i)*besselj(0,k(i)*r/R);
end
yy=r<=b;
hold on
plot(r,yy);
plot(r,y);
hold off
grid on
xlabel('r')
ylabel('f(r)')
title('Serie Fourier-Bessel')

Como podemos apreciar, son necesarios bastantes términos del desarrollo en serie para aproximarnos a la función f(r). Los valores de los primeros cinco coeficientes an son los siguientes

>> a(1:5)
ans =    0.7698    0.6615   -0.2830   -0.4643    0.1987

Funciones de Bessel de segunda especie

Representamos Y0(x),Y1(x),Y2(x) llamando a la función bessely(n,x) de MATLAB, dando valores a su primer parámetro n=0,1,2,

color=['b','r','g'];
x=0.02:0.02:15;
hold on
for n=0:length(color)-1
    y=bessely(n,x);
    plot(x,y,color(n+1),'displayName',num2str(n));
end
ylim([-2,1])
legend('-DynamicLegend','location','Best')
xlabel('x')
ylabel('Y(x)')
title('Funciones Y_\nu(x) de Bessel')
grid on
hold off

Cuando x se hace grande la función Yν(x), tiende hacia

Y 0 (x) 2 πx sin( x π 4 ) Y ν (x) 2 πx sin( x νπ 2 π 4 )

Ejemplos en el curso de Física

El problema de Kepler

Difracción Fraunhofer producida por una abertura circular

Modos normales de un cable en movimiento de rotación, suspendido verticalmente.

Modos de vibración de una membrana circular

Conducción del calor en un cilindro muy largo (I)

Conducción del calor en un cilindro muy largo (II)