Задача
В статья использованы возможности пакета SymPy совместно с пакетом NumPy. Всё сводиться к преобразованию символьных выражений в функции способные работать с другими модулями Python.
Процесс решения дифференциальных уравнений становиться наглядным и хорошо контролируемым на каждом этапе вычислений. Следует отметить, что колебательное звено в разных интерпретациях обсуждается в сетях [1,2]. Например, в [3] приводиться модель колебательного звена с подробным исследованием переходных процессов.
Надеюсь, что подобные исследования колебательного звена на Python найдут своих сторонников.
Код программы
Чтобы не утомлять читателя приведу сразу код программы с пояснением каждой значимой строки.
Код
import numpy as np from sympy import * from IPython.display import * import matplotlib.pyplot as plt init_printing(use_latex=True) var('t C1 C2') u = Function("u")(t) # Это переменная, но не функция. m=20 #Показатель массы. w=10.0#Показатель демпфирования колебаний. c=0.3#Показатель жёсткости. a=1#Бесконечный импульс силы. #Все показатели условные(только для исследования характера зависимостей). t#Текущее время. de = Eq(m*u.diff(t,t)+w*u.diff(t)+c*u,a) #-Дифференциальное уравнение колебаний. display(de)#-Вывод на дисплей. des = dsolve(de,u)#Символьное решение уравнения методом Коши в общем виде. display(des)#Вывод на дисплей. eq1=des.rhs.subs(t,0);#Условие равенства нулю перемещения в момент времени t=0. display(eq1)#Вывод на дисплей. eq2=des.rhs.diff(t).subs(t,0)#Условие равенства нулю скорости перемещения в момент # времени t=0. display(eq2)#Вывод на дисплей. seq=solve([eq1,eq2],C1,C2)#Решение системы для определения коэффициентов C1,C2. display(seq)#Вывод на дисплей. rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2])])#Вид решения дифференциального #уравнения с численными значениями коэффициентов. display(rez)#Вывод на дисплей. f=lambdify(t, rez, "numpy")#Перевод символьного решения в численное для работы #с модулем numpy . x = np.arange(0.0,500,0.01) plt.plot(x,f(x),color='b', linewidth=3) plt.xlabel('Time t seconds',fontsize=12) plt.ylabel('$f(t)$',fontsize=16) plt.grid(True) plt.show()
Настроим демпфирование и получим апериодическое движение и все этапы решения уравнения.
Eq(0.3*u(t) + 10.0*Derivative(u(t), t) + 20*Derivative(u(t), t, t), 1)
Eq(u(t), C1*exp(t*(-5 — sqrt(19))/20) + C2*exp(t*(-5 + sqrt(19))/20) + 3.33333333333333)
C1 + C2 + 3.33333333333333
C1*(-1/4 — sqrt(19)/20) + C2*(-1/4 + sqrt(19)/20)
{C1: 0.245131115588015, C2: -3.57846444892135}
0.245131115588015*exp(t*(-5 — sqrt(19))/20) — 3.57846444892135*exp(t*(-5 + sqrt(19))/20) + 3.33333333333333
Изменим параметры и перепишем листинг программы для отображения трёх движений с возрастающим демпфированием на одном графике.
Код программы для трех различных значений показателя демпфирования
import numpy as np from sympy import * from IPython.display import * import matplotlib.pyplot as plt init_printing(use_latex=True) var('t C1 C2') u = Function("u")(t) # Это переменная, но не функция. m=200 #Показатель массы. w=1.8#Показатель демпфирования колебаний. c=1.3#Показатель жёсткости. a=1#Бесконечный импульс силы. #Все показатели условные (только для исследования характера зависимостей). t#Текущее время. man=[0.8,2.12,5]#Список условных значений демпфирования (только для построения графиков). for w in man: de = Eq(m*u.diff(t,t)+w*u.diff(t)+c*u,a) #-Дифференциальное уравнение колебаний. display(de)#-Вывод на дисплей. des = dsolve(de,u)#Символьное решение уравнения методом Коши в общем виде. display(des)#Вывод на дисплей. eq1=des.rhs.subs(t,0);#Условие равенства нулю перемещения в момент времени t=0. display(eq1)#Вывод на дисплей. eq2=des.rhs.diff(t).subs(t,0)#Условие равенства нулю скорости перемещения в момент времени t=0. display(eq2)#Вывод на дисплей. seq=solve([eq1,eq2],C1,C2)#Решение системы для определения коэффициентов C1,C2. display(seq)#Вывод на дисплей. rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2])])#Вид решения дифференциального уравнения #с численными значениями коэффициентов. display(rez)#Вывод на дисплей. f=lambdify(t, rez, "numpy")#Перевод символьного решения в численное для работы с модулем numpy . x = np.arange(0.0,500,0.01) if w==man[0]:#Три перемещения на одном графике. plt.plot(x,f(x),color='r', linewidth=3) elif w==man[1]: plt.plot(x,f(x),color='g', linewidth=3) elif w==man[2]: plt.plot(x,f(x),color='b', linewidth=3) plt.xlabel('Time t seconds',fontsize=12) plt.ylabel('$f(t)$',fontsize=16) plt.grid(True) plt.show()
Установим значение параметра демпфирования и получим график периодического затухающего движения и все этапы решения трех уравнений.
Eq(1.3*u(t) + 0.8*Derivative(u(t), t) + 200*Derivative(u(t), t, t), 1)
Eq(u(t), (C1*sin(sqrt(406)*t/250) + C2*cos(sqrt(406)*t/250))/exp(t)**(1/500) + 0.769230769230769)
C2 + 0.769230769230769
sqrt(406)*C1/250 — C2/500
{C2: -0.769230769230769, C1: -0.0190881410379025}
(-0.0190881410379025*sin(sqrt(406)*t/250) — 0.769230769230769*cos(sqrt(406)*t/250))/exp(t)**(1/500) + 0.769230769230769
Eq(1.3*u(t) + 2.12*Derivative(u(t), t) + 200*Derivative(u(t), t, t), 1)
Eq(u(t), (C1*sin(sqrt(647191)*t/10000) + C2*cos(sqrt(647191)*t/10000))/exp(t)**(53/10000) + 0.769230769230769)
C2 + 0.769230769230769
sqrt(647191)*C1/10000 — 53*C2/10000
{C2: -0.769230769230769, C1: -0.0506776284001243}
(-0.0506776284001243*sin(sqrt(647191)*t/10000) — 0.769230769230769*cos(sqrt(647191)*t/10000))/exp(t)**(53/10000) + 0.769230769230769
Eq(1.3*u(t) + 5*Derivative(u(t), t) + 200*Derivative(u(t), t, t), 1)
Eq(u(t), (C1*sin(sqrt(1015)*t/400) + C2*cos(sqrt(1015)*t/400))/exp(t)**(1/80) + 0.769230769230769)
C2 + 0.769230769230769
sqrt(1015)*C1/400 — C2/80
{C2: -0.769230769230769, C1: -0.120724003956605}
(-0.120724003956605*sin(sqrt(1015)*t/400) — 0.769230769230769*cos(sqrt(1015)*t/400))/exp(t)**(1/80) + 0.769230769230769
Уменьшим показатель массы и получим графики движения (решение уравнений для сокращения упускаем).
Рассмотрим важный случай отрицательного демпфирования и получим соответствующий график (решение уравнений для сокращения упускаем).
Вывод
Полученная графическая иллюстрация работы колебательного звена полностью соответствует теории.
Символьное и численное решений дифференциального уравнения на SymPy и NumPy позволило создать адекватную модель колебательного звена с управление массой демпфированием и жёсткостью и при этом контролировать ход ращения уравнения колебаний. Кроме этого Python условно бесплатное программное обеспечение, что значительно расширяет возможности применение приведенной программы в исследовательских и учебных целях.
