Implementación de un péndulo esférico

El problema parece ser que no está teniendo en cuenta el tamaño del paso a tiempo al integrar. Esto debería ser obvio cuando estás haciendo cosas como

phi_v += phi_f

en el código. ¹ La aceleración no se puede simplemente sumar a la velocidad (¡las unidades no coinciden!). La relación esperada es,

$$a_\phi =\frac{\mathrm d v_\phi}{\mathrm dt}\approx\frac{\Delta v_\phi}{\Delta t}\implies v_\phi\approx v_\phi'+ a_\phi \Delta t$$ dónde$v_\phi'$es el valor anterior es la velocidad.

Lo que deberías estar haciendo es usar las 4 ecuaciones (2 posiciones, 2 velocidades),

$$\mathbf{v}=\frac{\mathrm d\mathbf x}{\mathrm dt} \text{ and } \mathbf{a}=\frac{1}{m}\frac{\mathrm d\mathbf f}{\mathrm dt} ,$$ dónde$\mathbf v=[\dot{\phi},\,\dot{\theta}]^T$y de manera similar para$\mathbf x$, y resuélvalos a través de un integrador de salto de rana, como la velocidad Verlet . Sin embargo, dado que tiene una aceleración dependiente de la velocidad (es decir,$a= f(\mathbf x,\,\mathbf v)$), necesita usar una versión modificada del integrador (del que hablo en esta otra respuesta mía ),

$$\begin{align} a_1(t)&=\text{compute from $x(t)$ and $v(t)$}\tag{1}\\ x(t+\Delta t)&=x(t)+v(t)\Delta t+\frac{1}{2}a_1(t)\Delta t^2\tag{2}\\ v(t+\Delta t)&=v(t)+a_1(t)\Delta t\tag{3}\\ a_2(t+\Delta t)&=\text{compute from $x(t+\Delta t)$, $v(t+\Delta t)$}\tag{4}\\ v(t+\Delta t)&=v(t+\Delta t)+\frac{1}{2}\left(a_2(t+\Delta t)-a_1(t)\right)\Delta t\tag{5} \end{align}$$ que luego simplemente recorre hasta que esté satisfecho (por ejemplo,$t\geq t_\text{end}$). Proporcionó$\Delta t$es lo suficientemente pequeño, esto debería proporcionar un solucionador más estable que la integración de Verlet de velocidad sola.

^{1. OP agregó el timestepstérmino a las ecuaciones en la versión 6 ; esta respuesta se publicó antes de esa edición.}

Me parece que estás tratando de resolver un sistema de ecuaciones diferenciales ordinarias estimando la solución en pequeños pasos. ¿Estoy en lo correcto? Si es así, hay una serie de cosas a considerar.

Primero, debe tener una buena estimación del operador derivado para aproximar el estado futuro dadas las condiciones iniciales. Este tipo de problema se llama problema de valor inicial (ivp).

Dado que la ecuación involucra derivadas temporales de primer y segundo orden, un paso común es definir las variables de impulso y propagar un sistema más grande. Esto ayudará a evitar errores a medida que aumenta el tiempo total. He visto ese problema en los códigos de trazado de rayos y lo que sugiero es tan común que la mayoría de las personas nunca propagan una ecuación de segundo orden. La idea es:

p_theta = d(theta)/dt

p_phi = d(phi)/dt

entonces tu ecuacion es

d(p_theta)/dt = (p_phi)^2*cos(theta)*sin(theta) - g/l*sin(theta)

d(p_phi)/dt = -2*p_theta*p_phi*cot(theta)

más las dos ecuaciones que definen las p.

Un simple paso de Euler se implementaría como,

p_theta(t0+dt) = p_theta(t0) + (d(p_theta)/dt)(t0)*dt
p_phi(t0+dt) = p_phi(t0) + (d(p_phi)/dt)(t0)*dt

más las ecuaciones para theta y phi en las definiciones de p.

Esta es la configuración típica. Ahora, el paso de Euler es muy pobre y nunca recomendado. Sería mejor implementar un método de orden superior como RK4 o RK5 (4), etc., con control de tamaño de paso.

Aparte de eso, el uso de ángulos a veces es un problema, ya que lleva a la posibilidad de que en algún paso no pueda determinar de manera única el siguiente valor debido a que el sistema es singular. En las simulaciones aeroespaciales usan cuarterones para arreglar esto. Puede consultar Goldstein Classical Mechanics para obtener detalles sobre las matemáticas o Zipfel Modeling and Simulation of Aerospace Vehicle Dynamics. Creo que no necesitas esta maquinaria en este momento. Debe escribir las ecuaciones correctamente y probar un algoritmo de paso simple antes de volverse demasiado sofisticado. Creo que Python tiene un paquete de resolución de ODE, por lo que no necesita escribir el suyo propio, solo configure correctamente y llame. Espero que eso ayude.