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

【问题标题】：Plotting pendulum motion in python using trapezoid rule使用梯形规则在 python 中绘制钟摆运动
【发布时间】：2018-04-13 19:09:46
【问题描述】：

我正在尝试绘制角度 θ 和角速度 ω 如何随时间 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 循环手动计算了前几次迭代的值，它似乎表现得应该做，但情节什么也没给我：

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

任何帮助将不胜感激。

【问题讨论】：

将所有绘图命令移出循环。您正在每次迭代中创建一个新图表。
感谢您提供带有可重复示例的完整问题。很少能从新用户那里看到这一点，而且真的很高兴。我马上给你写一个完整的答案。
干杯感谢您的帮助，也许它没有在策划任何事情哈哈

标签： python differential-equations pendulum

【解决方案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 有很多很棒的数学库等待探索。

【讨论】：

【解决方案2】：

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

将绘图命令更改为

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

打开标记，然后会绘制一些东西。

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

此方法相当慢，因为对绘图库的调用相当多且冗余。

【讨论】：