python – Pendulo Simple: Como puedo obtener la evolución exacta en el tiempo del ángulo?


Estoy trabajando en un problema del pendulo simple para obtener la evolución del ángulo θ en un intervalo de tiempo luego de conocer un conjunto de datos. Estoy llevando mi solucion en funciones elementales cuyo resultado es:

mas

donde el valor de la frecuencia angular ω0 se obtiene de la expresión:

Frecuencia

y los de A y ϕ ajustando la solución aproximada a las condiciones iniciales.

Para llegar a la solucion, me estoy basando en el trabajo de Meléndez et alt, llegando a la siguiente expresión exacta para la evolución temporal de θ(t):

Formula

Donde en R lo tengo de esta forma:

PosicionAng <- function(t,AngIn,FrecAng)
{
  2 * acos(sin(AngIn/2) * sn(K.fun((sin(AngIn/2)^2)-FrecAng * t), (sin(AngIn/2)^2)))
}

donde K la integral elíptica completa de primera clase:
K(m)

El cálculo de los valores de la función θ(t) se me esta haciendo complicado y al correr mi programa obtengo una grafica muy distinta basado en el trabajo de Melendez.

Mi grafica es:

grafica

Quisiera saber si tengo un error en el codigo o puedo enfocar la solucion de otra manera.

Partiendo de una situacion donde el angulo es: pi/4 (45 grados) y la longuitud de la cuerda 0.8m y la gravedad en 9.8 m/s^2.

Lo he trabajado con la library(elliptic) poder encontrar tal resultado pero no se que le pasa a mi grafica.

He construido el siguiente codigo:

library(elliptic)

L <- 0.8
g <- 9.8
AngIn <- pi/4

#Frecuencia angular
FrecuenciaAng <- function(g,L)
{
  2*pi*(sqrt(g/L))
}

#Formula para obtener el valor del angulo en determinado segundo t.
PosicionAng <- function(t,AngIn,FrecAng)
{
  2 * acos(sin(AngIn/2) * sn(K.fun((sin(AngIn/2)^2)-FrecAng * t), (sin(AngIn/2)^2)))
}

#'Funcion que aplica la formula anterior utilizando los diferentes intervalos (t).
Prueba <- function(L,g,AngIn)
{
  FrecAng <- FrecuenciaAng
  t = seq(from= 0, to= 10, by= 0.1)
  PosAng <- PosicionAng(t,AngIn,FrecAng)
  return (PosAng)
}

#' Evolucion del angulo en funcion del tiempo
EvolAng <- Prueba(L,g,AngIn)

plot(EvolAng,
     type = "l",
     main = "Evolucion del angulo en el tiempo",
     ylab = "Angulo",
     xlab = "Tiempo"
)