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

I am experimenting with oridnary differential equations and for that I programmed a little code that plots the trajectory of a pendulum in phase space. However, the trajectory isn't really tangential to the vector field given by the right hand side of the ODE as it is supposed to be. Is this a plotting issue? As I use the same function of the ODE for both the numerical integration and plotting of the vector field I dont know how this can happen...

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])

CodePudding user response：

This looks like an aspect issue. You could add ax.set_aspect('equal') to get the result you want.

See code below:

# 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')

And the output gives:

CodePudding user response：

Thank you @jylls for pointing me in the right direction!

I have found the solution for having the vector field point in the right directions when the aspect ratio isn't 1 in this post: Quiver plot arrow aspect ratio

What you have to do:

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