dz.py - master - AlexParkes/Group_Homework_Python

dz.py
211 строк9 KB
Mikhail Ershov
Оформил первую строчку
09 ноя 2025, 22:29
09 ноя 2025, 22:29f5d67d5
#-----------------БИБЛИОТЕКИ--------------------------------------------------------#
from scipy.integrate import solve_ivp    # Чтобы можно было решить ДУ
import matplotlib.pyplot as plt          # Для построения графиков
import numpy as np                       # Для разных вычислений
from scipy.fft import fft                # Для преобразования Фурье
from scipy import signal                 # Нужно для работы с сигналами и фильтрами

# Задаём переменные со значениями по умолчанию
default_m1 = 1
default_m2 = 2
default_c1 = 200
default_c2 = 300
default_c3 = 500

# Проверка введенных пользователем значений ИД
def err01(arg_m1):
    if arg_m1 == '' or abs(float(arg_m1)) <= 10 ** (-6): 
        return default_m1
    else: 
        return arg_m1

def err02(arg_m2):
    if arg_m2 == '' or abs(float(arg_m2)) <= 10 ** (-6): 
        return default_m2
    else: 
        return arg_m2

def err03(arg_c1):
    if arg_c1 == '' or abs(float(arg_c1)) <= 10 ** (-6):
        return default_c1
    else: 
        return arg_c1 

def err04(arg_c2):
    if arg_c2 == '' or abs(float(arg_c2)) <= 10 ** (-6): 
        return default_c2
    else: 
        return arg_c2 

def err05(arg_c3):
    if arg_c3 == '' or abs(float(arg_c3)) <= 10 ** (-6): 
        return default_c3
    else: 
        return arg_c3

# Ввод ИД
m1 = float(err01(input()))
m2 = float(err02(input()))
c1 = float(err03(input()))
c2 = float(err04(input()))
c3 = float(err05(input()))

# Вывод ИД
print(f"Масса тела 1 составляет {m1} кг")
print(f"Масса тела 2 составляет {m2} кг")
print(f"Жёсткость пружины 1 составляет {c1} Н/м")
print(f"Жёсткость пружины 2 составляет {c2} Н/м")
print(f"Жёсткость пружины 3 составляет {c3} Н/м")

# Коэффициенты a11, a12, a22, c11, c12, c22
a11 = 1.5 * m1
a12 = 0
a22 = m2
c11 = c1
c12 = c1
c22 = c1 + c2 + c3

print(f"Коэффициент a11 = {a11}")
print(f"Коэффициент a12 = {a12}")
print(f"Коэффициент a22 = {a22}")
print(f"Коэффициент c11 = {c11}")
print(f"Коэффициент c12 = {c12}")
print(f"Коэффициент c22 = {c22}")

# Параметры интегрирования
del_t = 0.001           # шаг интегрирования
tk = 5                  # время интегрирования
t_span = (0, tk)
t_eval = np.arange(0, tk, del_t)

# Учёт гармонической силы в МКС без фильтрации
def system_with_F(t, y):
    dydt = np.zeros(4)
    dydt[0] = y[1]
    dydt[1] = -c11*y[0]/a11 - c12*y[2]/a11 - F*np.cos(p*t)/a11 - Ap*np.random.randn()/a11
    dydt[2] = y[3]
    dydt[3] = -c12*y[0]/a22 - c22*y[2]/a22
    return dydt

# Выводы из аналитического расчёта из НИРС
D = (c1 + c2 + c3)**2 / m2**2 - (4*c1*(c2 + c3) - 4*c1**2)/(3*m1*m2) + 4*c1**2/(9*m1**2) # расчёт дискриминанта
w1 = ((3*m1*(c1+c2+c3)+2*m2*c1)/(2*3*m1*m2) - D**(1/2)/2)**(1/2)
w2 = ((3*m1*(c1+c2+c3)+2*m2*c1)/(2*3*m1*m2) + D**(1/2)/2)**(1/2)
print(f"{w1: 0.3f} рад/с")
print(f"{w2: 0.3f} рад/с")

# НИРС - Без учёта гармонической силы и без шумовой составляющей.
F = 0
p = 0
Ap = 0

# Функция системы без фильтрации
def system_nirs(t, y):
    dydt = np.zeros(4)
    dydt[0] = y[1]  # скорость тела 1
    dydt[1] = -c11*y[0]/a11 - c12*y[2]/a11 - F*np.cos(p*t)/a11 - Ap*np.random.randn()/a11
    dydt[2] = y[3]  # скорость тела 2
    dydt[3] = -c12*y[0]/a22 - c22*y[2]/a22
    return dydt

# Решение системы без фильтрации
sol1 = solve_ivp(system_nirs, t_span, [1, 0, 0, 0], t_eval=t_eval)
T = sol1.t
Y1_no_filter = sol1.y[0]  # перемещение тела 1
Y2_no_filter = sol1.y[2]  # перемещение тела 2

# Графики колебаний без фильтрации
plt.figure(figsize=(10, 6))
plt.plot(T, Y1_no_filter, label='Колебания тела 1')
plt.plot(T, Y2_no_filter, label='Колебания тела 2')
plt.title('Колебания механической системы (без фильтрации) 1 случай')
plt.xlabel('Время t, с')
plt.ylabel('qi, м')
plt.legend()
plt.grid(True)
plt.show()

# НИРС - Учёт гармонической силы без шумовой составляющей.
F = 30          # амплитуда гармонической возмущающей силы
p = 3           # частота возмущающей силы
Ap = 0      # амплитуда шумовой составляющей к гармонической силе

# Решение системы без фильтрации
sol1 = solve_ivp(system_nirs, t_span, [1, 0, 0, 0], t_eval=t_eval)
T = sol1.t
Y1_no_filter = sol1.y[0]  # перемещение тела 1
Y2_no_filter = sol1.y[2]  # перемещение тела 2

# Графики колебаний без фильтрации
plt.figure(figsize=(10, 6))
plt.plot(T, Y1_no_filter, label='Колебания тела 1')
plt.plot(T, Y2_no_filter, label='Колебания тела 2')
plt.title('Колебания механической системы (без фильтрации) 2 случай')
plt.xlabel('Время t, с')
plt.ylabel('qi, м')
plt.legend()
plt.grid(True)
plt.show()

# НИРС - Учёт гармонической силы с шумовой составляющей.
F = 30          # амплитуда гармонической возмущающей силы
p = 3           # частота возмущающей силы
Ap = 30      # амплитуда шумовой составляющей к гармонической силе

# Решение системы без фильтрации
sol1 = solve_ivp(system_nirs, t_span, [1, 0, 0, 0], t_eval=t_eval)
T = sol1.t
Y1_no_filter = sol1.y[0]  # перемещение тела 1
Y2_no_filter = sol1.y[2]  # перемещение тела 2

# Графики колебаний без фильтрации
plt.figure(figsize=(10, 6))
plt.plot(T, Y1_no_filter, label='Колебания тела 1')
plt.plot(T, Y2_no_filter, label='Колебания тела 2')
plt.title('Колебания механической системы (без фильтрации) 3 случай')
plt.xlabel('Время t, с')
plt.ylabel('qi, м')
plt.legend()
plt.grid(True)
plt.show()

# Построение АЧХ без фильтрации
Fs = 1000            # Частота дискретизации
L = 5000             # Длина сигнала
t_fft = np.arange(0, L)/Fs    # Период дискретизации
Y1_fft = Y1_no_filter[:L] if len(Y1_no_filter) >= L else Y1_no_filter
Y2_fft = Y2_no_filter[:L] if len(Y2_no_filter) >= L else Y2_no_filter

# Преобразование Фурье
Aw1 = fft(Y1_fft)    # Преобразование Фурье для первого тела
Aw2 = fft(Y2_fft)    # Преобразование Фурье для второго тела
P2Aw1 = np.abs(Aw1 / L)    # Вычисление двустороннего спектра (тело 1)
P2Aw2 = np.abs(Aw2 / L)    # Вычисление двустороннего спектра (тело 2)
P1Aw1 = P2Aw1[:L//2 + 1]     # Вычисление одностороннего спектра (тело 1)
P1Aw2 = P2Aw2[:L//2 + 1]     # Вычисление одностороннего спектра (тело 2)
P1Aw1[1:-1] *= 2     # Удвоение амплитуды кроме нулевой и последней частот
P1Aw2[1:-1] *= 2     # Удвоение амплитуды кроме нулевой и последней частот

# Массив частот на АЧХ
f = Fs * np.arange(0, L//2 + 1)/L

# Вывод графиков АЧХ
plt.figure(figsize=(10, 6))
plt.plot(f, P1Aw1, label='АЧХ тела 1')
plt.plot(f, P1Aw2, label='АЧХ тела 2')
plt.title('АЧХ без фильтрации')
plt.xlabel('w, Гц')
plt.ylabel('A, м')
plt.legend()
plt.grid(True)
plt.show()

# Первый коммит Жамили

# Создаем ФНЧ - фильтр нижних частот
fsrez = 50                        # Частота среза
shag = del_t / 10                 # Шаг
n = 4                             # Порядок фильтра в виде целочисленного скаляра
Rp = 0.8                          # Неравномерность в полосе пропускания от пика к пику
Wp = fsrez/(1 / (shag * 2))       # Частота ребра полосы пропускания
bc, ac = signal.cheby1(n, Rp, Wp, 'low')      #Фильтр Чебышева