Nonlinear pendulum in R: Problem to obtain the exact evolution in time of the angle

Question

I am working on a simple pendulum problem to obtain the evolution of the angle θ in a time interval after knowing a data set. To arrive at the solution, I am basing myself on the work of Meléndez et alt, arriving at the following exact expression for the temporal evolution of θ(t).

Where in R I have it this way:

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

The calculation of the values of the function θ(t) is getting complicated and when I run my program I get a very different graph based on Melendez's work.

My graph it has no elliptical curvature and is linear.

I would like to know if I have an error in the code or I can approach the solution in another way.

Starting from a situation where the angle is: pi/4 (45 degrees) and the length of the rope is 0.8m and the gravity is 9.8 m/s^2.

I have worked with the library (elliptic) to find such a result but I do not know what happens to my graph.

I have built the following code:

library(elliptic)

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

#Angular frequency
FrequencyAng <- function(g,L)
{
  2*pi*(sqrt(g/L))
}

#Obtain the value of the angle in a certain second t.
PositionAng <- function(t,AngIn,FreAng)
{
  2 * acos(sin(AngIn/2) * sn(K.fun((sin(AngIn/2)^2)-FrecAng * t), (sin(AngIn/2)^2))
}

#Function that applies the previous formula using the different intervals (t).
Test <- function(L,g,AngIn)
{
  FrecAng <- 0
  t = seq(from= 0, to= 10, by= 0.1)
  PosAng <- PositionAng(t,AngIn,FrecAng)
  return (PosAng)
}

#Evolution of the angle as a function of time
EvolAng <- Test(L,g,AngIn)

plot(EvolAng,
     type = "l",
     main = "Evolution of the angle in time",
     ylab = "Angle",
     xlab = "Time
)

sn is an elliptical function of Jacobi.

I want to obtain something like the one shown below

.

Answer 1

As Oliver points out, your problem is that you are setting FrecAng to 0 inside Test. Since PositionAng only uses t once, where it is multiplied by FrecAng (i.e. 0), this term is always 0, so this function produces a constant. I think it should be:

Test <- function(L,g,AngIn)
{
  FrecAng <- FrequencyAng(g, L)
  t = seq(from = 0, to = 10, by = 0.1)
  PosAng <- PositionAng(t, AngIn, FrecAng)
  return (PosAng)
}

Which produces:

EvolAng <- Test(L, g, AngIn)

plot(EvolAng,
     type = "l",
     main = "Evolution of the angle in time",
     ylab = "Angle",
     xlab = "Time")

This seems to give sensible results - for a longer pendulum we have what looks like simple harmonic motion:

plot(Test(20, 9.8, pi/4),
     type = "l",
     main = "Evolution of the angle in time",
     ylab = "Angle",
     xlab = "Time")

But for shorter lengths we see the non-linear effects:

plot(Test(0.4, 9.8, pi/4),
     type = "l",
     main = "Evolution of the angle in time",
     ylab = "Angle",
     xlab = "Time")