Encontrar el periodo de una EDO no lineal

Question

Estoy estudiando un sistema físico periódico con una EDO no lineal $$x''=f(x)+g(x)x'^2$$

Creo que la periodicidad viene de la $x'^2$ porque esto proporciona dos números posibles para dar un mismo valor del lado derecho.

A continuación se muestran tres curvas numéricas de esta ecuación con $f(x)=x-x^3$ y $g(x)=2/x-x$ .

Podemos ver que la curva oscila alrededor del punto fijo (haciendo $x''=0$ y $x'=0$ Aquí el punto fijo es $x^*=1$ )

Puedo resolver el período para la solución de punto casi fijo $x(t)=x^*+\epsilon \cdot \cos{\omega t}$ y esta perturbación me da una frecuencia $\omega=\sqrt{-f'(x^*)}$ . Para el ejemplo concreto que he dado, $\omega=\sqrt{-f'(1)}=\sqrt{-(1-3\cdot 1^2)}=\sqrt{2}$ así que $T=2\pi/\omega=\sqrt{2}\pi \approx 4.44$ y coincide con la curva roja bastante bien.

Mi pregunta es ¿cómo puedo resolver analíticamente el periodo de las curvas alejadas de la solución del punto fijo?

Gracias por su atención.

Answer 1

Si tienes una función de Hamilton sin términos separables como $$ H(x,p)=\frac1{2m(x)}p^2+V(x) $$ la dinámica resultante es \begin {align} \dot x &= ~~~H_p= \frac {p}{m(x)} \\ \dot p &= -H_x= \frac {m'(x)}{2m(x)^2}p^2-V'(x) \end {align} Ahora elimina $p,\dot p$ para conseguir $$ \ddot x=-\frac{m'(x)}{m(x)^2}\dot xp+\frac{\dot p}{m(x)} =-\frac{m'(x)}{m(x)}\dot x^2+\frac12\frac{m'(x)}{m(x)}\dot x^2-\frac{V'(x)}{m(x)} \\~\\ \ddot x+\frac{m'(x)}{2m(x)}\dot x^2+\frac{V'(x)}{m(x)}=0 $$ Ahora para el escalar $x$ las ecuaciones $(\ln|m(x)|)'=-2g(x)$ y $V'(x)=-m(x)f(x)$ son siempre integrables, lo que significa que su EDO siempre tiene una primera integral. Como ahora todas las soluciones tienen que permanecer en las curvas de nivel de esta primera integral, tienen que ser periódicas siempre que esa curva de nivel no contenga un punto estacionario. El periodo se puede calcular como $$ \frac T2=\int_{x_1}^{x_2}\frac{dx}{\sqrt{2(V(x_1)-V(x))/m(x)}} $$ donde $x_1<x_2$ , $V(x_2)=V(x_1)$ son los puntos extremos en $x$ dirección de una curva de nivel.

---

En su ejemplo obtengo $m(x)=e^{x^2}/x^4$ , $V'(x)=-e^{x^2}\frac{1-x^2}{x^3}=e^{x^2}(x^{-1}-x^{-3})$ . Como $$ \frac d{dx}e^{x^2}x^{-2} = e^{x^2}(2x^{-1}-2x^{-3}) $$ obtenemos $$ V(x)=\frac{e^{x^2}}{2x^2} $$ para que el punto de inflexión inferior pueda calcularse mediante la función Lambert-W a partir del superior, $$ -x_1^2e^{-x_1^2}=-x_2^2e^{-x_2^2}\implies x_1=\sqrt{-W_{0}(-x_2^2e^{-x_2^2})} $$

La integración numérica de la integral anterior para el periodo da la gráfica

from scipy.special import lambertw
from scipy.integrate import quad

E = 1.001+np.linspace(0,30,150+1); V0s = E*np.exp(1)

def integrand(V0): return lambda x: 1/(x*(V0*x**2*np.exp(-x**2)-1)**0.5)
def x1(V0): return (-lambertw(-1/V0    ).real)**0.5 
def x2(V0): return (-lambertw(-1/V0, -1).real)**0.5 
T = np.array([ 2*quad(integrand(V0), x1(V0), x2(V0))[0] for V0 in V0s])

plt.plot(E,T/(2**0.5*np.pi)); plt.grid(); 
plt.xlabel("$V_0=V(x_{1/2})$ in multiples of $e/2$"); 
plt.ylabel("$T$ in multiples of $\sqrt{2}\pi$"); plt.show()