TechOverflow Logo
TechOverflow

Затухаючий пружинний осцилятор: аналітичні та чисельні рішення на Python

English Deutsch Français 简体中文 Español Čeština Português Русский Українська

Ця публікація демонструє два підходи до симуляції затухаючого пружинного осцилятора: аналітичне та чисельне рішення.

Методи симуляції

Аналітичне рішення

Затухаючий гармонічний осцилятор описується диференціальним рівнянням:

$$\frac{d^2x}{dt^2} + 2\zeta\omega_0\frac{dx}{dt} + \omega_0^2x = 0$$

Для недемпфованої системи ($\zeta < 1$) аналітичне рішення має вигляд:

$$x(t) = e^{-\zeta\omega_0 t}\left[x_0\cos(\omega_d t) + \frac{v_0 + \zeta\omega_0 x_0}{\omega_d}\sin(\omega_d t)\right]$$

Де:

Чисельне рішення

Чисельний підхід використовує метод Рунге-Кутти 4-го порядку для інтегрування системи диференціальних рівнянь першого порядку:

$$\frac{dx}{dt} = v$$

$$\frac{dv}{dt} = -\omega_0^2x - 2\zeta\omega_0v$$

Цей метод забезпечує високу точність і дозволяє порівняти з аналітичним рішенням для перевірки коректності.

Затухаючий пружинний осцилятор

Вихідний код

damped_spring_oscillator.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# SPDX-License-Identifier: CC0 1.0 Universal

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")

# Параметри симуляції
fs = 1000  # Частота дискретизації (Гц)
t_duration = 5  # Тривалість (секунди)
f_resonant = 4  # Резонансна частота (Гц)
x0 = 1.0  # Початкове зміщення
v0 = 0.0  # Початкова швидкість

# Похідні параметри
omega0 = 2 * np.pi * f_resonant  # Власна кутова частота
dt = 1.0 / fs  # Крок часу
t = np.arange(0, t_duration, dt)  # Масив часу
n_samples = len(t)

# Коефіцієнт затухання (можна змінювати для керування затуханням)
# zeta < 1: недемпфований, zeta = 1: критично демпфований, zeta > 1: передемпфований
zeta = 0.03  # Легке затухання для коливального руху
gamma = 2 * zeta * omega0  # Коефіцієнт затухання

# Аналітичне рішення для недемпфованого осцилятора
omega_d = omega0 * np.sqrt(1 - zeta**2)  # Демпфована частота

# Для початкових умов x(0) = x0, v(0) = v0
# Рішення: x(t) = e^(-ζω₀t)[x0*cos(ωd*t) + ((v0 + ζω₀x0)/ωd)*sin(ωd*t)]
x_analytical = np.exp(-zeta * omega0 * t) * (
    x0 * np.cos(omega_d * t) +
    ((v0 + zeta * omega0 * x0) / omega_d) * np.sin(omega_d * t)
)

print(f"Власна частота: {f_resonant:.2f} Гц")
print(f"Демпфована частота: {omega_d / (2 * np.pi):.2f} Гц")
print(f"Коефіцієнт затухання: {zeta:.2f}")
print(f"Коефіцієнт затухання gamma: {gamma:.3f}")
print(f"Власна кутова частота ω₀: {omega0:.3f} рад/с")

# Чисельна симуляція з використанням методу Рунге-Кутти
def spring_ode(state, t):
    """
    Диференціальне рівняння для затухаючого гармонічного осцилятора
    d²x/dt² + 2ζω₀(dx/dt) + ω₀²x = 0
    state = [положення, швидкість]
    повертає [швидкість, прискорення]
    """
    x, v = state
    # прискорення = -ω₀²x - 2ζω₀v
    acceleration = -omega0**2 * x - 2 * zeta * omega0 * v
    return np.array([v, acceleration])

# Ініціалізація масивів для чисельного рішення
x_numerical = np.zeros(n_samples)
v_numerical = np.zeros(n_samples)

# Початкові умови
x_numerical[0] = x0
v_numerical[0] = v0

# Інтегрування методом Рунге-Кутти 4-го порядку
for i in range(n_samples - 1):
    state = np.array([x_numerical[i], v_numerical[i]])

    k1 = dt * spring_ode(state, t[i])
    k2 = dt * spring_ode(state + k1/2, t[i] + dt/2)
    k3 = dt * spring_ode(state + k2/2, t[i] + dt/2)
    k4 = dt * spring_ode(state + k3, t[i] + dt)

    next_state = state + (k1 + 2*k2 + 2*k3 + k4) / 6
    x_numerical[i+1] = next_state[0]
    v_numerical[i+1] = next_state[1]

# Побудова графіків результатів
plt.figure(figsize=(12, 8))

# Графік положення
plt.subplot(2, 1, 1)
plt.plot(t, x_analytical, 'b-', label='Аналітичне рішення', linewidth=2)
plt.plot(t[::10], x_numerical[::10], 'ro', markersize=2, label='Чисельне рішення', alpha=0.7)
plt.xlabel('Час (с)')
plt.ylabel('Положення (м)')
plt.title(f'Затухаючий пружинний осцилятор (f₀={f_resonant}Гц, ζ={zeta}, fs={fs}Гц)')
plt.legend()
plt.grid(True, alpha=0.3)

# Аналітична швидкість (похідна положення)
v_analytical = np.exp(-zeta * omega0 * t) * (
    (-zeta * omega0 * x0 + (v0 + zeta * omega0 * x0) * omega_d / omega_d) * np.cos(omega_d * t) +
    (-omega_d * x0 - zeta * omega0 * (v0 + zeta * omega0 * x0) / omega_d) * np.sin(omega_d * t)
)

# Графік швидкості
plt.subplot(2, 1, 2)
plt.plot(t, v_analytical, 'g-', label='Аналітична швидкість', linewidth=2)
plt.plot(t[::10], v_numerical[::10], 'mo', markersize=2, label='Чисельна швидкість', alpha=0.7)
plt.xlabel('Час (с)')
plt.ylabel('Швидкість (м/с)')
plt.title('Швидкість пружинного маса')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Обчислення похибки між аналітичним та чисельним рішеннями
error = np.abs(x_analytical - x_numerical)
max_error = np.max(error)
rms_error = np.sqrt(np.mean(error**2))

print(f"\nЧисельна точність:")
print(f"Максимальна похибка: {max_error:.6f} м")
print(f"Середньоквадратична похибка: {rms_error:.6f} м")
print(f"Відносна середньоквадратична похибка: {rms_error/np.max(np.abs(x_analytical))*100:.4f}%")
plt.savefig("dampened_spring_oscillator.svg", dpi=300, bbox_inches='tight')
Дивіться схожі статті за категоріями: Python, Physics