Spring-Mass System Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# System parameters
m = 1.0      # Mass (kg)
k = 10.0     # Spring constant (N/m)
c = 0.5      # Damping coefficient (N⋅s/m)

# Initial conditions
x0 = 0.1     # Initial displacement (m)
v0 = 0.0     # Initial velocity (m/s)

# Time parameters
dt = 0.01
t_max = 10
time = np.arange(0, t_max, dt)

# Calculate natural frequency
omega_n = np.sqrt(k/m)
print(f"Natural frequency: {omega_n:.2f} rad/s")
print(f"Natural period: {2*np.pi/omega_n:.2f} s")

def spring_mass_ode(x, v, t, force=0):
    """Spring-mass system equations"""
    dx_dt = v
    dv_dt = (force - c*v - k*x) / m
    return dx_dt, dv_dt

# Solve the system
x = np.zeros(len(time))
v = np.zeros(len(time))

x[0] = x0
v[0] = v0

for i in range(len(time)-1):
    # External force (can be modified)
    F_external = 0  # Free vibration

    dx, dv = spring_mass_ode(x[i], v[i], time[i], F_external)
    x[i+1] = x[i] + dx * dt
    v[i+1] = v[i] + dv * dt

# Create visualization
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# Plot displacement vs time
ax1.plot(time, x)
ax1.set_ylabel('Displacement (m)')
ax1.set_title('Spring-Mass System Response')
ax1.grid(True)

# Plot phase portrait (velocity vs displacement)
ax2.plot(x, v)
ax2.set_xlabel('Displacement (m)')
ax2.set_ylabel('Velocity (m/s)')
ax2.set_title('Phase Portrait')
ax2.grid(True)

plt.tight_layout()
plt.show()