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

Question

我正在尝试绘制角度 θ 和角速度 ω 相对于线性和非线性摆的时间 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。这每次都会创建一个包含一对点的新图。您可能想要做的是累积包含 1000 个 t、omega 和 theta 值的列表或数组。然后，您可以在循环外一次绘制所有这些。

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

1. pylab 在这一点上几乎已被弃用。请改用 pyplot。导入 from matplotlib import pyplot as plt 和 import matplotlib.pyplot as plt 是等价的。
2. 您在循环中递增 nsteps。这是完全没有必要的。 for 循环的工作方式是在每次迭代中分配一个新值 nsteps，取自您正在循环的可迭代对象中的下一个值，在本例中为 范围(1000)。当循环将 100 分配给 nsteps，并且您在循环体中将 999 分配给它时，下一次迭代将看到 nsteps 设置为 101 不管你喜不喜欢。
3. 您在调用 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 有很多很棒的数学库有待探索。