Я пытаюсь построить график изменения угла θ и угловой скорости ω относительно времени 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 вручную, и, кажется, он ведет себя так, как должен, но график просто ничего мне не дает:
Я знаю, что у него должна быть синусоида, поэтому я думаю, что это может быть проблема с тем, как я пытаюсь сгенерировать график.
Любая помощь приветствуется.
Спасибо за полный вопрос с воспроизводимым примером. Редко можно увидеть это от нового пользователя, и это настоящее удовольствие. Я напишу вам полный ответ на мгновение.
Ура, спасибо за помощь, может быть, это так же хорошо, что он ничего не замышлял, ха-ха
Добро пожаловать в MatPlotLib (и кажется, в Python в целом). Этот ответ будет скорее списком указателей, чем чем-либо еще. Надеюсь, это поможет вам исследовать правильные темы в следующий раз, когда вы решите сделать что-то подобное.
Основная проблема с вашим кодом заключается в том, что вы вызываете plt.plot
для каждой отдельной точки вашей временной эволюции. Это каждый раз создает новый график с парой точек. Вероятно, вы захотите собрать списки или массивы, содержащие 1000 значений t
, omega
и theta
. Затем вы можете построить все сразу вне цикла.
Прежде чем я покажу, как это сделать, я рассмотрю несколько более мелких проблем:
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 есть множество отличных математических библиотек, которые ждут своего изучения.
В принципе ваша версия тоже работает. Только режим рисования по умолчанию plot
заключается в рисовании линий без маркеров для точек выборки. Поскольку одна точка не образует линию, ничто не может быть нарисовано.
Измените команды сюжета на
plt.plot(t, theta,'.b')
plt.plot(t, omega,'.r')
включить маркеры, то будет что-то нарисовано.
Вам все равно придется удалить или изменить настройки axis
.
Этот метод довольно медленный, поскольку существует довольно много, а также много избыточных вызовов библиотеки построения графиков.
Уберите все команды сюжета из цикла. Вы создаете новую диаграмму на каждой итерации.