如何在Python中得到哈密顿谐振子的相图？

import numpy as np
import matplotlib.pyplot as plt

def hamiltonian_equations(q, p, omega):
    dqdt = p
    dpdt = -omega**2 * q
    return dqdt, dpdt

q0 = 1.0  # 初始位置
p0 = 0.0  # 初始动量
omega = 1.0  # 振子的频率
dt = 0.01  # 积分步长
num_steps = 1000  # 积分步数

q = np.zeros(num_steps)
p = np.zeros(num_steps)

q[0] = q0
p[0] = p0

for i in range(num_steps-1):
    dqdt, dpdt = hamiltonian_equations(q[i], p[i], omega)
    q[i+1] = q[i] + dqdt * dt
    p[i+1] = p[i] + dpdt * dt

plt.plot(q, p)
plt.xlabel('Position')
plt.ylabel('Momentum')
plt.title('Phase Portrait of Harmonic Oscillator')
plt.show()