2.4 Tebranishlar — Amaliyot

Laboratoriya ishi: Python da tebranish va rezonans simulyatsiyasi

---

1. Oddiy Garmonik Harakat (SHM)

import numpy as np
import matplotlib.pyplot as plt

def simple_harmonic_motion():
    """SHM simulyatsiyasi va vizualizatsiya"""
    
    # Parametrlar
    A = 0.1      # m (amplituda)
    omega = 5    # rad/s
    phi = 0      # rad (boshlang'ich faza)
    
    T = 2 * np.pi / omega
    t = np.linspace(0, 4*T, 500)
    
    # Harakat
    x = A * np.cos(omega * t + phi)
    v = -A * omega * np.sin(omega * t + phi)
    a = -A * omega**2 * np.cos(omega * t + phi)
    
    # Faza fazosi
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # Vaqt bo'yicha
    axes[0, 0].plot(t, x, 'b-', linewidth=2, label='x(t)')
    axes[0, 0].plot(t, v/omega, 'r--', linewidth=2, label='v(t)/ω')
    axes[0, 0].set_xlabel('Vaqt (s)')
    axes[0, 0].set_ylabel('Kattalik')
    axes[0, 0].set_title('Pozitsiya va Tezlik')
    axes[0, 0].legend()
    axes[0, 0].grid(True)
    
    # Faza fazosi (x vs v)
    axes[0, 1].plot(x, v, 'g-', linewidth=2)
    axes[0, 1].set_xlabel('Pozitsiya (m)')
    axes[0, 1].set_ylabel('Tezlik (m/s)')
    axes[0, 1].set_title('Faza Fazosi')
    axes[0, 1].axis('equal')
    axes[0, 1].grid(True)
    
    # Tezlanish
    axes[1, 0].plot(t, a, 'm-', linewidth=2)
    axes[1, 0].set_xlabel('Vaqt (s)')
    axes[1, 0].set_ylabel('Tezlanish (m/s²)')
    axes[1, 0].set_title('Tezlanish')
    axes[1, 0].grid(True)
    
    # a vs x (chiziqli munosabat)
    axes[1, 1].plot(x, a, 'c-', linewidth=2)
    axes[1, 1].set_xlabel('Pozitsiya (m)')
    axes[1, 1].set_ylabel('Tezlanish (m/s²)')
    axes[1, 1].set_title('a = -ω²x')
    axes[1, 1].grid(True)
    
    plt.tight_layout()
    plt.savefig('shm_analysis.png', dpi=150)
    plt.show()

simple_harmonic_motion()

---

2. So'nuvchi Tebranishlar

import numpy as np
import matplotlib.pyplot as plt

def damped_oscillation(zeta, omega_n, x0=1, v0=0, t_max=10, dt=0.01):
    """
    So'nuvchi tebranish simulyatsiyasi
    
    zeta < 1: Underdamped
    zeta = 1: Critically damped
    zeta > 1: Overdamped
    """
    t = np.arange(0, t_max, dt)
    
    if zeta < 1:
        # Underdamped
        omega_d = omega_n * np.sqrt(1 - zeta**2)
        A = np.sqrt(x0**2 + ((v0 + zeta*omega_n*x0)/omega_d)**2)
        phi = np.arctan2(x0 * omega_d, v0 + zeta*omega_n*x0)
        x = A * np.exp(-zeta*omega_n*t) * np.cos(omega_d*t - phi + np.pi/2)
        
    elif zeta == 1:
        # Critically damped
        x = (x0 + (v0 + omega_n*x0)*t) * np.exp(-omega_n*t)
        
    else:
        # Overdamped
        s1 = -omega_n * (zeta - np.sqrt(zeta**2 - 1))
        s2 = -omega_n * (zeta + np.sqrt(zeta**2 - 1))
        A = (v0 - s2*x0) / (s1 - s2)
        B = (s1*x0 - v0) / (s1 - s2)
        x = A * np.exp(s1*t) + B * np.exp(s2*t)
    
    return t, x

# Turli damping
omega_n = 5
zetas = [0.1, 0.3, 0.7, 1.0, 2.0]
labels = ['ζ=0.1 (underdamped)', 'ζ=0.3', 'ζ=0.7', 
          'ζ=1.0 (critical)', 'ζ=2.0 (overdamped)']

plt.figure(figsize=(12, 6))

for zeta, label in zip(zetas, labels):
    t, x = damped_oscillation(zeta, omega_n)
    plt.plot(t, x, linewidth=2, label=label)

plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.xlabel('Vaqt (s)')
plt.ylabel('Pozitsiya')
plt.title('So\'nuvchi Tebranishlar')
plt.legend()
plt.grid(True)
plt.savefig('damped_oscillations.png', dpi=150)
plt.show()

---

3. Majburiy Tebranish va Rezonans

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

def forced_oscillation(omega_d, omega_n=10, zeta=0.1, F0=1, m=1, t_max=30):
    """
    Majburiy tebranish simulyatsiyasi
    """
    k = omega_n**2 * m
    c = 2 * zeta * omega_n * m
    
    def dynamics(y, t):
        x, v = y
        F = F0 * np.cos(omega_d * t)
        a = (F - c*v - k*x) / m
        return [v, a]
    
    t = np.linspace(0, t_max, 3000)
    y0 = [0, 0]
    sol = odeint(dynamics, y0, t)
    
    return t, sol[:, 0]

# Chastota javob (frequency response)
omega_n = 10
zeta = 0.1
omega_range = np.linspace(0.1, 20, 100)

amplitudes = []
for omega_d in omega_range:
    r = omega_d / omega_n
    X = 1 / np.sqrt((1 - r**2)**2 + (2*zeta*r)**2)
    amplitudes.append(X)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Chastota javob
axes[0, 0].plot(omega_range/omega_n, amplitudes, 'b-', linewidth=2)
axes[0, 0].axvline(x=1, color='r', linestyle='--', label='ω = ωₙ')
axes[0, 0].set_xlabel('ω/ωₙ')
axes[0, 0].set_ylabel('Amplituda (normalized)')
axes[0, 0].set_title('Chastota Javob (Frequency Response)')
axes[0, 0].legend()
axes[0, 0].grid(True)

# Turli chastotalarda vaqt javob
omegas = [5, 10, 15]  # Below, at, above resonance
colors = ['green', 'red', 'blue']
labels = ['ω < ωₙ', 'ω ≈ ωₙ (rezonans)', 'ω > ωₙ']

for omega, color, label in zip(omegas, colors, labels):
    t, x = forced_oscillation(omega, omega_n=10, zeta=0.1, t_max=10)
    axes[0, 1].plot(t, x, color=color, linewidth=1.5, label=label)

axes[0, 1].set_xlabel('Vaqt (s)')
axes[0, 1].set_ylabel('Pozitsiya')
axes[0, 1].set_title('Vaqt Javob (turli ω)')
axes[0, 1].legend()
axes[0, 1].grid(True)

# Turli damping
zetas = [0.05, 0.1, 0.3, 0.7]
for zeta in zetas:
    amps = []
    for omega_d in omega_range:
        r = omega_d / omega_n
        X = 1 / np.sqrt((1 - r**2)**2 + (2*zeta*r)**2)
        amps.append(X)
    axes[1, 0].plot(omega_range/omega_n, amps, linewidth=2, label=f'ζ={zeta}')

axes[1, 0].set_xlabel('ω/ωₙ')
axes[1, 0].set_ylabel('Amplituda')
axes[1, 0].set_title('Damping Ta\'siri')
axes[1, 0].legend()
axes[1, 0].grid(True)
axes[1, 0].set_ylim([0, 12])

# Faza kechikishi
phases = []
for omega_d in omega_range:
    r = omega_d / omega_n
    phi = np.arctan2(2*zeta*r, 1 - r**2)
    phases.append(np.degrees(phi))

axes[1, 1].plot(omega_range/omega_n, phases, 'm-', linewidth=2)
axes[1, 1].axhline(y=90, color='r', linestyle='--')
axes[1, 1].set_xlabel('ω/ωₙ')
axes[1, 1].set_ylabel('Faza (°)')
axes[1, 1].set_title('Faza Kechikishi')
axes[1, 1].grid(True)

plt.tight_layout()
plt.savefig('forced_resonance.png', dpi=150)
plt.show()

---

4. Bog'langan Mayatniklar

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

def coupled_pendulums(k_couple=5, L=1, m=1, g=10):
    """Ikki bog'langan mayatnik"""
    
    omega_0 = np.sqrt(g/L)  # Tabiiy chastota
    omega_c = np.sqrt(k_couple/m)  # Ulanish chastota
    
    def dynamics(y, t):
        theta1, omega1, theta2, omega2 = y
        
        # Kichik burchaklar uchun linearized
        alpha1 = -omega_0**2 * theta1 - omega_c**2 * (theta1 - theta2)
        alpha2 = -omega_0**2 * theta2 - omega_c**2 * (theta2 - theta1)
        
        return [omega1, alpha1, omega2, alpha2]
    
    t = np.linspace(0, 30, 3000)
    
    # Boshlang'ich shart: faqat birinchi mayatnik
    y0 = [0.2, 0, 0, 0]
    
    sol = odeint(dynamics, y0, t)
    
    # Normal modlar chastotalari
    omega_1 = omega_0  # In-phase
    omega_2 = np.sqrt(omega_0**2 + 2*omega_c**2)  # Out-of-phase
    
    print(f"Normal mod 1 (in-phase): f = {omega_1/(2*np.pi):.2f} Hz")
    print(f"Normal mod 2 (out-of-phase): f = {omega_2/(2*np.pi):.2f} Hz")
    print(f"Beat chastotasi: f = {(omega_2-omega_1)/(2*np.pi):.2f} Hz")
    
    # Plot
    fig, axes = plt.subplots(2, 1, figsize=(12, 8))
    
    axes[0].plot(t, sol[:, 0], 'b-', linewidth=1.5, label='θ₁')
    axes[0].plot(t, sol[:, 2], 'r-', linewidth=1.5, label='θ₂')
    axes[0].set_ylabel('Burchak (rad)')
    axes[0].set_title('Bog\'langan Mayatniklar - Energiya Almashinuvi')
    axes[0].legend()
    axes[0].grid(True)
    
    # Energiyalar
    E1 = 0.5 * m * (L * sol[:, 1])**2 + 0.5 * m * g * L * sol[:, 0]**2
    E2 = 0.5 * m * (L * sol[:, 3])**2 + 0.5 * m * g * L * sol[:, 2]**2
    
    axes[1].plot(t, E1, 'b-', linewidth=1.5, label='E₁')
    axes[1].plot(t, E2, 'r-', linewidth=1.5, label='E₂')
    axes[1].plot(t, E1+E2, 'k--', linewidth=1, label='E_total')
    axes[1].set_xlabel('Vaqt (s)')
    axes[1].set_ylabel('Energiya (J)')
    axes[1].set_title('Energiya Almashinuvi')
    axes[1].legend()
    axes[1].grid(True)
    
    plt.tight_layout()
    plt.savefig('coupled_pendulums.png', dpi=150)
    plt.show()

coupled_pendulums()

---

5. Dron Vibratsiya Filtri

import numpy as np
import matplotlib.pyplot as plt

def vibration_filter_demo():
    """Dron IMU vibratsiya filtrlash"""
    
    dt = 0.001  # 1 kHz sampling
    t = np.arange(0, 1, dt)
    
    # Signal: yaxshi signal + propeller vibratsiya
    f_signal = 2  # Hz (dron harakati)
    f_vibration = 100  # Hz (propeller)
    
    signal_clean = np.sin(2*np.pi*f_signal*t)
    vibration = 0.3 * np.sin(2*np.pi*f_vibration*t)
    signal_noisy = signal_clean + vibration
    
    # Low-pass filtrlar
    def low_pass_ema(x, alpha):
        """Exponential Moving Average"""
        y = np.zeros_like(x)
        y[0] = x[0]
        for i in range(1, len(x)):
            y[i] = alpha * x[i] + (1 - alpha) * y[i-1]
        return y
    
    def butterworth_2nd(x, fc, fs):
        """2nd order Butterworth (simplified)"""
        omega_c = 2 * np.pi * fc / fs
        alpha = np.sin(omega_c) / (2 * 0.707)  # Q = 0.707
        
        b0 = (1 - np.cos(omega_c)) / 2
        b1 = 1 - np.cos(omega_c)
        b2 = b0
        a0 = 1 + alpha
        a1 = -2 * np.cos(omega_c)
        a2 = 1 - alpha
        
        y = np.zeros_like(x)
        for i in range(2, len(x)):
            y[i] = (b0*x[i] + b1*x[i-1] + b2*x[i-2] - a1*y[i-1] - a2*y[i-2]) / a0
        return y
    
    # Filtrlar qo'llash
    filtered_ema = low_pass_ema(signal_noisy, alpha=0.1)
    filtered_butter = butterworth_2nd(signal_noisy, fc=20, fs=1/dt)
    
    # Plot
    fig, axes = plt.subplots(3, 1, figsize=(12, 10))
    
    axes[0].plot(t[:500], signal_clean[:500], 'b-', linewidth=2, label='Original')
    axes[0].plot(t[:500], signal_noisy[:500], 'r-', alpha=0.5, label='+ Vibratsiya')
    axes[0].set_ylabel('Signal')
    axes[0].set_title('IMU Signal + Propeller Vibratsiya (100 Hz)')
    axes[0].legend()
    axes[0].grid(True)
    
    axes[1].plot(t[:500], signal_clean[:500], 'b-', linewidth=2, label='Original')
    axes[1].plot(t[:500], filtered_ema[:500], 'g-', linewidth=2, label='EMA Filter')
    axes[1].set_ylabel('Signal')
    axes[1].set_title('EMA Low-Pass Filter')
    axes[1].legend()
    axes[1].grid(True)
    
    axes[2].plot(t[:500], signal_clean[:500], 'b-', linewidth=2, label='Original')
    axes[2].plot(t[:500], filtered_butter[:500], 'm-', linewidth=2, label='Butterworth 2nd')
    axes[2].set_xlabel('Vaqt (s)')
    axes[2].set_ylabel('Signal')
    axes[2].set_title('Butterworth 2nd Order Filter')
    axes[2].legend()
    axes[2].grid(True)
    
    plt.tight_layout()
    plt.savefig('vibration_filter.png', dpi=150)
    plt.show()

vibration_filter_demo()

---

6. PID Tebranish Tahlili

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

def pid_oscillation_analysis():
    """PID tuning va tebranish xarakteristikasi"""
    
    # Plant: 2nd order system
    def closed_loop(y, t, Kp, Ki, Kd, setpoint=1):
        x, v, integral = y
        
        error = setpoint - x
        
        # PID
        P = Kp * error
        I = Ki * integral
        D = Kd * (-v)  # derivative on output
        
        u = P + I + D
        
        # Plant: m*a + c*v + k*x = u
        m, c, k = 1, 0.5, 4
        a = (u - c*v - k*x) / m
        
        return [v, a, error]
    
    t = np.linspace(0, 15, 1500)
    
    # Turli PID parametrlar
    configs = [
        (2, 0, 0, 'P only (oscillates)'),
        (4, 0, 0, 'High P (more oscillation)'),
        (2, 1, 0, 'PI'),
        (2, 0, 1, 'PD'),
        (2, 1, 0.5, 'PID (tuned)'),
    ]
    
    fig, axes = plt.subplots(2, 1, figsize=(12, 8))
    
    for Kp, Ki, Kd, label in configs:
        y0 = [0, 0, 0]
        sol = odeint(closed_loop, y0, t, args=(Kp, Ki, Kd))
        axes[0].plot(t, sol[:, 0], linewidth=2, label=label)
    
    axes[0].axhline(y=1, color='k', linestyle='--', alpha=0.5)
    axes[0].set_ylabel('Pozitsiya')
    axes[0].set_title('PID Tuning - Tebranish Tahlili')
    axes[0].legend(loc='right')
    axes[0].grid(True)
    
    # Optimal tuning: zeta = 0.7
    Kp, Ki, Kd = 2, 1, 0.5
    y0 = [0, 0, 0]
    sol = odeint(closed_loop, y0, t, args=(Kp, Ki, Kd))
    
    axes[1].plot(t, sol[:, 0], 'b-', linewidth=2)
    axes[1].axhline(y=1, color='r', linestyle='--')
    
    # Overshoot va settling time
    overshoot = (np.max(sol[:, 0]) - 1) * 100
    settling_idx = np.where(np.abs(sol[:, 0] - 1) < 0.02)[0]
    settling_time = t[settling_idx[0]] if len(settling_idx) > 0 else None
    
    axes[1].set_xlabel('Vaqt (s)')
    axes[1].set_ylabel('Pozitsiya')
    axes[1].set_title(f'Optimal PID: Overshoot={overshoot:.1f}%, Settling={settling_time:.2f}s')
    axes[1].grid(True)
    
    plt.tight_layout()
    plt.savefig('pid_oscillation.png', dpi=150)
    plt.show()

pid_oscillation_analysis()

---

7. Topshiriq: Robot Qo'li Tebranishi

def robot_arm_vibration():
    """
    Robot qo'li tebranish tahlili
    
    TODO:
    1. 2-bo'g'inli robot qo'li modeli
    2. Tabiiy chastotalarni toping
    3. Rezonansdan qochish strategiyasi
    """
    
    # Parametrlar
    m1, m2 = 2, 1  # kg
    L1, L2 = 0.5, 0.3  # m
    k1, k2 = 1000, 500  # N/m (bo'g'in elastikligi)
    
    # Sizning kodingiz...
    pass

# robot_arm_vibration()

---

✅ Laboratoriya Tekshirish Ro'yxati

• SHM faza fazosi
• So'nuvchi tebranishlar (3 holat)
• Rezonans simulyatsiyasi
• Bog'langan mayatniklar
• Vibratsiya filtri
• PID tebranish tahlili
• Robot qo'li topshiriq

---

Keyingi Mavzu

📖 2.5 Qattiq Jism Mexanikasi