使用梯形法则在 python 中绘制钟摆运动

Question

我正在尝试绘制角度 θ 和 angular 速度 ω 如何使用梯形规则求解微分方程的线性和非线性摆相对于时间 t 的变化，但我在生成实际情节时遇到问题。

这是我尝试为没有摩擦阻尼或驱动力的线性摆实现的代码，其中 θ 初始化为 0.2，ω 初始化为 0.0

import matplotlib.pylab as plt
import math

theta = 0.2 #angle
omega = 0.0 #angular velocity
t = 0.0 #time
dt = 0.01
nsteps = 0
k = 0.0 #dampening coefficient
phi = 0.66667 #angular frequency of driving force
A = 0.0 #amplitude of driving force

#2nd order ODE for linear pendulum
def f(theta, omega, t):
    return -theta - k*omega + A*math.cos(phi*t) 

#trapezoid rule
for nsteps in range(0,1000):
    k1a= dt*omega 
    k1b = dt*f(theta, omega, t) 
    k2a = dt*(omega + k1b)
    k2b = dt*f(theta + k1a, omega + k1b, t + dt)

    theta = theta + (k1a + k2a)/2 
    omega = omega + (k1b + k2b)/2 
    t = t + dt 
    nsteps = nsteps + 1

    plt.plot(t, theta)
    plt.plot(t, omega)
    plt.axis([0, 500, -math.pi, math.pi])
    plt.title('theta = 0.2, omega = 0.0')

plt.show()

我已经通过 for 循环手动计算出前几次迭代的值，它似乎表现得像它应该做的那样，但情节什么也没给我：

我知道它应该得到一个正弦曲线，所以我认为这可能只是我尝试生成绘图的方式的问题。

如有任何帮助，我们将不胜感激。

Answer 1

欢迎使用 MatPlotLib（Python 总的来说）。这个答案将是一个指针列表，而不是其他任何东西。希望它能在您下次决定做这样的事情时引导您研究正确的主题。

您的代码的主要问题是您为时间演变中的每个点调用 plt.plot。这每次都会创建一个包含一对点的新图。您可能想要做的是累积包含 t、omega 和 theta 的 1000 个值的列表或数组。然后，您可以在循环外一次绘制所有这些。

在我展示如何做到这一点之前，这里有一些我看到的小问题：

pylab 在这一点上几乎已被弃用。请改用 pyplot。导入 from matplotlib import pyplot as plt 和 import matplotlib.pyplot as plt 是等价的。
你在循环中递增 nsteps。这是完全没有必要的。 for 循环的工作方式是在每次迭代中分配一个新值 nsteps，取自您正在循环的可迭代对象中的下一个值，在本例中为 range(1000)。当循环将 100 分配给 nsteps，并且您在循环体中将 999 分配给它时，下一次迭代将看到 nsteps 设置为 101不管你喜不喜欢
您在调用 axis 时没有使用正确的限制：dt 为 0.01，1000 步将使您的时间范围为 10 秒，而不是 500。Matplotlib 很漂亮擅长自己设置轴，所以我更愿意完全省略该调用，或者至少将其固定为 plt.axis([0, dt * nsteps, -math.pi, math.pi]).

考虑到所有这些，下面是我将如何重写代码的循环部分：

from matplotlib import pyplot as plt
...

t_list = [t]
omega_list = [omega]
theta_list = [theta]

#trapezoid rule
for nsteps in range(0,1000):
    k1a = dt * omega 
    k1b = dt * f(theta, omega, t) 
    k2a = dt * (omega + k1b)
    k2b = dt * f(theta + k1a, omega + k1b, t + dt)

    theta = theta + (k1a + k2a) / 2 
    omega = omega + (k1b + k2b) / 2 
    t = t + dt 
    t_list.append(t)
    theta_list.append(theta)
    omega_list.append(omega)

plt.plot(t_list, theta_list)
plt.plot(t_list, omega_list)
plt.title('theta = 0.2, omega = 0.0')

plt.show()

上面的代码有效，但效率不高。我会研究 numpy 库作为 Python 中数值计算的开始。它提供了一个数组类型和大量的数学函数。 Numpy 数组也是 MatPlotLib 构建的基础。您传入的所有列表都在内部转换为数组。当您达到 numpy 可以为您做的事情的极限时，请查看 scipy 和其他库。 Python 有很多很棒的数学库等待探索。

Answer 2

原则上你的版本也可以。只有plot的默认绘制模式是为样本点绘制不带标记的线条。一点不成线，什么也画不出来

将绘图命令更改为

plt.plot(t, theta,'.b')
plt.plot(t, omega,'.r')

打开标记，然后就会画一些东西。

您仍然需要删除或调整 axis 设置。

这种方法相当慢，因为对绘图库的调用相当多，而且还有很多冗余。

使用梯形法则在 python 中绘制钟摆运动

Plotting pendulum motion in python using trapezoid rule

python

differential-equations

pendulum