r - R中的非线性摆：获得角度的精确时间演化的问题

Question

我正在研究一个简单的钟摆问题，以获得在知道数据集后的时间间隔内角度 θ 的演变。为了得出解决方案，我基于Meléndez 等人的工作，得出了 θ(t) 时间演化的以下精确表达式。

在 RI 的哪里有这种方式：

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

函数 θ(t) 值的计算变得越来越复杂，当我运行我的程序时，我得到一个基于 Melendez 工作的非常不同的图表。

我的图表没有椭圆曲率并且是线性的。

我想知道我的代码是否有错误，或者我可以通过其他方式解决问题。

从角度为：pi/4（45度）、绳索长度为0.8m、重力为9.8m/s^2的情况开始。

我与图书馆（椭圆）合作找到了这样的结果，但我不知道我的图表会发生什么。

我已经构建了以下代码：

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是 Jacobi 的椭圆函数。

我想获得如下所示的东西

.

Answer 1

正如奥利弗指出的那样，您的问题是您将FrecAng0 inside 设置为Test。由于PositionAng只使用t一次，乘以FrecAng（即 0），此项始终为 0，因此该函数产生一个常数。我认为应该是：

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

产生：

EvolAng <- Test(L, g, AngIn)

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

这似乎给出了合理的结果——对于更长的钟摆，我们有看起来像简谐运动：

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

但是对于较短的长度，我们会看到非线性效应：

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