La ecuación diferencial de segundo orden

${\displaystyle x^2\frac{d^{2}y}{d{x}^2}+x\frac{dy}{dx}+\left(x^2-n^2\right)y=0}$

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

${\displaystyle 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 J_n(x) y Y_n(x) para valores enteros n=0,1,2,3....

${\displaystyle \begin{array}{l}{J}_n(x)={\displaystyle \sum _{k=0}^{\infty }\frac{{\left(-1\right)}^k{\left(x/2\right)}^{n+2k}}{k!\left(n+k\right)!}}\\ Y_n(x)=\frac{J_n(x)\cos \left(nx\right)-J_{-n}(x)}{\sin \left(nx\right)}\end{array}}$

Para n entero se cumple

${\displaystyle \begin{array}{l}{J}_{-n}(x)={\left(-1\right)}^n{J}_n(x)\\ Y_{-n}(x)={\left(-1\right)}^n{Y}_n(x)\end{array}}$

Funciones de Bessel de primera especie

Representamos J₀(x), J₁(x), J₂(x) y J₃(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 J₀(x), J₁(x) y J₂(x). Utilizamos la función raices para calcular las raíces múltiples de una ecuación transcendente

function bessel_2
    x=linspace(0,30,30);
    f=@(x) besselj(0,x);
    rr=raices(f,x); 
    disp('J_0(x)')
    disp( rr)
    
    f=@(x) besselj(1,x);
     rr=raices(f,x); 
    disp('J_1(x)')
    disp (rr)

    f=@(x) besselj(2,x);
    rr=raices(f,x); 
    disp('J_2(x)')
    disp(rr)

    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
end

Aproximaciones

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

${\displaystyle \begin{array}{l}{J}_0(x)\cong \sqrt{\frac{2}{\pi x}}\cos \left(x-\frac{\pi }{4}\right)\\ J_n(x)\cong \sqrt{\frac{2}{\pi x}}\cos \left(x-\frac{n\pi }{2}-\frac{\pi }{4}\right)\end{array}}$

Representamos la función J_n(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 J_n(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

${\displaystyle {\displaystyle \underset{0}{\overset{1}{\int }}x{J}_n(k_{m}x)J_n(k_{l}x)dx}=\{\begin{array}{l}\mathrm{0,}m\ne l\\ \frac{1}{2}{J}_{n+1}^2(k_m),m=l\end{array}}$

donde k_m y k_l son raíces distintas de J_n(k)=0. Por ejemplo, para n=0

${\displaystyle \underset{0}{\overset{1}{{\displaystyle \int }}}x{J}_0(k_{m}x)J_0(k_{l}x)dx=\{\begin{array}{l}\mathrm{0,}m\ne l\\ \frac{1}{2}{J}_1^2(k_m),m=l\end{array}}$

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

${\displaystyle {\displaystyle \underset{0}{\overset{1}{\int }}{x}^{n+1}{J}_n(kx)dx=\frac{1}{k}}{J}_{n+1}(k)}$

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

Derivadas

${\displaystyle \begin{array}{l}\frac{d{J}_n(k)}{dk}=J_{n-1}(k)-\frac{n}{k}{J}_n(k)\\ \frac{d{J}_n(k)}{dk}=\frac{n}{k}{J}_n(k)-J_{n+1}(k)\end{array}}$

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

${\displaystyle f(r)={\displaystyle \sum _{n=1}^{\infty }{a}_n{J}_0\left(k_n\frac{r}{R}\right)}}$

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

Para determinar los coeficientes a_n se utiliza las relaciones de ortogonalidad. Multiplicamos ambos miembros por J₀(k_nr/R) e integramos entre 0 y R

${\displaystyle {\displaystyle \underset{0}{\overset{R}{\int }}r·f(r)J_0\left(k_m\frac{r}{R}\right)}dr={\displaystyle \sum _{n=1}^{\infty }{a}_n{\displaystyle \underset{0}{\overset{R}{\int }}r·J_0\left(k_m\frac{r}{R}\right)J_0\left(k_n\frac{r}{R}\right)dr}}}$

Cuando m≠n la integral se hace cero, resultando

${\displaystyle \begin{array}{l}{\displaystyle \underset{0}{\overset{R}{\int }}r·f(r)J_0\left(k_n\frac{r}{R}\right)}dr=a_n{\displaystyle \underset{0}{\overset{R}{\int }}r·J_0^2\left(k_n\frac{r}{R}\right)dr}\\ a_n=\frac{2}{R^2{J}_1^2\left(k_m\right)}{\displaystyle \underset{0}{\overset{R}{\int }}r·f(r)J_0\left(k_m\frac{r}{R}\right)}dr\end{array}}$

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

Calculamos los coeficientes a_n

${\displaystyle a_n=\frac{2}{R^2{J}_1^2\left(k_n\right)}{\displaystyle \underset{0}{\overset{R}{\int }}r\left(1-\left|r-1\right|\right)J_0\left(k_n\frac{r}{R}\right)}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

Consideremos la función

${\displaystyle f(r)=\{\begin{array}{l}10\le r<b\\ 0b\le r\le R\end{array}}$

Desarrollo en serie en términos de J₁(x)

${\displaystyle f(r)={\displaystyle \sum _{n=1}^{\infty }{a}_n{J}_1\left(k_n\frac{r}{R}\right)}}$

Los coeficientes se calculan del siguiente modo

${\displaystyle a_n=\frac{2}{R^2{J}_2^2(k_n)}{\displaystyle \underset{0}{\overset{R}{\int }}r·f(r)J_1\left(k_n\frac{r}{R}\right)dr}=\frac{2}{R^2{J}_2^2(k_n)}{\displaystyle \underset{0}{\overset{b}{\int }}r·J_1\left(k_n\frac{r}{R}\right)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 Y₀(x),Y₁(x),Y₂(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 Y_n(x), tiende hacia

${\displaystyle \begin{array}{l}{Y}_0(x)\cong \sqrt{\frac{2}{\pi x}}\sin \left(x-\frac{\pi }{4}\right)\\ Y_n(x)\cong \sqrt{\frac{2}{\pi x}}\sin \left(x-\frac{n\pi }{2}-\frac{\pi }{4}\right)\end{array}}$

Derivadas

${\displaystyle \frac{d{Y}_n(x)}{dx}=n\frac{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 I_n y K_n

Son soluciones de la ecuación diferencial

${\displaystyle \begin{array}{l}{x}^2\frac{d^{2}y}{d{x}^2}+x\frac{dy}{dx}-(x^2+n^2)y=0\\ y=C{I}_n(x)+D{K}_n(x),x>0\end{array}}$

Representamos I₀(x), I₁(x), I₂(x) y I₃(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 K₀(x), K₁(x), K₂(x) y K₃(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

${\displaystyle \begin{array}{l}\frac{d{I}_n(x)}{dx}=I_{n+1}(x)+\frac{n}{x}{I}_n(x),\frac{d{I}_n(x)}{dx}=I_{n-1}(x)-\frac{n}{x}{I}_n(x)\\ \frac{d{K}_n(x)}{dx}=-K_{n+1}(x)+\frac{n}{x}{K}_n(x),\frac{d{K}_n(x)}{dx}=-K_{n-1}(x)-\frac{n}{x}{K}_n(x)\end{array}}$

>> 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, j_n(x) y y_n(x) son las dos soluciones independientes de la ecuación diferencial

${\displaystyle x^2\frac{d^{2}y}{d{x}^2}+2x\frac{dy}{dx}+\left(x^2-n(n+1)\right)y=0}$

donde n es un número entero

${\displaystyle 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, j_n(x) y y_n(x), están relacionadas con las funciones de Bessel J_n(x) e Y_n(x)

${\displaystyle j_n(x)=\sqrt{\frac{\pi }{2x}}{J}_{n+\frac{1}{2}}(x),y_n(x)=\sqrt{\frac{\pi }{2x}}{Y}_{n+\frac{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

• Representamos j_n(x) para n=0,1,2,3

hold on
for l=0:3
    f=@(x) sqrt(pi./(2*x)).*besselj(l+1/2,x);
    fplot(f,[0,14])
end
xlabel('x')
ylabel('j_n(x)')
legend('n=0','n=1','n=2','n=3','Location', 'best')
grid on
title('Funciones esféricas de Bessel')



• Representamos y_n(x) para n=0,1,2,3

hold on
for l=0:3
    f=@(x) sqrt(pi./(2*x)).*bessely(l+1/2,x);
    fplot(f,[1,14])
end
ylim([-0.4,0.4])
xlabel('x')
ylabel('y_n(x)')
legend('n=0','n=1','n=2','n=3','Location', 'best')
grid on
title('Funciones esféricas de Neumann')

Funciones esféricas de Hankel

Las funciones esféricas de Hankel son

${\displaystyle \begin{array}{l}{h}_n^{(1)}=j_n(x)+i{y}_n(x)\\ h_n^{(2)}=j_n(x)-i{y}_n(x)\end{array}}$

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

${\displaystyle \begin{array}{l}{h}_{-1}^{(1)}(x)=\frac{1}{x}{e}^{ix}\\ h_0^{(1)}(x)=-i\frac{1}{x}{e}^{ix}\\ h_1^{(1)}(x)=-\frac{x+i}{x^2}{e}^{ix}\\ h_2^{(1)}(x)=i\frac{x^2+3ix-3}{x^3}{e}^{ix}\\ h_3^{(1)}(x)=\frac{x^3+6i{x}^2-15x-15i}{x^4}{e}^{ix}\end{array}}$

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

Condensador concectado a una fuente de corriente alterna

Ondas guiadas. Guía de onda circular