🔬 Amaliyot: Integral Hisob Simulyatsiyasi

Laboratoriya maqsadi

Bu laboratoriyada integral hisobni Python yordamida amalda qo'llaymiz: raqamli integrallash, yuzalar, hajmlar va fizik masalalar.

Kerakli kutubxonalar

pip install numpy matplotlib scipy

1-laboratoriya: Raqamli Integrallash

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

def numerical_integration_comparison():
    """Turli raqamli integrallash usullarini solishtirish"""
    
    # Test funksiya: sin(x), [0, pi] - aniq integral = 2
    f = np.sin
    a, b = 0, np.pi
    exact = 2.0
    
    def riemann_sum(f, a, b, n, method='midpoint'):
        """Riemann yig'indisi"""
        h = (b - a) / n
        if method == 'left':
            x = np.linspace(a, b - h, n)
        elif method == 'right':
            x = np.linspace(a + h, b, n)
        else:  # midpoint
            x = np.linspace(a + h/2, b - h/2, n)
        return h * np.sum(f(x))
    
    def trapezoidal(f, a, b, n):
        """Trapetsiya qoidasi"""
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        return h * (0.5 * y[0] + np.sum(y[1:-1]) + 0.5 * y[-1])
    
    def simpson(f, a, b, n):
        """Simpson qoidasi (n juft bo'lishi kerak)"""
        if n % 2 == 1:
            n += 1
        x = np.linspace(a, b, n + 1)
        y = f(x)
        h = (b - a) / n
        return h/3 * (y[0] + 4*np.sum(y[1::2]) + 2*np.sum(y[2:-1:2]) + y[-1])
    
    # Turli n qiymatlar uchun xato
    ns = [4, 8, 16, 32, 64, 128, 256, 512]
    
    errors = {
        'Riemann (midpoint)': [],
        'Trapetsiya': [],
        'Simpson': []
    }
    
    for n in ns:
        errors['Riemann (midpoint)'].append(abs(riemann_sum(f, a, b, n) - exact))
        errors['Trapetsiya'].append(abs(trapezoidal(f, a, b, n) - exact))
        errors['Simpson'].append(abs(simpson(f, a, b, n) - exact))
    
    # Vizualizatsiya
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Xato vs n
    ax1 = axes[0]
    for name, err in errors.items():
        ax1.loglog(ns, err, 'o-', label=name)
    
    # Reference lines
    ax1.loglog(ns, 1/np.array(ns)**2, 'k--', alpha=0.5, label='O(1/n²)')
    ax1.loglog(ns, 1/np.array(ns)**4, 'k:', alpha=0.5, label='O(1/n⁴)')
    
    ax1.set_xlabel('n (bo\'linmalar soni)', fontsize=12)
    ax1.set_ylabel('Xato', fontsize=12)
    ax1.set_title('Raqamli integrallash xatolari', fontsize=14)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Trapetsiya vizualizatsiyasi
    ax2 = axes[1]
    n = 8
    x = np.linspace(a, b, 200)
    ax2.plot(x, f(x), 'b-', linewidth=2, label='sin(x)')
    
    x_trap = np.linspace(a, b, n + 1)
    y_trap = f(x_trap)
    
    for i in range(n):
        ax2.fill([x_trap[i], x_trap[i], x_trap[i+1], x_trap[i+1]], 
                [0, y_trap[i], y_trap[i+1], 0], alpha=0.3, color='green')
        ax2.plot([x_trap[i], x_trap[i+1]], [y_trap[i], y_trap[i+1]], 'g-', linewidth=1)
    
    ax2.set_xlabel('x', fontsize=12)
    ax2.set_ylabel('y', fontsize=12)
    ax2.set_title(f'Trapetsiya qoidasi (n={n})', fontsize=14)
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()


def scipy_integration():
    """SciPy bilan integrallash"""
    
    # 1D integral
    f = lambda x: np.exp(-x**2)
    result, error = integrate.quad(f, 0, 2)
    print(f"∫₀² e^(-x²) dx = {result:.6f} ± {error:.2e}")
    
    # 2D integral
    f2 = lambda y, x: np.sin(x) * np.cos(y)
    result2, error2 = integrate.dblquad(f2, 0, np.pi, 0, np.pi/2)
    print(f"∫∫ sin(x)cos(y) dA = {result2:.6f} ± {error2:.2e}")
    
    # ODE yechish (integrallash orqali)
    from scipy.integrate import odeint
    
    # dv/dt = a(t) = 10 - 0.5*v
    def acceleration(v, t):
        return 10 - 0.5 * v
    
    t = np.linspace(0, 20, 200)
    v0 = 0
    v = odeint(acceleration, v0, t)
    
    # Pozitsiya (tezlikni integrallash)
    from scipy.integrate import cumulative_trapezoid
    s = cumulative_trapezoid(v.flatten(), t, initial=0)
    
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    
    axes[0].plot(t, v, 'b-', linewidth=2)
    axes[0].set_xlabel('Vaqt (s)')
    axes[0].set_ylabel('Tezlik (m/s)')
    axes[0].set_title('Tezlik (ODE yechimi)')
    axes[0].grid(True, alpha=0.3)
    
    axes[1].plot(t, s, 'r-', linewidth=2)
    axes[1].set_xlabel('Vaqt (s)')
    axes[1].set_ylabel('Pozitsiya (m)')
    axes[1].set_title('Pozitsiya (tezlik integrali)')
    axes[1].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    numerical_integration_comparison()
    scipy_integration()

2-laboratoriya: Yuzalar Hisoblash

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

def area_between_curves():
    """Egri chiziqlar orasidagi yuza"""
    
    # y = x va y = x²
    f1 = lambda x: x
    f2 = lambda x: x**2
    
    x = np.linspace(-0.5, 1.5, 200)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Grafik
    ax1 = axes[0]
    ax1.plot(x, f1(x), 'b-', linewidth=2, label='y = x')
    ax1.plot(x, f2(x), 'r-', linewidth=2, label='y = x²')
    
    # Yuza
    x_fill = np.linspace(0, 1, 100)
    ax1.fill_between(x_fill, f2(x_fill), f1(x_fill), alpha=0.3, color='green', label='Yuza')
    
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_title('Egri chiziqlar orasidagi yuza')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(-0.5, 1.5)
    ax1.set_ylim(-0.2, 1.2)
    
    # Hisoblash
    area, _ = integrate.quad(lambda x: f1(x) - f2(x), 0, 1)
    ax1.text(0.5, 0.3, f'Yuza = {area:.4f}', fontsize=14, ha='center')
    
    # Sin va cos
    ax2 = axes[1]
    x2 = np.linspace(0, 2*np.pi, 200)
    ax2.plot(x2, np.sin(x2), 'b-', linewidth=2, label='sin(x)')
    ax2.plot(x2, np.cos(x2), 'r-', linewidth=2, label='cos(x)')
    
    # Kesishish nuqtalari: sin(x) = cos(x) -> x = pi/4, 5pi/4
    x_int = [np.pi/4, 5*np.pi/4]
    
    # Yuzalar
    x_fill1 = np.linspace(0, np.pi/4, 50)
    ax2.fill_between(x_fill1, np.sin(x_fill1), np.cos(x_fill1), alpha=0.3, color='orange')
    
    x_fill2 = np.linspace(np.pi/4, 5*np.pi/4, 100)
    ax2.fill_between(x_fill2, np.sin(x_fill2), np.cos(x_fill2), alpha=0.3, color='green')
    
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax2.set_title('sin(x) va cos(x) orasidagi yuza')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # Umumiy yuza
    area1, _ = integrate.quad(lambda x: abs(np.sin(x) - np.cos(x)), 0, 2*np.pi)
    ax2.text(np.pi, 0, f'Umumiy yuza = {area1:.4f}', fontsize=12, ha='center')
    
    plt.tight_layout()
    plt.show()


def parametric_area():
    """Parametrik egri chiziq yuzasi"""
    
    # Ellips: x = a*cos(t), y = b*sin(t)
    a, b = 3, 2
    t = np.linspace(0, 2*np.pi, 200)
    
    x = a * np.cos(t)
    y = b * np.sin(t)
    
    # Yuza = (1/2) * |∮ (x dy - y dx)|
    # = (1/2) * |∫ x(t) y'(t) - y(t) x'(t) dt|
    
    def area_parametric():
        dx_dt = -a * np.sin(t)
        dy_dt = b * np.cos(t)
        integrand = x * dy_dt - y * dx_dt
        return 0.5 * np.abs(np.trapz(integrand, t))
    
    area_num = area_parametric()
    area_exact = np.pi * a * b
    
    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(x, y, 'b-', linewidth=2)
    ax.fill(x, y, alpha=0.3, color='blue')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(f'Ellips yuzasi: a={a}, b={b}')
    ax.set_aspect('equal')
    ax.grid(True, alpha=0.3)
    
    ax.text(0, 0, f'Raqamli: {area_num:.4f}\nAniq: πab = {area_exact:.4f}', 
            fontsize=12, ha='center', va='center')
    
    plt.show()


if __name__ == "__main__":
    area_between_curves()
    parametric_area()

3-laboratoriya: Hajmlar va Aylanish Jismlari

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy import integrate

def volume_of_revolution():
    """Aylanish jismi hajmi"""
    
    # y = sqrt(x), [0, 4] ni x-o'qi atrofida aylantirish
    f = lambda x: np.sqrt(x)
    a, b = 0, 4
    
    # Disk usuli: V = π ∫ [f(x)]² dx
    volume, _ = integrate.quad(lambda x: np.pi * f(x)**2, a, b)
    
    # 3D vizualizatsiya
    fig = plt.figure(figsize=(14, 5))
    
    # 2D profil
    ax1 = fig.add_subplot(131)
    x = np.linspace(a, b, 100)
    ax1.plot(x, f(x), 'b-', linewidth=2, label='y = √x')
    ax1.plot(x, -f(x), 'b-', linewidth=2)
    ax1.fill_between(x, -f(x), f(x), alpha=0.3, color='blue')
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_title('Profil')
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)
    
    # 3D surface
    ax2 = fig.add_subplot(132, projection='3d')
    
    u = np.linspace(a, b, 50)
    v = np.linspace(0, 2*np.pi, 50)
    U, V = np.meshgrid(u, v)
    
    X = U
    Y = f(U) * np.cos(V)
    Z = f(U) * np.sin(V)
    
    ax2.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)
    ax2.set_xlabel('X')
    ax2.set_ylabel('Y')
    ax2.set_zlabel('Z')
    ax2.set_title(f'Aylanish jismi\nV = {volume:.4f}')
    
    # Disk kesimlar
    ax3 = fig.add_subplot(133)
    
    # Har bir x uchun disk
    for xi in [0.5, 1, 2, 3]:
        r = f(xi)
        circle = plt.Circle((0, 0), r, fill=False, color='blue', linewidth=2)
        ax3.add_patch(circle)
        ax3.text(0, r + 0.1, f'x={xi}', ha='center', fontsize=10)
    
    ax3.set_xlim(-2.5, 2.5)
    ax3.set_ylim(-2.5, 2.5)
    ax3.set_aspect('equal')
    ax3.set_title('Disk kesimlari (YZ tekislik)')
    ax3.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    print(f"Hajm (aniq): V = 8π = {8*np.pi:.4f}")
    print(f"Hajm (raqamli): V = {volume:.4f}")


def shell_method():
    """Silindr qobig'i usuli"""
    
    # y = x² ni y-o'qi atrofida aylantirish, [0, 2]
    # Shell: V = 2π ∫ x * f(x) dx
    
    f = lambda x: x**2
    a, b = 0, 2
    
    volume, _ = integrate.quad(lambda x: 2*np.pi * x * f(x), a, b)
    
    fig = plt.figure(figsize=(12, 5))
    
    # 2D
    ax1 = fig.add_subplot(121)
    x = np.linspace(a, b, 100)
    ax1.plot(x, f(x), 'b-', linewidth=2, label='y = x²')
    ax1.fill_between(x, 0, f(x), alpha=0.3, color='blue')
    ax1.axvline(x=0, color='k', linewidth=1)
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_title('Profil')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # 3D
    ax2 = fig.add_subplot(122, projection='3d')
    
    u = np.linspace(a, b, 30)
    v = np.linspace(0, 2*np.pi, 50)
    U, V = np.meshgrid(u, v)
    
    # Paraboloid
    X = U * np.cos(V)
    Y = U * np.sin(V)
    Z = f(U)
    
    ax2.plot_surface(X, Y, Z, cmap='plasma', alpha=0.7)
    ax2.set_xlabel('X')
    ax2.set_ylabel('Y')
    ax2.set_zlabel('Z')
    ax2.set_title(f'y = x² aylanishi\nV = {volume:.4f}')
    
    plt.tight_layout()
    plt.show()
    
    print(f"Hajm = 2π ∫₀² x·x² dx = 2π · x⁴/4 |₀² = 8π = {8*np.pi:.4f}")


if __name__ == "__main__":
    volume_of_revolution()
    shell_method()

4-laboratoriya: Fizik Qo'llanmalar

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import cumulative_trapezoid, quad

def rocket_motion():
    """Raketa harakati simulyatsiyasi"""
    
    # Tezlanish modeli: a(t) = a0 * exp(-t/tau)
    a0 = 30  # m/s² (boshlang'ich)
    tau = 20  # s (vaqt konstantasi)
    
    t = np.linspace(0, 60, 1000)
    
    # Tezlanish
    a = a0 * np.exp(-t / tau)
    
    # Tezlik (integrallash)
    v = cumulative_trapezoid(a, t, initial=0)
    
    # Pozitsiya (ikki marta integrallash)
    s = cumulative_trapezoid(v, t, initial=0)
    
    # Vizualizatsiya
    fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True)
    
    axes[0].plot(t, a, 'r-', linewidth=2)
    axes[0].set_ylabel('Tezlanish (m/s²)', fontsize=12)
    axes[0].set_title('Raketa Harakati', fontsize=14)
    axes[0].grid(True, alpha=0.3)
    axes[0].fill_between(t, 0, a, alpha=0.3, color='red')
    
    axes[1].plot(t, v, 'g-', linewidth=2)
    axes[1].set_ylabel('Tezlik (m/s)', fontsize=12)
    axes[1].grid(True, alpha=0.3)
    
    axes[2].plot(t, s/1000, 'b-', linewidth=2)
    axes[2].set_xlabel('Vaqt (s)', fontsize=12)
    axes[2].set_ylabel('Masofa (km)', fontsize=12)
    axes[2].grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Aniq hisoblash
    # v = ∫ a dt = a0*tau * (1 - exp(-t/tau))
    v_final_exact = a0 * tau * (1 - np.exp(-60/tau))
    print(f"Final tezlik (raqamli): {v[-1]:.2f} m/s")
    print(f"Final tezlik (aniq): {v_final_exact:.2f} m/s")


def work_calculation():
    """Ish hisoblash"""
    
    # O'zgaruvchan kuch: F(x) = 10 + 5*sin(2*pi*x)
    F = lambda x: 10 + 5 * np.sin(2 * np.pi * x)
    
    x = np.linspace(0, 3, 200)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Kuch grafigi
    ax1 = axes[0]
    ax1.plot(x, F(x), 'b-', linewidth=2)
    ax1.fill_between(x, 0, F(x), alpha=0.3, color='blue')
    ax1.set_xlabel('Masofa x (m)', fontsize=12)
    ax1.set_ylabel('Kuch F (N)', fontsize=12)
    ax1.set_title('O\'zgaruvchan kuch', fontsize=14)
    ax1.grid(True, alpha=0.3)
    
    # Kumulyativ ish
    work_cumulative = cumulative_trapezoid(F(x), x, initial=0)
    
    ax2 = axes[1]
    ax2.plot(x, work_cumulative, 'r-', linewidth=2)
    ax2.set_xlabel('Masofa x (m)', fontsize=12)
    ax2.set_ylabel('Ish W (J)', fontsize=12)
    ax2.set_title('Kumulyativ ish', fontsize=14)
    ax2.grid(True, alpha=0.3)
    
    # Umumiy ish
    total_work, _ = quad(F, 0, 3)
    ax2.axhline(y=total_work, color='g', linestyle='--', label=f'Total: {total_work:.2f} J')
    ax2.legend()
    
    plt.tight_layout()
    plt.show()


def pid_integral_term():
    """PID kontroller integral termi"""
    
    t = np.linspace(0, 10, 500)
    dt = t[1] - t[0]
    
    # Xato signali (exponensial pasayish)
    e = 5 * np.exp(-0.3 * t) * np.sin(2 * t)
    
    # Integral term
    Ki = 2.0
    integral_term = cumulative_trapezoid(e, t, initial=0) * Ki
    
    # Vizualizatsiya
    fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
    
    axes[0].plot(t, e, 'b-', linewidth=2, label='Xato e(t)')
    axes[0].fill_between(t, 0, e, alpha=0.3, where=(e >= 0), color='green', label='Ijobiy')
    axes[0].fill_between(t, 0, e, alpha=0.3, where=(e < 0), color='red', label='Manfiy')
    axes[0].set_ylabel('Xato', fontsize=12)
    axes[0].set_title('PID Kontroller Integral Termi', fontsize=14)
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    axes[1].plot(t, integral_term, 'r-', linewidth=2)
    axes[1].set_xlabel('Vaqt (s)', fontsize=12)
    axes[1].set_ylabel(f'Integral term (Ki={Ki})', fontsize=12)
    axes[1].grid(True, alpha=0.3)
    axes[1].axhline(y=0, color='k', linestyle='-', alpha=0.3)
    
    plt.tight_layout()
    plt.show()


def energy_consumption():
    """Energiya sarfi hisoblash"""
    
    t = np.linspace(0, 100, 1000)
    
    # Quvvat modeli (dron hovering + manevrlar)
    P_base = 200  # W
    P = P_base + 50 * np.sin(0.2 * t) + 30 * np.sin(0.5 * t)
    
    # Energiya (kumulyativ)
    E = cumulative_trapezoid(P, t, initial=0)
    
    fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
    
    axes[0].plot(t, P, 'b-', linewidth=1)
    axes[0].axhline(y=P_base, color='r', linestyle='--', label=f'Bazaviy: {P_base} W')
    axes[0].set_ylabel('Quvvat (W)', fontsize=12)
    axes[0].set_title('Dron Energiya Sarfi', fontsize=14)
    axes[0].legend()
    axes[0].grid(True, alpha=0.3)
    
    axes[1].plot(t, E / 1000, 'g-', linewidth=2)
    axes[1].set_xlabel('Vaqt (s)', fontsize=12)
    axes[1].set_ylabel('Energiya (kJ)', fontsize=12)
    axes[1].grid(True, alpha=0.3)
    
    # Batareya qoldig'i (15 kJ batareya)
    battery_capacity = 15  # kJ
    remaining = battery_capacity - E / 1000
    
    ax3 = axes[1].twinx()
    ax3.plot(t, remaining / battery_capacity * 100, 'r--', linewidth=1, alpha=0.7)
    ax3.set_ylabel('Batareya (%)', color='red', fontsize=12)
    ax3.tick_params(axis='y', labelcolor='red')
    
    plt.tight_layout()
    plt.show()
    
    # Flight time estimation
    total_energy, _ = quad(lambda t: 200 + 50*np.sin(0.2*t) + 30*np.sin(0.5*t), 0, 100)
    flight_time = battery_capacity * 1000 / (total_energy / 100)
    print(f"Taxminiy parvoz vaqti: {flight_time:.1f} s ({flight_time/60:.1f} min)")


if __name__ == "__main__":
    rocket_motion()
    work_calculation()
    pid_integral_term()
    energy_consumption()

5-laboratoriya: Monte Carlo Integrallash

import numpy as np
import matplotlib.pyplot as plt

def monte_carlo_integration():
    """Monte Carlo usuli bilan integrallash"""
    
    # ∫∫ f(x,y) dA, bu yerda f(x,y) = x² + y², [0,1] x [0,1]
    f = lambda x, y: x**2 + y**2
    
    # Aniq qiymat: ∫₀¹∫₀¹ (x² + y²) dx dy = 2/3
    exact = 2/3
    
    # Monte Carlo
    def monte_carlo(n):
        x = np.random.uniform(0, 1, n)
        y = np.random.uniform(0, 1, n)
        return np.mean(f(x, y))  # Area = 1
    
    # Turli n uchun
    ns = [100, 500, 1000, 5000, 10000, 50000, 100000]
    estimates = []
    errors = []
    
    for n in ns:
        trials = [monte_carlo(n) for _ in range(50)]
        estimates.append(np.mean(trials))
        errors.append(np.std(trials))
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Xato vs n
    ax1 = axes[0]
    ax1.errorbar(ns, estimates, yerr=errors, fmt='bo-', capsize=3)
    ax1.axhline(y=exact, color='r', linestyle='--', label=f'Aniq = {exact:.4f}')
    ax1.set_xscale('log')
    ax1.set_xlabel('n (nuqtalar soni)', fontsize=12)
    ax1.set_ylabel('Taxmin', fontsize=12)
    ax1.set_title('Monte Carlo Integrallash', fontsize=14)
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Vizualizatsiya
    ax2 = axes[1]
    n_vis = 1000
    x = np.random.uniform(0, 1, n_vis)
    y = np.random.uniform(0, 1, n_vis)
    z = f(x, y)
    
    scatter = ax2.scatter(x, y, c=z, cmap='viridis', s=10, alpha=0.6)
    plt.colorbar(scatter, ax=ax2, label='f(x,y) = x² + y²')
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax2.set_title(f'Monte Carlo nuqtalari (n={n_vis})')
    ax2.set_aspect('equal')
    
    plt.tight_layout()
    plt.show()


def pi_estimation():
    """Monte Carlo bilan π ni hisoblash"""
    
    # Birlik doira ichidagi nuqtalar / umumiy nuqtalar = π/4
    
    def estimate_pi(n):
        x = np.random.uniform(-1, 1, n)
        y = np.random.uniform(-1, 1, n)
        inside = (x**2 + y**2) <= 1
        return 4 * np.sum(inside) / n
    
    ns = np.logspace(2, 6, 30).astype(int)
    pi_estimates = [estimate_pi(n) for n in ns]
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))
    
    # Taxminlar
    ax1 = axes[0]
    ax1.semilogx(ns, pi_estimates, 'bo-', markersize=4)
    ax1.axhline(y=np.pi, color='r', linestyle='--', label=f'π = {np.pi:.6f}')
    ax1.set_xlabel('n')
    ax1.set_ylabel('π taxmini')
    ax1.set_title('Monte Carlo bilan π hisoblash')
    ax1.legend()
    ax1.grid(True, alpha=0.3)
    
    # Vizualizatsiya
    ax2 = axes[1]
    n_vis = 5000
    x = np.random.uniform(-1, 1, n_vis)
    y = np.random.uniform(-1, 1, n_vis)
    inside = (x**2 + y**2) <= 1
    
    ax2.scatter(x[inside], y[inside], c='blue', s=1, alpha=0.5, label='Ichida')
    ax2.scatter(x[~inside], y[~inside], c='red', s=1, alpha=0.5, label='Tashqarida')
    
    theta = np.linspace(0, 2*np.pi, 100)
    ax2.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=2)
    
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax2.set_title(f'π ≈ 4 × {np.sum(inside)}/{n_vis} = {4*np.sum(inside)/n_vis:.4f}')
    ax2.set_aspect('equal')
    ax2.legend(loc='upper right')
    
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    monte_carlo_integration()
    pi_estimation()

Xulosa

Bu laboratoriyalarda biz:

• ✅ Raqamli integrallash usullarini solishtirdik
• ✅ Yuzalar va hajmlarni hisobladik
• ✅ Fizik masalalarni yechdik (raketa, ish, energiya)
• ✅ Monte Carlo integrallashni o'rgandik

Keyingi qadamlar

1. Differensial tenglamalar
2. Ko'p o'lchovli integrallar
3. Chiziqli algebra qo'llanmalari