GNUPLOT 错误拉格朗日双摆 GIF 生成器的最后一行没有以前的文件名

Question

晚安。我花了 5 个多小时尝试制作拉格朗日双摆的简单 GIF。数学没问题，一切都很好，但我无法生成 gif，因为每次我尝试加载时都会遇到同样的错误，而且在互联网上我找不到解决方案。下面是我的代码：

# Constants
m1 = 1.0
m2 = 1.0
L1 = 1.0
L2 = 1.0
g = 9.8

# Initial conditions
theta1_0 = 0.0
theta2_0 = 0.0
theta1_dot_0 = 0.0
theta2_dot_0 = 0.0

# Time
dt = 0.01
t_max = 10.0

# Initial values of theta1, theta2, theta1_dot, and theta2_dot
theta1 = theta1_0
theta2 = theta2_0
theta1_dot = theta1_dot_0
theta2_dot = theta2_dot_0

# Set output file
set term gif animate delay 10 optimize optimize
set output "double_pendulum.gif"

# Numerical integration using the 4th order Runge-Kutta method
do for [i = 1:int(t_max / dt)] {
    # Calculate derivatives
    dtheta1_dt = theta1_dot
    dtheta2_dt = theta2_dot
    dtheta1_dot_dt = (m2 * g * sin(theta2) * cos(theta1 - theta2) - m2 * sin(theta1 - theta2) * (L1 * theta1_dot**2 * cos(theta1 - theta2) + L2 * theta2_dot**2) - (m1 + m2) * g * sin(theta1)) / (L1 * (m1 + m2 * sin(theta1 - theta2)**2))
    dtheta2_dot_dt = (( (m1 + m2) * (L1 * theta1_dot**2 * sin(theta1 - theta2) - g * sin(theta2) + g * sin(theta1) * cos(theta1 - theta2)) + m2 * L2 * theta2_dot**2 * sin(theta1 - theta2) * cos(theta1 - theta2)) / (L2 * (m1 + m2 * sin(theta1 - theta2)**2)))

    # Update the values using the 4th order Runge-Kutta method
    k1_theta1 = dt * dtheta1_dt
    k1_theta2 = dt * dtheta2_dt
    k1_theta1_dot = dt * dtheta1_dot_dt
    k1_theta2_dot = dt * dtheta2_dot_dt

    k2_theta1 = dt * (dtheta1_dt + 0.5 * k1_theta1_dot)
    k2_theta2 = dt * (dtheta2_dt + 0.5 * k1_theta2_dot)
    k2_theta1_dot = dt * (dtheta1_dot_dt + 0.5 * dtheta1_dot_dt)
    k2_theta2_dot = dt * (dtheta2_dot_dt + 0.5 * dtheta2_dot_dt)

    k3_theta1 = dt * (dtheta1_dt + 0.5 * k2_theta1_dot)
    k3_theta2 = dt * (dtheta2_dt + 0.5 * k2_theta2_dot)
    k3_theta1_dot = dt * (dtheta1_dot_dt + 0.5 * dtheta1_dot_dt)
    k3_theta2_dot = dt * (dtheta2_dot_dt + 0.5 * dtheta2_dot_dt)

    k4_theta1 = dt * (dtheta1_dt + k3_theta1_dot)
    k4_theta2 = dt * (dtheta2_dt + k3_theta2_dot)
    k4_theta1_dot = dt * (dtheta1_dot_dt + dtheta1_dot_dt)
    k4_theta2_dot = dt * (dtheta2_dot_dt + dtheta2_dot_dt)

    theta1_dot = theta1_dot + (1.0 / 6.0) * (k1_theta1_dot + 2.0 * k2_theta1_dot + 2.0 * k3_theta1_dot + k4_theta1_dot)
    theta2_dot = theta2_dot + (1.0 / 6.0) * (k1_theta2_dot + 2.0 * k2_theta2_dot + 2.0 * k3_theta2_dot + k4_theta2_dot)
    theta1 = theta1 + (1.0 / 6.0) * (k1_theta1 + 2.0 * k2_theta1 + 2.0 * k3_theta1 + k4_theta1)
    theta2 = theta2 + (1.0 / 6.0) * (k1_theta2 + 2.0 * k2_theta2 + 2.0 * k3_theta2 + k4_theta2)

    # Update time
    t = i * dt

    # Plot the double pendulum motion
    if (i == 1) {
        set xrange [-2:2]
        set yrange [-2:2]
        plot L1 * sin(theta1), -L1 * cos(theta1), L1 * sin(theta1) + L2 * sin(theta2), -L1 * cos(theta1) - L2 * cos(theta2) with lines title 'Pendulum 1', '' using 3:4 with lines title 'Pendulum 2'
    } else {
        replot
    }
    pause 0.01
}

我正在使用 gnuplot 5.4 Patch 8 并使用四阶 Runge-Kutta 方法来求解二维模型的微分方程。

这是错误的图片。我不知道该怎么办了。我不是 gnuplot 方面的专家。但我还做了一些其他的事情。

Answer 1

错误来自于您的初始绘图命令，该命令是在循环中第一次执行的。为了便于阅读，将命令拆分为多行，它是：

plot L1 * sin(theta1), \
    -L1 * cos(theta1), \
     L1 * sin(theta1) + L2 * sin(theta2), \
    -L1 * cos(theta1) - L2 * cos(theta2) with lines title 'Pendulum 1', \
    '' using 3:4 with lines title 'Pendulum 2'

问题是最后一个情节条款

'' using 3:4

。空引号表示“先前的文件名”，但实际上没有先前的文件名。我不清楚你在那里的意图，但是

''

无法工作。

您是否可能在某个文件中保存了以前的拟合，并且这就是您想要绘制进行比较的内容？将该文件名放在引号内。

Answer 2

我试图让你的脚本运行。 实际上，在我的评论中我错了：你的 thetas 是 时间的函数，但是你的线

t=i*dt

毫无用处，因为

t

没有被使用。

一些评论：

• 使用

  set size ratio -1

  使 x 和 y 绘图尺寸相同
• 如果您的初始值

  theta1_0, theta2_0, theta1_dot_0, theta2_dot_0

  全部为零，则什么都不会移动。因此，至少有一个值应该不等于零。
• 如果你想绘制钟摆，你最好使用绘图风格

  with vectors

  ，如果你也想绘制质量，另外，绘图风格

  with points

• plot '+' every ::::0

  用于在循环的每次迭代中绘制单个点（或向量）。

不保证数学正确，我没有检查这一点，只是从你的脚本中接管它。至少，输出看起来有些合理。

脚本：

### gif of double pendulum
reset session

# Constants
m1 = 1.0
m2 = 1.0
L1 = 1.0
L2 = 1.0
g  = 9.81

# Initial conditions
theta1_0 = 30.0
theta2_0 = 30.0
theta1_dot_0 = 1.0
theta2_dot_0 = 1.0

# Time
dt    = 0.02
t_max = 10.0

# Initial values of theta1, theta2, theta1_dot, and theta2_dot
theta1 = theta1_0
theta2 = theta2_0
theta1_dot = theta1_dot_0
theta2_dot = theta2_dot_0

# Set output file
set term gif size 400,400 animate delay 1 optimize
set output "SO76561252.gif"

set xrange [-2:2]
set yrange [-2:2]
set size ratio -1
set key noautotitle

# Numerical integration using the 4th order Runge-Kutta method
do for [i = 0:int(t_max / dt)] {

    # Calculate derivatives
    dtheta1_dt = theta1_dot
    dtheta2_dt = theta2_dot
    dtheta1_dot_dt = (m2 * g * sin(theta2) * cos(theta1 - theta2) - m2 * sin(theta1 - theta2) * (L1 * theta1_dot**2 * cos(theta1 - theta2) + L2 * theta2_dot**2) - (m1 + m2) * g * sin(theta1)) / (L1 * (m1 + m2 * sin(theta1 - theta2)**2))
    dtheta2_dot_dt = (( (m1 + m2) * (L1 * theta1_dot**2 * sin(theta1 - theta2) - g * sin(theta2) + g * sin(theta1) * cos(theta1 - theta2)) + m2 * L2 * theta2_dot**2 * sin(theta1 - theta2) * cos(theta1 - theta2)) / (L2 * (m1 + m2 * sin(theta1 - theta2)**2)))

    # Update the values using the 4th order Runge-Kutta method
    k1_theta1 = dt * dtheta1_dt
    k1_theta2 = dt * dtheta2_dt
    k1_theta1_dot = dt * dtheta1_dot_dt
    k1_theta2_dot = dt * dtheta2_dot_dt

    k2_theta1 = dt * (dtheta1_dt + 0.5 * k1_theta1_dot)
    k2_theta2 = dt * (dtheta2_dt + 0.5 * k1_theta2_dot)
    k2_theta1_dot = dt * (dtheta1_dot_dt + 0.5 * dtheta1_dot_dt)
    k2_theta2_dot = dt * (dtheta2_dot_dt + 0.5 * dtheta2_dot_dt)

    k3_theta1 = dt * (dtheta1_dt + 0.5 * k2_theta1_dot)
    k3_theta2 = dt * (dtheta2_dt + 0.5 * k2_theta2_dot)
    k3_theta1_dot = dt * (dtheta1_dot_dt + 0.5 * dtheta1_dot_dt)
    k3_theta2_dot = dt * (dtheta2_dot_dt + 0.5 * dtheta2_dot_dt)

    k4_theta1 = dt * (dtheta1_dt + k3_theta1_dot)
    k4_theta2 = dt * (dtheta2_dt + k3_theta2_dot)
    k4_theta1_dot = dt * (dtheta1_dot_dt + dtheta1_dot_dt)
    k4_theta2_dot = dt * (dtheta2_dot_dt + dtheta2_dot_dt)

    theta1_dot = theta1_dot + (1.0 / 6.0) * (k1_theta1_dot + 2.0 * k2_theta1_dot + 2.0 * k3_theta1_dot + k4_theta1_dot)
    theta2_dot = theta2_dot + (1.0 / 6.0) * (k1_theta2_dot + 2.0 * k2_theta2_dot + 2.0 * k3_theta2_dot + k4_theta2_dot)
    theta1 = theta1 + (1.0 / 6.0) * (k1_theta1 + 2.0 * k2_theta1 + 2.0 * k3_theta1 + k4_theta1)
    theta2 = theta2 + (1.0 / 6.0) * (k1_theta2 + 2.0 * k2_theta2 + 2.0 * k3_theta2 + k4_theta2)

    # Plot the double pendulum motion
    plot '+' every ::::0 u (0):(0):(x1=L1*sin(theta1)):(y1=-L1*cos(theta1))     w vec nohead lw 2 lc "red"  ti 'Pendulum 1', \
         '+' every ::::0 u (x1):(y1):(dx2=L2*sin(theta2)):(dy2=-L2*cos(theta2)) w vec nohead lw 2 lc "blue" ti 'Pendulum 2', \
         '+' every ::::0 u (x1):(y1)         w p pt 7 ps 2 lc "red", \
         '+' every ::::0 u (x1+dx2):(y1+dy2) w p pt 7 ps 2 lc "blue", \
}
set output
### end of script

结果：

SO76561252.gif

我不知道为什么 GIF 移动得这么慢（尽管最小

delay 1

）。在不同的查看器中，GIF 的移动速度比在浏览器中快得多。

编辑：看起来延迟是否为<20 ms, the animation speed in the browser will be slow, but with 20 ms it will be much faster. So, maybe 20 ms is some minimum delay?