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:
Calculamos la longitud de arco ds de una porción infinitesimal de la lemniscata de Bernoulli
Integramos, respecto del ángulo, θ
Hacemos el cambio de variable
La integral se convierte en
>> 4*ellipke(1/2)/sqrt(2) ans = 5.2441
La integral elíptica completa de primera especie es
Una integral similar, la encontramos al calcular el periodo de un péndulo para cualquier amplitud θ0<π
Una generalización de esta integral
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
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
La longitud de la elipse es
Habitualmente, el semieje mayor, a>b
donde ε es la excentricidad de la elipse
La integral elíptica completa de segunda especie es
Una expresión aproximada debida a Ramanujan es
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
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
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:
Obtenemos el valor de las funciones sn(x|k), cn(x|k) y dn(x|k), mediante la función
Representamos la función elíptica de Jacobi sn(x|k) para dos valores de k2<1: 0.1 y 0.9
x=linspace(0,10,200); [sn,~,~]=ellipj(x,0.1); hold on plot(x,sn) [sn,~,~]=ellipj(x,0.9); plot(x,sn) hold off ylim([-1.05,1.05]) xlabel('x') ylabel('sn(x)') legend ('0.1','0.9', 'location', 'southwest') title('Funciones elípticas de Jacobi, sn') grid on
Representamos la función elíptica de Jacobi cn(x|k) para dos valores de k2<1: 0.1 y 0.9
x=linspace(0,10,200); [~,cn,~]=ellipj(x,0.1); hold on plot(x,cn) [~,cn,~]=ellipj(x,0.9); plot(x,cn) hold off ylim([-1.05,1.05]) xlabel('x') ylabel('cn(x)') legend ('0.1','0.9', 'location', 'southwest') title('Funciones elípticas de Jacobi, cn') grid on
Cuando x=K(k), integral elíptica completa de primera especie, sn(x|k) vale 1. Cuando x=2K(k), sn(x|k)=0. Por lo que 2K(k) es medio periodo, tal como podemos comprobar en las gráficas de las funciones sn(x|k) y cn(x|k)
>> K=ellipke(0.1); >> [sn,cn,~]=ellipj(2*K,0.1); >> sn,cn sn = 1.2246e-16 cn = -1 >> 2*K ans = 3.2249 >> K=ellipke(0.9); >> [sn,cn,dn]=ellipj(2*K,0.9); >> sn,cn sn = 1.2246e-16 cn = -1 >> 2*K ans = 5.1562
Las funciones sn(x|k) y cn(x|k) están relacionadas, como sus análogas trigonométricas
sn2(x|k)+cn2(x|k)=1
Representamos la función elíptica de Jacobi dn(x|k) para dos valores de k2<1: 0.1 y 0.9
x=linspace(0,10,200); [~,~,dn]=ellipj(x,0.1); hold on plot(x,dn) [~,~,dn]=ellipj(x,0.9); plot(x,dn) hold off xlabel('x') ylabel('dn(x)') legend ('0.1','0.9', 'location', 'southwest') title('Funciones elípticas de Jacobi, dn') grid on
Se trata de una función periódica de periodo 2K(k)
>> K=ellipke(0.1); >> [~,~,dn]=ellipj(2*K,0.1); >> dn dn = 1 >> K=ellipke(0.9); >> [~,~,dn]=ellipj(2*K,0.9); >> dn dn = 1
Valores límite k2
Para k2=1, las funciones sn(x|1) y cn(x|1) dejan de ser periódicas
x=linspace(0,10,200); [sn,cn,~]=ellipj(x,1); hold on plot(x,sn) plot(x,cn) hold off ylim([-.05,1.05]) xlabel('x') ylabel('cn(x)') legend ('sn','cn', 'location', 'southwest') title('Funciones elípticas de Jacobi, k^2=1') grid on
Comprobamos la equivalencia:
sn(x|1)=tanh(x),
cn(x|1)=1/cosh(x),
dn(x|1)=1/cosh(x)
Añadiendo al script, las líneas de código
... fplot(@(x) tanh(x),[0,10]) fplot(@(x) 1./cosh(x),[0,10]) ...
Para k2=0, las funciones sn(x|k) y cn(x|k) son iguales a sus análogas trigonométricas
sn(x|0)=sin(x),
cn(x|0)=cos(x),
dn(x|0)=1
Cuando k2>1
Las funciones elípticas se calculan habitualmente para 0≤k2≤1 (es el intervalo que permite MATLAB). Cuando k2>1, es posible calcularlas mediante las relaciones:
Derivadas
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
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.
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
Para E>1, se descarta el signo menos, la partícula se mueve en el intervalo
Si E<1, la partícula se puede mover en uno o en el otro intervalo
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
Supongamos que el móvil de masa m=1 parte de la posición del extremo derecho del intervalo, , con velocidad inicial negativa (dx/dt)0<0
Hacemos el cambio de variable,
El resultado es
Para e>1, k2=(1+e)/(2e)≤1
Para e=1, k2=(1+e)/(2e)=1
Para e<1, k2=(1+e)/(2e)>1
Despejamos x en función del tiempo t, empleando la función elíptica sn de Jacobi
Se ha utilizado la relación, sn2(x|k)+cn2(x|k)=1
Derivando respecto del tiempo, obtenemos la velocidad del móvil
Para este valor límite de k2, la posición y velocidad toman las expresiones
Para estos valores de k2 fuera del intervalo [0,1], la posición y velocidad toman las expresiones
Representamos la trayectoria de un móvil en el espacio de las fases (x, v) para las energías:
E=1.25. La partícula se mueve en el intervalo , la trayectoria en el espacio de las fases es cerrada y el tiempo que tarda en completarla es
E=1. Como el caso anterior, el móvil parte de la posición inicial , tardando un tiempo infinito en alcanzar el origen x=0
E=0.5. La partícula se mueve en el intervalo alrededor del mínimo x=2, la trayectoria en el espacio de las fases es cerrada, tardando un tiempo P en completarla. También es posible que la partícula se mueva en el intervalo centrado en x=-2. Pero no es posible que una partícula con energía E<1, pase de un intervalo al otro
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.
Establecemos coordendas cilíndricas para describir el movimiento. Las coordendas del punto P son (φ, ρ, z), tal como vemos en la figura.
- φ es el ángulo que forma el meriano que pasa por P con el meridiano de referencia
- z es la altura de P sobre el ecuador
- ρ es la distancia del P al eje de rotación Z
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
-
Teniendo en cuenta que φ=ks y que ρ2+z2=1, escribimos ds en función de ρ y dρ
Teniendo en cuenta que φ=ks y que ρ2+z2=1, escribimos ds en función de z y dz
La longitud del arco s que recorre desde el polo Norte hasta un punto P que dista ρ del eje de rotación es
La inversa de esta función es ρ= sn(s|k)
La longitud del arco s que recorre desde el polo Norte hasta un punto P que dista z del plano del ecuador es
La inversa de esta función es z= cn(s|k)
Para trazar la espiral para un determinado valor de k2, se procede del siguiente modo:
- Dada la longitud del arco s, se obtiene el ángulo φ=s/k
- Se determinan las coordenadas ρ, z del punto P utilizando las funciones sn y cn
ρ= sn(s|k), z= cn(s|k)
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 (
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
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
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
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