Integrales elípticas

Lemniscata de Bernoulli

La ecuación de la Lemniscata de Bernouilli en coordenadas polares es

r2=a2cos(2θ), -π/4≤θ≤π/4

a=1;
th=linspace(-pi/4,pi/4,200);
r=a*sqrt(cos(2*th));
x=r.*cos(th);
y=r.*sin(th);
hold on
plot(x,y, 'b')
plot(-x,y,'b')
hold off
axis equal
grid on
xlabel('x')
ylabel('y')
title('Lemniscata de Bernoulli')

La longitud de la lemniscata de Bernoulli es:

l=4 0 π/4 ds

Calculamos la longitud de arco ds de una porción infinitesimal de la lemniscata de Bernoulli

{ x=rcosθ y=rsinθ { dx=cosθ·drrsinθ·dθ dy=sinθ·dr+rcosθ·dθ d s 2 =d x 2 +d y 2 =d r 2 + r 2 d θ 2 r 2 = a 2 cos(2θ) rdr= a 2 sin(2θ)·dθ d s 2 = a 2 d θ 2 cos(2θ)

Integramos, respecto del ángulo, θ

l=4 0 π/4 ds =4a 0 π/4 dθ cos(2θ)

Hacemos el cambio de variable

{ cos(2θ)= cos 2 φ 2sin(2θ)·dθ=2cosφsinφ·dφ

La integral se convierte en

l=4a 0 π/2 sinφdφ 1 cos 4 φ =4a 0 π/2 dφ 1+ cos 2 φ =4a 0 π/2 dφ 2 sin 2 φ = 4a 2 0 π/2 dφ 1 1 2 sin 2 φ = 4a 2 K( 1 2 )

>> 4*ellipke(1/2)/sqrt(2)
ans =    5.2441

La integral elíptica completa de primera especie es

K(k)= 0 π/2 dφ 1 k 2 sin 2 φ

Una integral similar, la encontramos al calcular el periodo de un péndulo para cualquier amplitud θ0

Una generalización de esta integral

F(θ,k)= 0 θ dφ 1 k 2 sin 2 φ 0θ π 2

se denomina integral elíptica incompleta de primera especie.

En la página titulada 'Caída del extremo libre de una cadena doblada' se utilizan las integrales elípticas imcompletas

Longitud de una elipse

La ecuación de una elipse de semiejes a y b es

x 2 a 2 + y 2 b 2 =1

o en forma paramétrica
x=a·cosφ
y=b·sinφ

Como podemos ver a la derecha de la figura, la longitud de un elemento diferencial ds del arco de la elipse es

ds= d x 2 +d y 2 = a 2 sin 2 φ+ b 2 cos 2 φ ·dφ=b 1 b 2 a 2 b 2 sin 2 φ ·dφ =b 1 k 2 sin 2 φ ·dφ(b>a)

La longitud de la elipse es

L=b 0 2π 1 k 2 sin 2 φ ·dφ=4 b 0 π/2 1 k 2 sin 2 φ ·dφ

Habitualmente, el semieje mayor, a>b

L=4a 0 π/2 1 k 2 sin 2 φ dφ, k 2 = a 2 b 2 a 2 = ε 2 ,a>b

donde ε es la excentricidad de la elipse

La integral elíptica completa de segunda especie es

E(k)= 0 π/2 1 k 2 sin 2 φ ·dφ

Una expresión aproximada debida a Ramanujan es

Lπ( 3(a+b) (3a+b)(a+3b) )

Calculamos la longitud del arco de una elipse de semiejes a=2 y b=1.

>> a=2; b=1;
>> e=sqrt(a^2-b^2)/a;
>> [K,E]=ellipke(e^2);
>> 4*a*E
ans =    9.6884
>> pi*(3*(a+b)-sqrt((3*a+b)*(a+3*b)))
ans =    9.6884

Una generalización de esta integral es

E(θ,k)= 0 θ 1 k 2 sin 2 φ ·dφ

se denomina integral elíptica incompleta de segunda especie.

En la páginas tituladas 'Caída de una varilla inclinada', 'Caída del extremo libre de una cadena doblada' y 'Oscilaciones del líquido contenido en dos vasos comunicantes' se utilizan integrales elípticas completas e incompletas

Periodicidad

Las funciones elípticas E(θ,k) y F(θ,k) tiene una periodicidad que se expresa de la siguiente forma

E(θ+nπ,k)=E(θ,k)+2n·E(k) F(θ+nπ,k)=F(θ,k)+2n·F(k)

Lo que nos permite calcular los valores de las funciones elípticas para cualquier ángulo conocidos los valores para dichas funciones en el intervalo (0, π/2) y la integral elíptica completa.

Por ejemplo, para la función E(θ,k), obtenemos el valor de E(7π/6,k), conocido el valor de E(π/6,k) y la integral elíptica completa E(k)

k2=3/4;
[~,E]=ellipke(k2);
th=7*pi/6;
n=floor(th/pi);
resto=rem(th,pi);
disp(ellipticE(th,k2))
if resto<pi/2
    disp(ellipticE(resto,k2)+2*n*E)
 else
     disp(ellipticE(resto-pi,k2)+2*(n+1)*E)
 end
  
fplot(@(x) ellipticE(x,k2), [0,2*pi])
line([th,th],[0,5],'lineStyle','--')
line([resto,resto],[0,5],'lineStyle','--')
grid on
xlabel('\theta')
ylabel('E(\theta,k)')
title('Función elíptica E')

    2.9282
    2.9282

Por ejemplo, para la función F(θ,k), obtenemos el valor de F(5π/6,k), conocido el valor de F(-π/6,k)=-F(π/6,k) y la integral elíptica completa F(k)

k2=3/4;
[K,~]=ellipke(k2);
th=5*pi/6;
n=floor(th/pi);
resto=rem(th,pi);
disp(ellipticF(th,k2))
if resto<pi/2
    disp(ellipticF(resto,k2)+2*n*K)
else
    disp(ellipticF(resto-pi,k2)+2*(n+1)*K)
end

fplot(@(x) ellipticF(x,k2), [-pi/2,2*pi])
line([th,th],[-2,5],'lineStyle','--')
line([resto-pi,resto-pi],[-2,5],'lineStyle','--')
grid on
xlabel('\theta')
ylabel('F(\theta,k)')
title('Función elíptica F')

    3.7708
    3.7708

Esta propiedad no tiene utilidad si se trabaja con MATLAB. Para algunas aplicaciones, la podría tener si se utiliza el código de funciones elípticas que han limitado el ángulo al intervalo, 0≤θ≤π/2

Funciones elípticas de Jacobi

Las funciones elípticas de Jacobi sn(x|k), cn(x|k), dn(x|k) se definen del siguiente modo:

x= 0 θ dφ 1 k 2 sin 2 φ 0θ π 2 snx=sinθ cnx=cosθ dnx= 1 k 2 sin 2 θ

Obtenemos el valor de las funciones sn(x|k), cn(x|k) y dn(x|k), mediante la función ellipj de MATLAB, k2≤1.

Valores límite k2

Derivadas

d dx sn(x|k)=cn(x|k)dn(x|k) d dx cn(x|k)=sn(x|k)dn(x|k) d dx dn(x|k)= k 2 ·sn(x|k)cn(x|k)

Movimiento en un pozo doble de potencial

Este ejemplo, ilustra las funciones elípticas de Jacobi. Estudiaremos el movimiento de una partícula de masa m=1 en el potencial dado por la función Ep(x)=1-x2/2+x4/16

La energía potencial presenta dos mínimos x=±2 y un máximo local, x=0. Véase la figura más abajo

Calculamos la derivada primera y la igualamos a cero. Calculamos la derivada segunda

d E p (x) dx =x+ x 3 4 =0{ x=0 x=±2 d 2 E p (x) d x 2 =1+ 3 4 x 2

Para x=0, la deriva segunda es negativa, (máximo). Para x=±2, la derivada segunda es positiva (mínimo)

Dado el valor de la energía total E, determinamos el intervalo o intervalos en los que se puede mover la partícula.

1 2 m ( dx dt ) 2 + E p (x)=E

Como la energía cinética es siempre positiva, la partícula se puede mover en aquellas posiciones en las que E≥Ep(x). Resolvemos la ecuación bicuadrada

x 4 16 x 2 2 +1=E x 2 =4( 1± E )

Representamos la energía potencial y los intervalos para dos valores de la energía total E=1.25, 0.5

f=@(x) 1-x.^2/2+x.^4/16; %energía potencial
fplot(f,[-3.1,3.1])
E=1.25;
x1=-2*sqrt(1+sqrt(E));
x2=-x1;
line([x1,x2],[f(x1),f(x2)],'color','b')
line([x1,x1],[0,f(x1)],'lineStyle','--','color','k')
line([x2,x2],[0,f(x2)],'lineStyle','--','color','k')

E=0.5;
x1=-2*sqrt(1+sqrt(E));
x2=-2*sqrt(1-sqrt(E));
line([x1,x2],[f(x1),f(x2)],'color','r')
line([x1,x1],[0,f(x1)],'lineStyle','--','color','k')
line([x2,x2],[0,f(x2)],'lineStyle','--')
x1=2*sqrt(1-sqrt(E));
x2=2*sqrt(1+sqrt(E));
line([x1,x2],[f(x1),f(x2)],'color','r')
line([x1,x1],[0,f(x1)],'lineStyle','--','color','k')
line([x2,x2],[0,f(x2)],'lineStyle','--','color','k')

grid on
xlabel('x')
ylabel('E_p(x)')
title('Energía potencial')

Como la energía total E es positiva, para simplificar la notación escribimos E=e2

En la ecuación de la energía despejamos dt e integramos

t=± x 0 x dx 2 m ( e 2 E p (x) )

Supongamos que el móvil de masa m=1 parte de la posición del extremo derecho del intervalo, x 0 =2 1+e , con velocidad inicial negativa (dx/dt)0<0

t= 2 1+e x dz 2( e 2 1)+ z 2 z 4 8 = 8 2 1+e x dz ( 4(e+1) z 2 )( z 2 +4(e1) )

Hacemos el cambio de variable,

z=2 1+e ·cosφ dz=2 1+e ·sinφ·dφ

El resultado es

t= 2 0 Φ(x) dφ ( (1+e) cos 2 φ+(e1) ) = 2 0 Φ(x) dφ ( 2e(1+e) sin 2 φ ) = 1 e 0 Φ(x) dφ ( 1 1+e 2e sin 2 φ )

Representamos la trayectoria de un móvil en el espacio de las fases (x, v) para las energías:

e=sqrt(1.25);
m2=(1+e)/(2*e);
K=ellipke(m2);
P=4*K/sqrt(e); %periodo
t=linspace(0,P,200);
[sn,cn,dn]=ellipj(sqrt(e)*t,m2);
x=2*sqrt(1+e)*cn;
v=-2*sqrt(1+e)*sn.*dn;
hold on
plot(x,v, 'r')

%e=1;
t=linspace(0,10,200);
x=2*sqrt(2)./cosh(t);
v=-2*sqrt(2)*sinh(t)./cosh(t).^2;
plot(x,v,'k')
plot(-x,v,'k','lineStyle','--')
plot(x,-v,'k','lineStyle','--')
plot(-x,-v,'k','lineStyle','--')

e=sqrt(0.5);
m2=(1+e)/(2*e);
K=ellipke(1/m2);
P=2*K/sqrt((1+e)/2); %periodo
t=linspace(0,P,200);
[sn,cn,dn]=ellipj(sqrt((1+e)/2)*t,1/m2);
x=2*sqrt(1+e)*dn;
v=-2*sqrt(2)*e*sn.*cn;
plot(x,v,'b')
plot(-x,v,'b','lineStyle','--')

hold off
xlabel('x')
ylabel('v')
title('Espacio de las fases')
grid on

Vuelta al mundo en un día

Supongamos que partimos del polo Norte y nos movemos sobre la superficie de la Tierra con velocidad constante v.

v=ds/dt, donde ds es la longitud del arco entre dos puntos P y Q sobre la superficie de la Tierra que recorremos en un tiempo dt.

Por otra parte, nos desplazamos con la misma velocidad angular ω de la Tierra, por ejemplo, siempre estamos bajo el Sol, al mediodía. El ángulo (longitud geográfica) φ=ωt. En un día, cruzamos todos los meridianos, dando una vuelta completa a la Tierra.

La trayectoria sobre la superficie de la Tierra que describimos, se denomina espiral de Seiffert, cuyas ecuación, eliminando el tiempo t, depende de un parámetro k.

φ=ksk= ω v

Establecemos coordendas cilíndricas para describir el movimiento. Las coordendas del punto P son (φ, ρ, z), tal como vemos en la figura.

Suponiendo que la Tierra es una esfera de radio unidad R=1, se cumple que ρ2+z2=1

En la figura, se ha dibujado en azul, el meridiano de referencia, el ecuador y el eje de rotación. En rojo, el meridiano que pasa por el punto P, cuyas coordendas cilíndricas (φ, ρ, z). El código empleado para dibujar parte de esta figura, es el siguiente:

hold on

%esfera
R=1; %radio de la Tierra
theta=linspace(0,pi,50);
phi=linspace(0,2*pi,50);
[phi,theta]=meshgrid(phi,theta);
x=R*sin(theta).*cos(phi);
y=R*sin(theta).*sin(phi);
z=R*cos(theta);
h1=surfl(x,y,z);
set(h1,'EdgeColor',[0.6,0.6,0.6],'EdgeAlpha',0.9,'FaceAlpha',0.9)
shading interp
colormap(gray);

%meridiano de referencia
phi=0;
theta=linspace(0,2*pi,100);
x=R*sin(theta).*cos(phi);
y=R*sin(theta).*sin(phi);
z=R*cos(theta);
%trayectoria
h1=line(x,y,z);
set(h1,'Color',[0,0,0.5],'LineWidth',1)

%meridiano de que pasa por el punto P
phi=pi/3;
theta=linspace(0,2*pi,100);
x=R*sin(theta).*cos(phi);
y=R*sin(theta).*sin(phi);
z=R*cos(theta);
%trayectoria
h1=line(x,y,z);
set(h1,'Color',[.5,0,0],'LineWidth',1.5)

%ecuador
phi=linspace(0,2*pi,100);
theta=pi/2;
x=R*sin(theta).*cos(phi);
y=R*sin(theta).*sin(phi);
z=R*cos(theta)*ones(length(x),1);
%trayectoria
h1=line(x,y,z);
set(h1,'Color',[0,0,0.5],'LineWidth',1)

%eje de la Tierra
h1=line([0,0],[0,0],[-1.1,1.1]);
set(h1,'Color',[0,0,0.5],'LineWidth',1)

%punto P
theta=pi/6;
phi=pi/3;
x=R*sin(theta).*cos(phi);
y=R*sin(theta).*sin(phi);
z=R*cos(theta);
%punto P y centro de la esfera
plot3(x,y,z,'ro','linewidth',1,'markersize',4,'markeredgecolor','r',
'markerfacecolor','r')
plot3(0,0,0,'bo','linewidth',1,'markersize',4,'markeredgecolor','b',
'markerfacecolor','b')

view (112,28)
hold off
axis equal
axis off

Supongamos que nos movemos sobre la superficie de la Tierra, desde el punto P (φ, ρ, z) al punto muy próximo Q (φ+dφ, ρ+dρ, z+dz), la longitud del arco ds es

d s 2 = ρ 2 d φ 2 +d ρ 2 +d z 2

Trazado de la espiral

Para trazar la espiral para un determinado valor de k2, se procede del siguiente modo:

Por ejemplo, para k2=1

%esfera
R=1; %radio de la Tierra
phi=linspace(0,pi,30);
theta=linspace(0,2*pi,40);
[phi,theta]=meshgrid(phi,theta);
x=R*sin(phi).*cos(theta);
y=R*sin(phi).*sin(theta);
z=R*cos(phi);
h1=surfl(x,y,z);
set(h1,'EdgeColor',[0.6,0.6,0.6],'EdgeAlpha',0.7,'FaceAlpha',0.7)
shading interp
colormap(gray);
axis equal

k2=1;  %cambiar este valor
s=linspace(0,3*pi,200);
[sn,cn,dn]=ellipj(s,k2);
phi=sqrt(k2)*s;
xp=sn.*cos(phi);
yp=sn.*sin(phi);
zp=cn;
%trayectoria
h1=line(xp,yp,zp);
set(h1,'Color',[.7,0,0],'LineWidth',1.5)
view (130,20)
hold off
axis equal
grid on
xlabel('x')
ylabel('y')
zlabel('z')
title('Espiral de Seiffert')
axis off

Para k2=1, las funciones sn(x|1) y cn(x|1) dejan de ser periódicas, tendiendo hacia el valor 1 y 0, respectivamente, por lo que la espiral parte del polo Norte y tiende asintóticamente, hacia el ecuador

Cuando k2<1, las funciones sn(x|k) y cn(x|k) son periódicas y hemos determinado su periodo 4K(k), siendo K (ellipke en MATLAB), la integral elíptica completa de primera especie

La espiral será una trayectoria cerrada si el punto inicial y final difieren en φ=2nπ, con n=1,2,3... e igual a m periodos de sn o cn para un valor dado de k2

2nπ k =m·4K( k )m,n=1,2,3...

Dados los números enteros m y n, resolvemos la ecuación transcendente para calcular el valor de k2 que hace que la espiral sea una trayectoria cerrada. Sea por ejemplo m=1 y n=1

n=1;
m=1;
f=@(x) 2*n*pi/sqrt(x)-4*m*ellipke(x);
k2=fzero(f, 0.5);
disp(k2)
   0.6284

En el script de la trayectoria cambiamos el valor de la variable k2 a 0.6284, obtiendo la trayectoria cerrada

Referencias

Alain J. Brizard. A primer on elliptic functions with applications in classical mechanics. November 26, 2007. https://arxiv.org/abs/0711.4064.

Paul Erdos. Spiraling the Earth with C. G. J. Jacobi. Am. J. Phys. 68 (10) October 2000, pp. 888-895

Ejemplos en el curso de Física

Funciones elípticas: ellipke, ellipticF, ellipticE, ellipticPi, ellipj

Movimiento sobre una cúpula semiesférica

Movimiento sobre una superficie semicircular cóncava

Modelo PREM del interior de la Tierra

Caída de una varilla inclinada

Flexión de una viga en voladizo (II)

Pandeo de una barra delgada empotrada en un extremo

Flexión de una regla

El péndulo simple (II)

El péndulo interrumpido

Oscilador no lineal

Oscilaciones de una esfera que flota en el agua

Campo eléctrico de un sistema de dos o más cargas eléctricas

Campo y potencial eléctrico de una distribución continua de carga

Campo magnético producido por una corriente circular en un punto fuera de su eje

Campo magnético producido por una espira de forma cualesquiera

Circuitos acoplados. Coeficiente de inducción mutua

Una partícula desliza sobre una pista circular que gira