Amaliyot — Differensial Tenglamalar

Python va SciPy yordamida ODE yechish va simulyatsiya qilish.

---

Loyiha: Tebranish Sistemalarini Modellashtirish

Maqsad

So'nuvchi tebranishlarni simulyatsiya qilish va turli so'nish koeffitsientlari ($\zeta$) ta'sirini vizualizatsiya qilish.

---

1-Qism: Muhitni Sozlash

# Kerakli kutubxonalar
import numpy as np
from scipy.integrate import odeint, solve_ivp
import matplotlib.pyplot as plt

# Grafika sozlamalari
plt.rcParams['figure.figsize'] = (12, 6)
plt.rcParams['font.size'] = 12
plt.rcParams['axes.grid'] = True

---

2-Qism: So'nuvchi Ossilyator

Matematik Model

$$\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 x = 0$$

State-Space Ko'rinishi

def damped_oscillator(t, y, zeta, omega_n):
    """
    So'nuvchi ossilyator state-space modeli.
    
    y = [x, v] - holat vektori
    zeta - so'nish koeffitsienti
    omega_n - tabiiy chastota
    """
    x, v = y
    dxdt = v
    dvdt = -2*zeta*omega_n*v - omega_n**2*x
    return [dxdt, dvdt]

Simulyatsiya

# Parametrlar
omega_n = 2*np.pi  # 1 Hz tabiiy chastota
t_span = (0, 5)
t_eval = np.linspace(0, 5, 500)
y0 = [1.0, 0.0]  # x(0)=1, v(0)=0

# Turli so'nish koeffitsientlari
zeta_values = [0.0, 0.1, 0.3, 0.7, 1.0, 2.0]
labels = ['So\'nishsiz (ζ=0)', 
          'Kam so\'nuvchi (ζ=0.1)',
          'Kam so\'nuvchi (ζ=0.3)', 
          'Kam so\'nuvchi (ζ=0.7)',
          'Kritik (ζ=1.0)', 
          'O\'ta so\'nuvchi (ζ=2.0)']

plt.figure(figsize=(14, 8))
for zeta, label in zip(zeta_values, labels):
    sol = solve_ivp(damped_oscillator, t_span, y0, 
                    args=(zeta, omega_n), t_eval=t_eval)
    plt.plot(sol.t, sol.y[0], label=label, linewidth=2)

plt.xlabel('Vaqt (s)')
plt.ylabel('Siljish x(t)')
plt.title('So\'nish Koeffitsientining Ta\'siri')
plt.legend()
plt.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
plt.tight_layout()
plt.savefig('damping_comparison.png', dpi=150)
plt.show()

Kutilgan Natija

• $\zeta = 0$: Cheksiz tebranadi
• $\zeta < 1$: Tebranib so'nadi
• $\zeta = 1$: Eng tez muvozanatga qaytadi (tebranishsiz)
• $\zeta > 1$: Sekin muvozanatga qaytadi

---

3-Qism: Majburiy Tebranishlar va Rezonans

Matematik Model

$$\ddot{x} + 2\zeta\omega_n\dot{x} + \omega_n^2 x = F_0\cos(\omega t)$$

def forced_oscillator(t, y, zeta, omega_n, F0, omega):
    """Majburiy tebranish modeli."""
    x, v = y
    dxdt = v
    dvdt = -2*zeta*omega_n*v - omega_n**2*x + F0*np.cos(omega*t)
    return [dxdt, dvdt]

# Parametrlar
omega_n = 2*np.pi
zeta = 0.1
F0 = 1.0
y0 = [0.0, 0.0]
t_span = (0, 20)
t_eval = np.linspace(0, 20, 2000)

# Turli majburiy chastotalar
omega_values = [0.5*omega_n, 0.9*omega_n, omega_n, 1.1*omega_n, 2*omega_n]

fig, axs = plt.subplots(len(omega_values), 1, figsize=(14, 12))
for i, omega in enumerate(omega_values):
    sol = solve_ivp(forced_oscillator, t_span, y0,
                    args=(zeta, omega_n, F0, omega), t_eval=t_eval)
    axs[i].plot(sol.t, sol.y[0], 'b-', linewidth=1.5)
    axs[i].set_ylabel(f'ω = {omega/omega_n:.1f}ωn')
    axs[i].set_xlim([0, 20])

axs[-1].set_xlabel('Vaqt (s)')
fig.suptitle('Majburiy Tebranishlar (ζ = 0.1)', fontsize=14)
plt.tight_layout()
plt.savefig('forced_oscillation.png', dpi=150)
plt.show()

Rezonans Egri Chizig'i

def amplitude_response(omega, omega_n, zeta, F0=1.0):
    """Amplituda-chastota xarakteristikasi."""
    r = omega / omega_n
    A = F0 / (omega_n**2 * np.sqrt((1 - r**2)**2 + (2*zeta*r)**2))
    return A

omega_range = np.linspace(0.01, 3, 500) * omega_n
zeta_values = [0.05, 0.1, 0.3, 0.5, 1.0]

plt.figure(figsize=(12, 6))
for zeta in zeta_values:
    A = amplitude_response(omega_range, omega_n, zeta)
    plt.plot(omega_range/omega_n, A, label=f'ζ = {zeta}', linewidth=2)

plt.xlabel('ω / ωn')
plt.ylabel('Amplituda')
plt.title('Rezonans Egri Chiziqlari')
plt.legend()
plt.axvline(x=1, color='r', linestyle='--', alpha=0.5)
plt.xlim([0, 3])
plt.ylim([0, 12])
plt.savefig('resonance_curves.png', dpi=150)
plt.show()

---

4-Qism: Raqamli Usullarni Taqqoslash

Euler vs RK4

def euler_method(f, y0, t_span, h):
    """Oddiy Euler usuli."""
    t = np.arange(t_span[0], t_span[1] + h, h)
    y = np.zeros((len(t), len(y0)))
    y[0] = y0
    
    for i in range(len(t) - 1):
        y[i+1] = y[i] + h * np.array(f(t[i], y[i]))
    
    return t, y

def rk4_method(f, y0, t_span, h):
    """Runge-Kutta 4-tartib usuli."""
    t = np.arange(t_span[0], t_span[1] + h, h)
    y = np.zeros((len(t), len(y0)))
    y[0] = y0
    
    for i in range(len(t) - 1):
        k1 = np.array(f(t[i], y[i]))
        k2 = np.array(f(t[i] + h/2, y[i] + h*k1/2))
        k3 = np.array(f(t[i] + h/2, y[i] + h*k2/2))
        k4 = np.array(f(t[i] + h, y[i] + h*k3))
        y[i+1] = y[i] + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
    
    return t, y

# Test funksiyasi: y' = -y, y(0) = 1, aniq yechim: y = e^(-t)
def simple_ode(t, y):
    return [-y[0]]

y0 = [1.0]
t_span = (0, 5)
h = 0.5  # Katta qadam (xatolikni ko'rish uchun)

# Usullarni qo'llash
t_euler, y_euler = euler_method(simple_ode, y0, t_span, h)
t_rk4, y_rk4 = rk4_method(simple_ode, y0, t_span, h)
t_exact = np.linspace(0, 5, 100)
y_exact = np.exp(-t_exact)

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t_exact, y_exact, 'k-', label='Aniq yechim', linewidth=2)
plt.plot(t_euler, y_euler[:, 0], 'ro-', label='Euler', markersize=8)
plt.plot(t_rk4, y_rk4[:, 0], 'bs-', label='RK4', markersize=8)
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title(f'Yechimlarni Taqqoslash (h = {h})')

plt.subplot(1, 2, 2)
error_euler = np.abs(y_euler[:, 0] - np.exp(-t_euler))
error_rk4 = np.abs(y_rk4[:, 0] - np.exp(-t_rk4))
plt.semilogy(t_euler, error_euler, 'ro-', label='Euler xatoligi')
plt.semilogy(t_rk4, error_rk4, 'bs-', label='RK4 xatoligi')
plt.xlabel('t')
plt.ylabel('Absolyut xatolik')
plt.legend()
plt.title('Xatoliklarni Taqqoslash')

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

---

5-Qism: Raketa Simulyatsiyasi

Tsiolkovsky Tenglamasi

def rocket_1d(t, y, ve, mdot_0, m_empty, g):
    """
    1D raketa modeli (vertikal uchish).
    
    y = [h, v, m] - balandlik, tezlik, massa
    ve - gazlar chiqish tezligi
    mdot_0 - yonilg'i sarfi (kg/s)
    m_empty - quruq massa
    g - gravitatsiya
    """
    h, v, m = y
    
    # Yonilg'i tugadimi?
    if m <= m_empty:
        mdot = 0
        thrust = 0
    else:
        mdot = mdot_0
        thrust = ve * mdot
    
    # Havo qarshiligi (soddalashtirilgan)
    rho = 1.225 * np.exp(-h / 8500)  # Atmosfera modeli
    Cd = 0.3
    A = 0.5  # m^2
    drag = 0.5 * rho * v**2 * Cd * A * np.sign(v)
    
    dhdt = v
    dvdt = thrust / max(m, m_empty) - g - drag / max(m, m_empty)
    dmdt = -mdot
    
    return [dhdt, dvdt, dmdt]

# Parametrlar (model raketa)
ve = 2500        # m/s (qattiq yonilg'i)
mdot_0 = 5       # kg/s
m_total = 100    # kg
m_empty = 20     # kg
g = 9.81

y0 = [0, 0, m_total]
t_span = (0, 60)
t_eval = np.linspace(0, 60, 600)

sol = solve_ivp(rocket_1d, t_span, y0, 
                args=(ve, mdot_0, m_empty, g),
                t_eval=t_eval, max_step=0.1)

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

# Balandlik
axs[0, 0].plot(sol.t, sol.y[0] / 1000, 'b-', linewidth=2)
axs[0, 0].set_ylabel('Balandlik (km)')
axs[0, 0].set_title('Balandlik vs Vaqt')
axs[0, 0].axhline(y=0, color='k', linewidth=0.5)

# Tezlik
axs[0, 1].plot(sol.t, sol.y[1], 'r-', linewidth=2)
axs[0, 1].set_ylabel('Tezlik (m/s)')
axs[0, 1].set_title('Tezlik vs Vaqt')
axs[0, 1].axhline(y=0, color='k', linewidth=0.5)

# Massa
axs[1, 0].plot(sol.t, sol.y[2], 'g-', linewidth=2)
axs[1, 0].set_xlabel('Vaqt (s)')
axs[1, 0].set_ylabel('Massa (kg)')
axs[1, 0].set_title('Massa vs Vaqt')
axs[1, 0].axhline(y=m_empty, color='r', linestyle='--', label='Quruq massa')
axs[1, 0].legend()

# Tezlanish
a = np.gradient(sol.y[1], sol.t) / g
axs[1, 1].plot(sol.t, a, 'm-', linewidth=2)
axs[1, 1].set_xlabel('Vaqt (s)')
axs[1, 1].set_ylabel('Tezlanish (g)')
axs[1, 1].set_title('Tezlanish vs Vaqt')
axs[1, 1].axhline(y=0, color='k', linewidth=0.5)

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

# Statistika
max_h = np.max(sol.y[0])
max_v = np.max(sol.y[1])
burn_time = (m_total - m_empty) / mdot_0

print(f"Maksimal balandlik: {max_h/1000:.2f} km")
print(f"Maksimal tezlik: {max_v:.1f} m/s")
print(f"Yonish vaqti: {burn_time:.1f} s")
print(f"Delta-v (nazariy): {ve * np.log(m_total/m_empty):.1f} m/s")

---

6-Qism: Dron Barqarorlik Tahlili

Quadcopter Pitch Dinamikasi

def quadcopter_pitch(t, y, Iyy, l, k_f, omega_sq_diff, k_d):
    """
    Quadcopter pitch (tanlanish) dinamikasi.
    
    y = [theta, omega] - burchak va burchak tezligi
    Iyy - inertsiya momenti
    l - motor-markaz masofasi
    k_f - itarish koeffitsienti
    omega_sq_diff - old va orqa motorlar omega^2 farqi
    k_d - aerodinamik so'nish
    """
    theta, omega = y
    
    torque = l * k_f * omega_sq_diff
    
    dthetadt = omega
    domegadt = (torque - k_d * omega) / Iyy
    
    return [dthetadt, domegadt]

# Parametrlar
Iyy = 0.02      # kg*m^2
l = 0.2         # m
k_f = 1e-5      # N/(rad/s)^2
k_d = 0.001     # aerodinamik so'nish

# Boshlang'ich burchak og'ishi
y0 = [np.radians(10), 0]  # 10 daraja
t_span = (0, 5)
t_eval = np.linspace(0, 5, 500)

# Boshqaruvchi kiritmagan holda (omega_sq_diff = 0)
sol_no_control = solve_ivp(quadcopter_pitch, t_span, y0,
                           args=(Iyy, l, k_f, 0, k_d), t_eval=t_eval)

# Oddiy P-kontrollerli (omega_sq_diff ~ -theta)
def quadcopter_with_control(t, y, Iyy, l, k_f, k_d, Kp):
    theta, omega = y
    omega_sq_diff = -Kp * theta  # P-kontroller
    torque = l * k_f * omega_sq_diff
    return [omega, (torque - k_d * omega) / Iyy]

Kp = 5000
sol_with_control = solve_ivp(quadcopter_with_control, t_span, y0,
                             args=(Iyy, l, k_f, k_d, Kp), t_eval=t_eval)

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(sol_no_control.t, np.degrees(sol_no_control.y[0]), 'r-', 
         label='Boshqaruvsiz', linewidth=2)
plt.plot(sol_with_control.t, np.degrees(sol_with_control.y[0]), 'b-',
         label=f'P-kontroller (Kp={Kp})', linewidth=2)
plt.xlabel('Vaqt (s)')
plt.ylabel('Pitch burchagi (°)')
plt.legend()
plt.title('Dron Pitch Barqarorligi')
plt.axhline(y=0, color='k', linewidth=0.5)

plt.subplot(1, 2, 2)
plt.plot(sol_no_control.t, np.degrees(sol_no_control.y[1]), 'r-',
         label='Boshqaruvsiz', linewidth=2)
plt.plot(sol_with_control.t, np.degrees(sol_with_control.y[1]), 'b-',
         label=f'P-kontroller', linewidth=2)
plt.xlabel('Vaqt (s)')
plt.ylabel('Burchak tezligi (°/s)')
plt.legend()
plt.title('Burchak Tezligi')

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

---

7-Qism: Topshiriqlar

Topshiriq 1: So'nish Tahlili

1. $\zeta = 0.2$ va $\omega_n = 5\,rad/s$ uchun so'nuvchi ossilyatorni simulyatsiya qiling
2. Amplituda 10% ga tushguncha nechta tebranish kerakligini hisoblang
3. Nazariy va simulyatsiya natijalarini taqqoslang

Topshiriq 2: Raketa Optimallashtirish

1. Yuqoridagi raketa simulyatsiyasida yonilg'i sarfini ($\dot{m}$) o'zgartirib, maksimal balandlikni tadqiq qiling
2. Optimal $\dot{m}$ qiymatini toping
3. Nima uchun juda katta va juda kichik $\dot{m}$ samarasiz?

Topshiriq 3: Dron PID Tuning

1. Quadcopter modeliga I va D komponentlarini qo'shing
2. Ziegler-Nichols usuli bilan PID parametrlarini tanlang
3. P, PI, va PID kontrollerlarni taqqoslang

---

Xulosa

Bu amaliyotda:

• So'nuvchi tebranishlarni modellashtirishni o'rgandik
• Euler va RK4 usullarini qo'lladik
• Raketa va dron dinamikasini simulyatsiya qildik
• Boshqarish tizimlarining asoslarini ko'rdik

Keyingi qadam: Boshqarish nazariyasi (Control Theory) bo'limida ushbu bilimlarni chuqurlashtirish.