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\underset{0}{\overset{\pi /4}{\int }}ds$

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

$\begin{array}{l}\left\{\begin{array}{l}x=r\mathrm{cos}\theta \\ y=r\mathrm{sin}\theta \end{array}\\ \left\{\begin{array}{l}dx=\mathrm{cos}\theta ·dr-r\mathrm{sin}\theta ·d\theta \\ dy=\mathrm{sin}\theta ·dr+r\mathrm{cos}\theta ·d\theta \end{array}\\ d{s}^{2}=d{x}^{2}+d{y}^{2}=d{r}^{2}+{r}^{2}d{\theta }^{2}\\ {r}^{2}={a}^{2}\mathrm{cos}\left(2\theta \right)\\ rdr=-{a}^{2}\mathrm{sin}\left(2\theta \right)·d\theta \\ d{s}^{2}={a}^{2}\frac{d{\theta }^{2}}{\mathrm{cos}\left(2\theta \right)}\end{array}$

Integramos, respecto del ángulo, θ

$l=4\underset{0}{\overset{\pi /4}{\int }}ds=4a\underset{0}{\overset{\pi /4}{\int }}\frac{d\theta }{\sqrt{\mathrm{cos}\left(2\theta \right)}}$

Hacemos el cambio de variable

$\left\{\begin{array}{l}\mathrm{cos}\left(2\theta \right)={\mathrm{cos}}^{2}\phi \\ -2\mathrm{sin}\left(2\theta \right)·d\theta =-2\mathrm{cos}\phi \mathrm{sin}\phi ·d\phi \end{array}$

La integral se convierte en

$\begin{array}{l}l=4a\underset{0}{\overset{\pi /2}{\int }}\frac{\mathrm{sin}\phi d\phi }{\sqrt{1-{\mathrm{cos}}^{4}\phi }}=4a\underset{0}{\overset{\pi /2}{\int }}\frac{d\phi }{\sqrt{1+{\mathrm{cos}}^{2}\phi }}=4a\underset{0}{\overset{\pi /2}{\int }}\frac{d\phi }{\sqrt{2-{\mathrm{sin}}^{2}\phi }}=\\ \frac{4a}{\sqrt{2}}\underset{0}{\overset{\pi /2}{\int }}\frac{d\phi }{\sqrt{1-\frac{1}{\sqrt{2}}{\mathrm{sin}}^{2}\phi }}=\frac{4a}{\sqrt{2}}K\left(\frac{1}{\sqrt{2}}\right)\end{array}$

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

La integral elíptica completa de primera especie es

$K\left(k\right)=\underset{0}{\overset{\pi /2}{\int }}\frac{d\phi }{\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }}$

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

Una generalización de esta integral

$F\left(\theta ,k\right)=\underset{0}{\overset{\theta }{\int }}\frac{d\phi }{\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }}\text{ }0\le \theta \le \frac{\pi }{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

$\frac{{x}^{2}}{{a}^{2}}+\frac{{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

$\begin{array}{l}ds=\sqrt{d{x}^{2}+d{y}^{2}}=\sqrt{{a}^{2}{\mathrm{sin}}^{2}\phi +{b}^{2}{\mathrm{cos}}^{2}\phi }·d\phi =b\sqrt{1-\frac{{b}^{2}-{a}^{2}}{{b}^{2}}{\mathrm{sin}}^{2}\phi }·d\phi \\ =b\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }·d\phi \text{ }\left(b>a\right)\end{array}$

La longitud de la elipse es

$L=b\underset{0}{\overset{2\pi }{\int }}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }·d\phi =4b\underset{0}{\overset{\pi /2}{\int }}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }·d\phi$

Habitualmente, el semieje mayor, a>b

$L=4a\underset{0}{\overset{\pi /2}{\int }}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }d\phi ,\text{ }{k}^{2}=\frac{{a}^{2}-{b}^{2}}{{a}^{2}}={\epsilon }^{2},\text{ }a>b$

donde ε es la excentricidad de la elipse

La integral elíptica completa de segunda especie es

$E\left(k\right)=\underset{0}{\overset{\pi /2}{\int }}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }·d\phi$

Una expresión aproximada debida a Ramanujan es

$L\approx \pi \left(3\left(a+b\right)-\sqrt{\left(3a+b\right)\left(a+3b\right)}\right)$

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\left(\theta ,k\right)=\underset{0}{\overset{\theta }{\int }}\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }·d\phi$

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

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

$\begin{array}{l}E\left(\theta +n\pi ,k\right)=E\left(\theta ,k\right)+2n·E\left(k\right)\\ F\left(\theta +n\pi ,k\right)=F\left(\theta ,k\right)+2n·F\left(k\right)\end{array}$

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:

$\begin{array}{l}x=\underset{0}{\overset{\theta }{\int }}\frac{d\phi }{\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\phi }}\text{ }0\le \theta \le \frac{\pi }{2}\\ \text{sn}\text{\hspace{0.17em}}x=\mathrm{sin}\theta \\ \text{cn}\text{\hspace{0.17em}}x=\mathrm{cos}\theta \\ \text{dn}\text{\hspace{0.17em}}x=\sqrt{1-{k}^{2}{\mathrm{sin}}^{2}\theta }\end{array}$

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.

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

$\begin{array}{l}\text{sn}\left(x|k\right)=\frac{1}{\sqrt{k}}\text{sn}\left(\sqrt{k}·x|\frac{1}{k}\right)\\ \text{cn}\left(x|k\right)=\text{dn}\left(\sqrt{k}·x|\frac{1}{k}\right)\\ \text{dn}\left(x|k\right)=c\text{n}\left(\sqrt{k}·x|\frac{1}{k}\right)\end{array}$

$\begin{array}{l}\frac{d}{dx}\text{sn}\left(x|k\right)=\text{cn}\left(x|k\right)\text{dn}\left(x|k\right)\hfill \\ \begin{array}{l}\frac{d}{dx}\text{cn}\left(x|k\right)=-\text{sn}\left(x|k\right)\text{dn}\left(x|k\right)\\ \frac{d}{dx}\text{dn}\left(x|k\right)=-{k}^{2}·\text{sn}\left(x|k\right)\text{cn}\left(x|k\right)\end{array}\hfill \end{array}$

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

$\begin{array}{l}\frac{d{E}_{p}\left(x\right)}{dx}=-x+\frac{{x}^{3}}{4}=0\left\{\begin{array}{l}x=0\\ x=±2\end{array}\\ \frac{{d}^{2}{E}_{p}\left(x\right)}{d{x}^{2}}=-1+\frac{3}{4}{x}^{2}\end{array}$

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.

$\frac{1}{2}m{\left(\frac{dx}{dt}\right)}^{2}+{E}_{p}\left(x\right)=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

$\begin{array}{l}\frac{{x}^{4}}{16}-\frac{{x}^{2}}{2}+1=E\\ {x}^{2}=4\left(1±\sqrt{E}\right)\end{array}$

• Para E>1, se descarta el signo menos, la partícula se mueve en el intervalo

• $\left[-2\sqrt{1+\sqrt{E}},\text{\hspace{0.17em}}2\sqrt{1+\sqrt{E}}\right]$

• Si E<1, la partícula se puede mover en uno o en el otro intervalo

• $\left[-2\sqrt{1+\sqrt{E}},\text{\hspace{0.17em}}-2\sqrt{1-\sqrt{E}}\right],\text{\hspace{0.17em}}\left[2\sqrt{1-\sqrt{E}},\text{\hspace{0.17em}}2\sqrt{1+\sqrt{E}}\right]$

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=±\underset{{x}_{0}}{\overset{x}{\int }}\frac{dx}{\sqrt{\frac{2}{m}\left({e}^{2}-{E}_{p}\left(x\right)\right)}}$

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=-\underset{2\sqrt{1+e}}{\overset{x}{\int }}\frac{dz}{\sqrt{2\left({e}^{2}-1\right)+{z}^{2}-\frac{{z}^{4}}{8}}}=-\sqrt{8}\underset{2\sqrt{1+e}}{\overset{x}{\int }}\frac{dz}{\sqrt{\left(4\left(e+1\right)-{z}^{2}\right)\left({z}^{2}+4\left(e-1\right)\right)}}$

Hacemos el cambio de variable,

$\begin{array}{l}z=2\sqrt{1+e}·\mathrm{cos}\phi \\ dz=-2\sqrt{1+e}·\mathrm{sin}\phi ·d\phi \end{array}$

$t=\sqrt{2}\underset{0}{\overset{\Phi \left(x\right)}{\int }}\frac{d\phi }{\sqrt{\left(\left(1+e\right){\mathrm{cos}}^{2}\phi +\left(e-1\right)\right)}}=\sqrt{2}\underset{0}{\overset{\Phi \left(x\right)}{\int }}\frac{d\phi }{\sqrt{\left(2e-\left(1+e\right){\mathrm{sin}}^{2}\phi \right)}}=\frac{1}{\sqrt{e}}\underset{0}{\overset{\Phi \left(x\right)}{\int }}\frac{d\phi }{\sqrt{\left(1-\frac{1+e}{2e}{\mathrm{sin}}^{2}\phi \right)}}$

• Para e>1, k2=(1+e)/(2e)≤1

• $\sqrt{e}·t=\underset{0}{\overset{\Phi \left(x\right)}{\int }}\frac{d\phi }{\sqrt{\left(1-{k}^{2}{\mathrm{sin}}^{2}\phi \right)}}\left\{\begin{array}{l}{k}^{2}=\frac{1+e}{2e}\\ \Phi \left(x\right)=\mathrm{arccos}\left(\frac{x}{2\sqrt{1+e}}\right)\end{array}$

Despejamos x en función del tiempo t, empleando la función elíptica sn de Jacobi

$\begin{array}{l}\text{sn}\left(\sqrt{e}·t|k\right)=\mathrm{sin}\left(\Phi \left(x\right)\right)\\ \text{sn}\left(\sqrt{e}·t|k\right)=\sqrt{1-{\mathrm{cos}}^{2}\left(\Phi \left(x\right)\right)}\\ \text{sn}\left(\sqrt{e}·t|k\right)=\sqrt{1-{\left(\frac{x}{2\sqrt{1+e}}\right)}^{2}}\\ {x}^{2}=4\left(1+e\right)·{\text{cn}}^{2}\left(\sqrt{e}·t|k\right)\\ x=2\sqrt{1+e}·\text{cn}\left(\sqrt{e}·t|k\right)\end{array}$

Se ha utilizado la relación, sn2(x|k)+cn2(x|k)=1

Derivando respecto del tiempo, obtenemos la velocidad del móvil

$\frac{dx}{dt}=-2\sqrt{e\left(1+e\right)}·\text{sn}\left(\sqrt{e}·t|k\right)\text{dn}\left(\sqrt{e}·t|k\right)$

• Para e=1, k2=(1+e)/(2e)=1

• Para este valor límite de k2, la posición y velocidad toman las expresiones

$\begin{array}{l}x=2\sqrt{2}·\frac{1}{\mathrm{cosh}t}\\ \frac{dx}{dt}=-2\sqrt{2}·\mathrm{tanh}t·\frac{1}{\mathrm{cosh}t}=-2\sqrt{2}\frac{\mathrm{sinh}t}{{\mathrm{cosh}}^{2}t}\end{array}$

• Para e<1, k2=(1+e)/(2e)>1

• Para estos valores de k2 fuera del intervalo [0,1], la posición y velocidad toman las expresiones

$\begin{array}{l}x=2\sqrt{1+e}·\text{dn}\left(\sqrt{\frac{1+e}{2e}}\sqrt{e}·t|\frac{1}{k}\right)=2\sqrt{1+e}\text{·dn}\left(\sqrt{\frac{1+e}{2}}·t|\frac{1}{k}\right)\\ \frac{dx}{dt}=-2\sqrt{e\left(1+e\right)}·\frac{1}{\sqrt{\frac{1+e}{2e}}}\text{sn}\left(\sqrt{\frac{1+e}{2}}·t|\frac{1}{k}\right)\text{cn}\left(\sqrt{\frac{1+e}{2}}·t|\frac{1}{k}\right)\\ =-2\sqrt{2e}\text{·sn}\left(\sqrt{\frac{1+e}{2}}·t|\frac{1}{k}\right)\text{cn}\left(\sqrt{\frac{1+e}{2}}·t|\frac{1}{k}\right)\end{array}$

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 $[ −2 1+ E , 2 1+ E ]$, la trayectoria en el espacio de las fases es cerrada y el tiempo que tarda en completarla es

• $P=\frac{4K\left(k\right)}{\sqrt{e}}$

• E=1. Como el caso anterior, el móvil parte de la posición inicial $x 0 =2 1+e$, tardando un tiempo infinito en alcanzar el origen x=0

• E=0.5. La partícula se mueve en el intervalo $[ 2 1− E , 2 1+ E ]$ 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

• $P=\frac{2K\left(k\right)}{\sqrt{\left(1+e\right)/2}}$

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.

$\phi =ks\text{ }k=\frac{\omega }{v}$

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

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}={\rho }^{2}d{\phi }^{2}+d{\rho }^{2}+d{z}^{2}$

• Teniendo en cuenta que φ=ks y que ρ2+z2=1, escribimos ds en función de ρ y

• $\begin{array}{l}d{s}^{2}={\rho }^{2}{k}^{2}d{s}^{2}+\frac{1}{{\rho }^{2}-1}d{\rho }^{2}\\ \left(1-{\rho }^{2}{k}^{2}\right)d{s}^{2}=\frac{1}{{\rho }^{2}-1}d{\rho }^{2}\\ ds=\frac{d\rho }{\sqrt{\left(1-{\rho }^{2}{k}^{2}\right)\left({\rho }^{2}-1\right)}}\end{array}$

La longitud del arco s que recorre desde el polo Norte hasta un punto P que dista ρ del eje de rotación es

$s\left(\rho ,{k}^{2}\right)=\underset{0}{\overset{\rho }{\int }}\frac{d\rho }{\sqrt{\left(1-{\rho }^{2}{k}^{2}\right)\left({\rho }^{2}-1\right)}}$

La inversa de esta función es ρ= sn(s|k)

• Teniendo en cuenta que φ=ks y que ρ2+z2=1, escribimos ds en función de z y dz

• $\begin{array}{l}d{s}^{2}={\rho }^{2}d{\phi }^{2}+d{\rho }^{2}+d{z}^{2}\\ d{s}^{2}={\rho }^{2}{k}^{2}d{s}^{2}+\frac{1}{{z}^{2}-1}d{z}^{2}\\ \left(1-{k}^{2}+{z}^{2}{k}^{2}\right)d{s}^{2}=\frac{1}{{z}^{2}-1}d{z}^{2}\\ ds=\frac{dz}{\sqrt{\left(1-{k}^{2}+{z}^{2}{k}^{2}\right)\left({z}^{2}-1\right)}}\end{array}$

La longitud del arco s que recorre desde el polo Norte hasta un punto P que dista z del plano del ecuador es

$s\left(z,{k}^{2}\right)=\underset{0}{\overset{z}{\int }}\frac{dz}{\sqrt{\left(1-{k}^{2}+{z}^{2}{k}^{2}\right)\left({z}^{2}-1\right)}}$

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

$\frac{2n\pi }{k}=m·4K\left(k\right)\text{ }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

Flexión de una viga en voladizo (II)

Flexión de una regla

El péndulo simple (II)

El péndulo interrumpido

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

Carga inducida en un conductor esférico