Funciones de Bessel

La ecuación diferencial de segundo orden

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

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

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

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

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

Para n entero se cumple

J n (x)= ( 1 ) n J n (x) Y n (x)= ( 1 ) n Y n (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], 'displayName',num2str(n));
end
legend('-DynamicLegend','location','Best')
xlabel('x')
ylabel('J(x)')
title('Funciones J_n(x) de Bessel')
grid on
hold off

En la siguiente tabla, se proporcionan los primeros ceros de las funciones de Bessel J0(x), J1(x) y J2(x). Utilizamos la función raices para calcular las raíces múltiples de una ecuación transcendente

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

Aproximaciones

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

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

Representamos la función Jn(x) y su aproximación asintótica

n=3;
f=@(x) sqrt(2./(pi*x)).*cos(x-n*pi/2-pi/4);
k=0:4;
hold on
fplot(@(x) besselj(n,x),[1,20])
fplot(f,[1,20])
r=(2*k+1)*pi/2+n*pi/2+pi/4; %ceros coseno
plot(r,0,'o','markersize',3,'markeredgecolor','r','markerfacecolor','r')
hold off
legend('J_n(x)','cos(x)')
xlabel('x')
ylabel('J_n(x)')
title('Funciones J_n(x) de Bessel')
grid on
hold off

En la figura, se ha representado los ceros de la función cos(x-nπ/2-π/4), que como vemos son próximos a los ceros de la función Jn(x). Este hecho nos permite calcular los ceros de la función de Bessel mediante el comando fzero de MATLAB. Recuérdese que los ceros de la función coseno son (2k+π)/2, k=0,1,2,3...

n=3;
for k=0:5
    r=(2*k+1)*pi/2+n*pi/2+pi/4; %ceros coseno
    x1=fzero(@(x) besselj(n,x), r);
    disp(x1)
end
    6.3802
    9.7610
   13.0152
   16.2235
   19.4094
   22.5827

Relaciones de ortogonalidad

0 1 x J n ( k m x) J n ( k l x)dx ={ 0,ml 1 2 J n+1 2 ( k m ),m=l

donde km y kl son raíces distintas de Jn(k)=0. Por ejemplo, para n=0

0 1 x J 0 ( k m x) J 0 ( k l x)dx={ 0,ml 1 2 J 1 2 ( k m ),m=l

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)); %raíces distintas
integral(f,0,1)
ans =   8.9311e-18
x=linspace(0,40,40);
f=@(x) besselj(0,x);
k=raices(f,x); 

f=@(x) x.*(besselj(0,k(2)*x).^2); %segunda raíz de J0
integral(f,0,1)
ans =    0.0579
>> besselj(1,k(2))^2/2
ans =    0.0579

Otras relaciones

0 1 x n+1 J n (ax)dx= 1 a J n+1 (a)

para n>-1. Por ejemplo, para n=0

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

x J n1 (x)+x J n+1 (x)=2n J n (x) J n1 (x) J n+1 (x)=2 d J n (x) dx

De la que se deduce

d J n (x) dx =n J n (x) x J n+1 (x)

>> syms x n;
>> diff(besselj(n,x))
ans =(n*besselj(n, x))/x - besselj(n + 1, x)

Referencia: I.S. Gradshteyn, I.M. Ryzhik. Table of Integrals, Series, and Products. 2007, Elsevier Inc., apartado 6.561, n° 5, pág. 676, apartado 8.472, n° 1, 2, pág. 926

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 n r R ) dr= a n 0 R r· J 0 2 ( k n r R )dr a n = 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:length(k)
    y=y+a(i)*besselj(0,k(i)*r/R);
end
hold on
fplot(@(r) 1-abs(r-1),[0,R]);
plot(r,y);
hold off
grid on
xlabel('r')
ylabel('f(r)')
title('Serie Fourier-Bessel')

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

Desarrollo en serie en términos de J1(x)

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

Los coeficientes se calculan del siguiente modo

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

x=linspace(0,40,40);
f=@(x) besselj(1,x);
k=raices(f,x); 
 
R=2;
b=1;
a=zeros(1,length(k));
for i=1:length(k)
    f=@(r) r.*besselj(1,k(i)*r/R);
    a(i)=2*integral(f,0,b)/(R*besselj(2,k(i)))^2;
end
r=linspace(0,R,200);
y=zeros(1,length(r));
for i=1:length(k)
    y=y+a(i)*besselj(1,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')

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,

hold on
for n=0:4
    fplot(@(x) bessely(n,x),[0.02,15],'displayName',num2str(n));
end
ylim([-2,1])
legend('-DynamicLegend','location','Best')
xlabel('x')
ylabel('Y(x)')
title('Funciones Y_n(x) de Bessel')
grid on
hold off

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

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

Derivadas

d Y n (x) dx =n Y n (x) x Y n+1 (x)

>> syms x n;
>> diff(bessely(n,x))
ans =(n*bessely(n, x))/x - bessely(n + 1, x)

Funciones de Bessel modificadas In y Kn

Son soluciones de la ecuación diferencial

x 2 d 2 y d x 2 +x dy dx ( x 2 + n 2 )y=0 y=C I n (x)+D K n (x),x>0

Representamos I0(x), I1(x), I2(x) y I3(x) llamando a la función besseli(n,x) de MATLAB, dando valores a su primer parámetro n=0,1,2,3,

hold on
for n=0:3
   fplot(@(x) besseli(n,x),[0.1,5])
end
hold off
ylim([0,5])
grid on
legend('I_0','I_1','I_2','I_3', 'location','NorthWest')
xlabel('x')
ylabel('I_n(x)')
title('Funciones de Bessel, modificadas')

Representamos K0(x), K1(x), K2(x) y K3(x) llamando a la función besselk(n,x) de MATLAB, dando valores a su primer parámetro n=0,1,2,3,

hold on
for n=0:3
   fplot(@(x) besselk(n,x),[0.1,5])
end
hold off
ylim([0,5])
grid on
legend('K_0','K_1','K_2','K_3')
xlabel('x')
ylabel('K_n(x)')
title('Funciones de Bessel, modificadas')

Derivadas

d I n (x) dx = I n+1 (x)+ n x I n (x), d I n (x) dx = I n1 (x) n x I n (x) d K n (x) dx = K n+1 (x)+ n x K n (x), d K n (x) dx = K n1 (x) n x K n (x)

>> syms x n;
>> diff(besseli(n,x))
ans =besseli(n + 1, x) + (n*besseli(n, x))/x
>> diff(besselk(n,x))
ans =(n*besselk(n, x))/x - besselk(n + 1, x)

Funciones esféricas de Bessel

Las funciones esféricas de Bessel, jn(x) y yn(x) son las dos soluciones independientes de la ecuación diferencial

x 2 d 2 y d x 2 +2x dy dx +( x 2 n(n+1) )y=0

donde n es un número entero

y(x)=A· j n (x)+B· y n (x)

A y B son coeficientes a determinar a partir de las condciones de contorno

Las funciones esféricas de Bessel, jn(x) y yn(x), están relacionadas con las funciones de Bessel Jn(x) e Yn(x)

j n (x)= π 2x J n+ 1 2 (x), y n (x)= π 2x Y n+ 1 2 (x)

Las primeras funciones esféricas de Bessel son

syms x;
for n=-1:3
    jn=sqrt(pi/(2*x))*besselj(n+1/2,x);
    yn=sqrt(pi/(2*x))*bessely(n+1/2,x);
    disp([simplify(jn),simplify(yn)])
end
jn(x)yn(x)
j 1 (x)= cosx x y 1 (x)= sinx x
j 0 (x)= sinx x y 0 (x)= cosx x
j 1 (x)= sinx x 2 cosx x y 1 (x)= cosx x 2 sinx x
j 2 (x)=( 3 x 2 1 ) sinx x 3cosx x 2 y 2 (x)=( 3 x 2 +1 ) cosx x 3sinx x 2
j 3 (x)=( 15 x 3 6 x ) sinx x ( 15 x 2 1 ) cosx x y 3 (x)=( 15 x 3 + 6 x ) cosx x ( 15 x 2 1 ) sinx x

Funciones esféricas de Hankel

Las funciones esféricas de Hankel son

h n (1) = j n (x)+i y n (x) h n (2) = j n (x)i y n (x)

Las primeras funciones h n (1) son

syms x;
for n=-1:3
    hn=sqrt(pi/(2*x))*besselh(n+1/2,1,x);
    disp(simplify(hn))
end

h 1 (1) (x)= 1 x e ix h 0 (1) (x)=i 1 x e ix h 1 (1) (x)= x+i x 2 e ix h 2 (1) (x)=i x 2 +3ix3 x 3 e ix h 3 (1) (x)= x 3 +6i x 2 15x15i x 4 e ix

Ejemplos en el curso de Física

La ecuación de Laplace, coordenadas cilíndricas (I)

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)

La ecuación de Schrödinger en coordenadas cilíndricas

La ecuación de Schrödinger en coordenadas esféricas (I)

La ecuación de Schrödinger en coordenadas esféricas (II)

Movimiento vertical de un cohete con rozamiento (I)

Movimiento vertical de un cohete con rozamiento (II)

Ecuaciones de Navier-Stokes. Estado transitorio (II)

Fluido en un cilindro en rotación, estado transitorio