phase portrait - Nonlinear damped pendulum

[Giovanni_62]

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

# Parameters
g = 9.81       # Acceleration due to gravity (m/s^2)
L = 1.0        # Length of the pendulum (m)
b = 0.1        # Damping coefficient

# Define the ODEs
def pendulum_ode(t, state):
    theta, omega = state
    dtheta_dt = omega
    domega_dt = - (g / L) * np.sin(theta) - b * omega
    return [dtheta_dt, domega_dt]

# Create a grid of initial conditions in the phase space
theta_range = np.linspace(-2 * np.pi, 2 * np.pi, 3)
omega_range = np.linspace(-10, 10, 7)

# Initialize lists to store trajectories
trajectories = []

# Integrate the ODEs for a subset of initial conditions
for theta_0 in theta_range:
    for omega_0 in omega_range:
        initial_state = [theta_0, omega_0]
        t_span = (0, 10)  # Time span for integration
        sol = solve_ivp(pendulum_ode, t_span, initial_state, t_eval=np.linspace(*t_span, 1000))
        trajectories.append(sol.y)

# Plot the phase portrait with restricted x and y limits
plt.figure(figsize=(8, 6))
for traj in trajectories:
    plt.plot(traj[0], traj[1], linewidth=0.5)
plt.xlim(-2 * np.pi, 2 * np.pi)  # Restrict x-axis limits
plt.ylim(-10, 10)  # Restrict y-axis limits
plt.xlabel('Angle (radians)')
plt.ylabel('Angular Velocity (radians/s)')
plt.title('Phase Portrait of Damped Nonlinear Pendulum')
plt.grid(True)
plt.show()

Someone help me improve the program, so that we have a visualization of the various trajectories, for example as in

Plot phase diagram

[deanhystad]

Looks interesting:

https://stackoverflow.com/questions/5870...oordinates

http://be150.caltech.edu/2018/handouts/l...raits.html

https://pysd-cookbook.readthedocs.io/en/...raits.html