ODE 解的 Matplotlib 图与 RHS 向量场不相切

Matplotlib plot of ODE solution is not tangential to RHS vector field

我正在试验常微分方程,为此我编写了一个小代码来绘制相位 space 中摆的轨迹。然而,轨迹并不像预期的那样真正与 ODE 右侧给出的矢量场相切。这是一个密谋问题吗?当我使用 ODE 的相同函数进行矢量场的数值积分和绘图时,我不知道这是怎么发生的...

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Constants
g = 9.81
l = 1
delta = 0.6

# Initial state vector: [alpha, dalpha_dot]
x0 = np.array([-6, 8])

# Right hand side od pendulum ODE
def pendulum_rhs(t, y):
    return np.array([y[1], g/l * np.sin(y[0]) - delta * y[1]])

# Solve ODE using runge kutta method
ode_sol = solve_ivp (pendulum_rhs, [0, 10], x0, dense_output=True, t_eval=np.linspace(0 ,10, 200))

ode_sol_y = ode_sol.y.T

# Get the vector field along the trajectory
draw_arrow_every_nth = 5
vector_field_at_ode_sol_y = np.array([pendulum_rhs(0, curr_y) for curr_y in ode_sol_y[::draw_arrow_every_nth, :]])

# create whole vector field
alpha = np.linspace(-2*np.pi, 2*np.pi, 21)
alpha_dot = np.linspace(-10, 10, 21)
# full coorindate arrays
alpha_2dgrid, alpha_dot_2dgrid = np.meshgrid(alpha, alpha_dot)

alpha_dt = alpha_dot_2dgrid
alpha_dot_dt = g/l * np.sin(alpha_2dgrid) - delta * alpha_dot_2dgrid

# Plotting
plt.figure()
plt.quiver(alpha_2dgrid, alpha_dot_2dgrid,
       alpha_dt, alpha_dot_dt)
plt.plot(ode_sol_y[:,0], ode_sol_y[:,1])
plt.quiver(ode_sol_y[::draw_arrow_every_nth, 0], ode_sol_y[::draw_arrow_every_nth, 1], vector_field_at_ode_sol_y[:,0], vector_field_at_ode_sol_y[:,1])

这看起来像是一个 aspect 问题。您可以添加 ax.set_aspect('equal') 以获得您想要的结果。

查看下面的代码:

# Plotting
fig,ax=plt.subplots(figsize=(8,8))
ax.plot(ode_sol_y[:,0], ode_sol_y[:,1])
plt.quiver(ode_sol_y[::draw_arrow_every_nth, 0], ode_sol_y[::draw_arrow_every_nth, 1], vector_field_at_ode_sol_y[:,0], vector_field_at_ode_sol_y[:,1])
plt.quiver(alpha_2dgrid, alpha_dot_2dgrid,alpha_dt, alpha_dot_dt)
ax.set_aspect('equal')

并且输出给出:

感谢@jylls 为我指明了正确的方向!

我找到了解决方案,当纵横比不是 1 时,矢量场指向正确的方向 post:Quiver plot arrow aspect ratio

你必须做的事情:

pylab.quiver(x, y, u, v, angles='xy', scale_units='xy', scale=10)