2

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

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

任何帮助将非常感激。

4

2 回答 2

2

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

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

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

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

于 2018-04-13T20:57:53.370 回答
1

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

将绘图命令更改为

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

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

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

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

于 2018-04-14T09:16:33.597 回答