Sep-08-2023, 07:42 PM
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