""" Hybrid Dynamical System: Bouncing Ball in 1-D This model demonstrates a hybrid system with: - Flow State 1: Free fall (when ball center of mass is above radius r) - Flow State 2: Spring-mass-damper (when ball is in contact with ground) - Discrete transitions between these states """ import numpy as np import matplotlib.pyplot as plt from scipy.integrate import solve_ivp class HybridBouncingBall: def __init__(self, m=0.1, r=0.1, g=9.81, k=5000.0, c=10.0): """ Parameters: ----------- m : float Mass of the ball (kg) r : float Radius of the ball (m) g : float Gravitational acceleration (m/s^2) k : float Spring constant when in contact with ground (N/m) c : float Damping coefficient when in contact with ground (N·s/m) """ self.m = m self.r = r self.g = g self.k = k self.c = c # For tracking state transitions self.state_history = [] self.transition_times = [] def free_fall_dynamics(self, t, y): """ Flow dynamics for free fall state. State: y = [position, velocity] """ pos, vel = y dpos = vel dvel = -self.g return [dpos, dvel] def spring_damper_dynamics(self, t, y): """ Flow dynamics for spring-mass-damper state. State: y = [position, velocity] When ball is compressed against ground: F = -k*(r - pos) - c*vel - m*g """ pos, vel = y dpos = vel # Spring force kicks in when position < r # Compression is (r - pos) compression = self.r - pos spring_force = self.k * compression damping_force = self.c * vel dvel = (spring_force - damping_force - self.m * self.g) / self.m return [dpos, dvel] def event_contact_ground(self, t, y): """Event: Ball contacts ground (transition to spring-damper)""" pos, vel = y return pos - self.r def event_leave_ground(self, t, y): """Event: Ball leaves ground (transition to free fall)""" pos, vel = y # Leave ground when position > r AND velocity > 0 if pos > self.r and vel > 0: return 0 return 1 # Make events terminal to stop integration event_contact_ground.terminal = True event_leave_ground.terminal = True def simulate(self, y0, t_span, max_transitions=20): """ Simulate the hybrid system. Parameters: ----------- y0 : list Initial state [position, velocity] t_span : tuple Time span (t_start, t_end) max_transitions : int Maximum number of state transitions to simulate Returns: -------- t_all : array Time points y_all : array State trajectory states : list State labels ('free_fall' or 'spring_damper') """ t_all = [] y_all = [] states = [] current_state = "free_fall" if y0[0] > self.r else "spring_damper" current_y = y0 current_t = t_span[0] t_end = t_span[1] transitions = 0 while current_t < t_end and transitions < max_transitions: if current_state == "free_fall": # Integrate free fall until contact with ground sol = solve_ivp( self.free_fall_dynamics, [current_t, t_end], current_y, events=self.event_contact_ground, dense_output=True, max_step=0.01, ) # Store results t_all.append(sol.t) y_all.append(sol.y.T) states.extend(["free_fall"] * len(sol.t)) # Check if event occurred if sol.t_events[0].size > 0: # Transition to spring-damper current_state = "spring_damper" current_t = sol.t[-1] current_y = sol.y[:, -1] self.transition_times.append(current_t) transitions += 1 else: break else: # spring_damper # Integrate spring-damper until leaving ground sol = solve_ivp( self.spring_damper_dynamics, [current_t, t_end], current_y, events=self.event_leave_ground, dense_output=True, max_step=0.01, ) # Store results t_all.append(sol.t) y_all.append(sol.y.T) states.extend(["spring_damper"] * len(sol.t)) # Check if event occurred if sol.t_events[0].size > 0: # Transition to free fall current_state = "free_fall" current_t = sol.t[-1] current_y = sol.y[:, -1] self.transition_times.append(current_t) transitions += 1 else: break # Concatenate all results t_all = np.concatenate(t_all) y_all = np.vstack(y_all) return t_all, y_all, states def plot_simulation(t, y, states, ball, show_phase=True): """ Plot the simulation results. Parameters: ----------- t : array Time points y : array State trajectory states : list State labels ball : HybridBouncingBall Ball object show_phase : bool Whether to show phase portrait """ # Convert states to numeric for coloring state_numeric = np.array([1 if s == "free_fall" else 2 for s in states]) if show_phase: fig, axes = plt.subplots(2, 1, figsize=(8, 10)) else: fig, axes = plt.subplots(2, 1, figsize=(12, 8)) axes = axes.reshape(-1, 1) # Plot 1: Position vs Time ax1 = axes[0] if show_phase else axes[0] scatter = ax1.scatter(t, y[:, 0], c=state_numeric, s=1, cmap="coolwarm", alpha=0.6) ax1.axhline( y=ball.r, color="k", linestyle="--", label=f"Ground contact (h={ball.r}m)", linewidth=1, ) ax1.axhline( y=0, color="gray", linestyle="-", label="Ground level", linewidth=1, alpha=0.5 ) # Mark transitions for t_trans in ball.transition_times: ax1.axvline(x=t_trans, color="green", linestyle=":", alpha=0.3, linewidth=1) ax1.set_xlabel("Time (s)", fontsize=11) ax1.set_ylabel("Position (m)", fontsize=11) ax1.set_title("Ball Position vs Time", fontsize=12, fontweight="bold") ax1.grid(True, alpha=0.3) ax1.legend(fontsize=9) if show_phase: # Plot 2: Vector Field / Phase Portrait ax2 = axes[1] # Create vector field grid pos_range = np.linspace(0, max(y[:, 0]) * 1.1, 15) vel_range = np.linspace(min(y[:, 1]) * 1.1, max(y[:, 1]) * 1.1, 15) Pos, Vel = np.meshgrid(pos_range, vel_range) # Calculate vector field dPos = np.zeros_like(Pos) dVel = np.zeros_like(Vel) for i in range(Pos.shape[0]): for j in range(Pos.shape[1]): pos_val = Pos[i, j] vel_val = Vel[i, j] # Determine which dynamics to use if pos_val > ball.r: # Free fall dynamics dPos[i, j] = vel_val dVel[i, j] = -ball.g else: # Spring-damper dynamics dPos[i, j] = vel_val compression = ball.r - pos_val spring_force = ball.k * compression damping_force = ball.c * vel_val dVel[i, j] = (spring_force - damping_force - ball.m * ball.g) / ball.m # Plot vector field with much smaller scale ax2.quiver(Pos, Vel, dPos, dVel, alpha=0.3, color='gray', scale=300, width=0.003) # Plot trajectory scatter2 = ax2.scatter( y[:, 0], y[:, 1], c=state_numeric, s=2, cmap="coolwarm", alpha=0.7 ) ax2.axvline( x=ball.r, color="k", linestyle="--", label=f"Contact threshold", linewidth=1.5 ) ax2.axhline(y=0, color="gray", linestyle="-", linewidth=1, alpha=0.5) ax2.set_xlabel("Position (m)", fontsize=11) ax2.set_ylabel("Velocity (m/s)", fontsize=11) ax2.set_title("Vector Field & Phase Portrait", fontsize=12, fontweight="bold") ax2.grid(True, alpha=0.3) ax2.legend(fontsize=9) plt.tight_layout() return fig if __name__ == "__main__": # Create ball with specific parameters ball = HybridBouncingBall( m=0.10, # 1 kg mass r=0.1, # 10 cm radius g=9.81, # Earth gravity k=500.0, # Spring constant c=0.89, # Damping coefficient ) # Initial conditions: drop from 2 meters with zero velocity y0 = [1.0, 0.0] # [position (m), velocity (m/s)] # Simulate for 5 seconds t_span = (0, 5.0) print("Simulating hybrid bouncing ball system...") print(f"Initial conditions: h0 = {y0[0]} m, v0 = {y0[1]} m/s") print( f"Ball parameters: m={ball.m} kg, r={ball.r} m, k={ball.k} N/m, c={ball.c} N·s/m" ) print() t, y, states = ball.simulate(y0, t_span, max_transitions=30) print(f"Simulation complete!") print(f"Total time simulated: {t[-1]:.3f} s") print(f"Number of state transitions: {len(ball.transition_times)}") print(f"Transition times: {[f'{tt:.3f}' for tt in ball.transition_times[:10]]}") print() # Count time in each state free_fall_count = states.count("free_fall") spring_damper_count = states.count("spring_damper") total_points = len(states) print(f"Time distribution:") print(f" Free fall: {free_fall_count/total_points*100:.1f}%") print(f" Spring-damper: {spring_damper_count/total_points*100:.1f}%") # Plot results fig = plot_simulation(t, y, states, ball, show_phase=True) plt.savefig( "/home/danesabo/Documents/Dane's Vault/Presentations/ERLM/bouncing_ball_hybrid.png", dpi=300, bbox_inches="tight", ) print(f"\nPlot saved to: bouncing_ball_hybrid.png") plt.show()