Построение движения маятника в питоне с использованием правила трапеции

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 в каждой итерации, взятой из следующего значения в итерируемом цикле, в этом случае range(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. Он предоставляет тип массива и множество математических функций. Матричные массивы - это основа того, на чем построен MatPlotLib. Все списки, которые вы передаете, конвертируются во внутренние массивы. Когда вы достигнете предела того, что numpy может сделать для вас, посмотрите на scipy и другие библиотеки. В Python есть много замечательных математических библиотек, ожидающих изучения.

Answer 2

В принципе ваша версия тоже работает. Только режим рисования по умолчанию plot рисовать линии без маркеров для точек выборки. Поскольку одна точка не делает линию, ничто не может быть нарисовано.

Измените команды графика на

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

чтобы включить маркеры, тогда будет что-то нарисовано.

Вам все равно придется удалить или адаптировать axis Настройки.

Этот метод довольно медленный, так как существует довольно много, а также много избыточных вызовов библиотеки заговоров.