#python #sympy #ode
Вопрос:
Я пытался написать программу, которая будет принимать набор значений и начальных условий для гармонического генератора, решать дифференциальное уравнение и отображать функцию phi(t).
import sympy as sp
from matplotlib import pyplot as plt
import numpy as np
t = sp.Symbol("t")
phi = sp.Function('phi')
#Values of the oscilator
beta = 0.4
omega_0 = 0.4
f0 = 0
omega = 0
#Initial conditions
phi_0 = 3.14/4
v_phi_0 = 0
#phi''(t) 2*beta*phi'(t) omega_0^2*phi(t) - f0*sin(omega*t) = 0
eq = sp.diff(phi(t),t,2) 2*beta*sp.diff(phi(t),t) omega_0**2*phi(t)-f0*sp.sin(omega*t)
print(sp.dsolve(eq))
equation = sp.dsolve(eq,ics={phi(0):phi_0,sp.diff(phi(t),t).subs(t,0):v_phi_0})
print(equation)
wzor = sp.lambdify(t,equation.rhs,"numpy")
xvals = np.arange(0,30,.1)
yvals = wzor(xvals)
# make figure
fig, ax = plt.subplots(1,1)
ax.plot(xvals, yvals)
ax.set_xlabel('t')
ax.set_ylabel('phi(t)')
plt.show()
Кажется, все работает нормально, но когда значения beta и omega_0 равны (так называемое критическое затухание) и меньше 1, программа, похоже, неправильно решает уравнение.
Eq(phi(t), (C1*sin(3.65002571393103e-9*t) C2*sin(3.6500259065127e-9*t) C3*cos(3.65002571393103e-9*t) C4*cos(3.6500259065127e-9*t))*exp(-0.4*t))
В других случаях, в том числе, когда beta и omega_0 равны и больше 1, есть только 2 константы C1 и C2 после решения дифференциального уравнения, но в данном случае их 4, чего не должно быть. Из-за этого программа не может продолжить и нарисовать сюжет, из-за этих 2 дополнительных констант, для которых я не могу решить.
ax.plot(xvals, yvals)
TypeError: can't convert expression to float
Редактировать: Кроме того, при приложении движущей силы, so f0 и omega отличается от 0, возникает следующая ошибка
NotImplementedError: Cannot find 2 solutions to the homogeneous equation necessary to apply undetermined coefficients to 0.16*phi(t) - 5*sin(0.4*t) 0.8*Derivative(phi(t), t) Derivative(phi(t), (t, 2)) (number of terms != order)
Возможно, это связано с основной проблемой, но я подумал, что об этом стоит упомянуть.
Правка 2: Другие проблемы начинают появляться, когда beta и omega_0 равны и больше 3, но также являются десятичными. Никаких проблем не возникает, когда это целые числа.
Комментарии:
1. Частоты по существу равны нулю, их ненулевые значения являются артефактами числового формата с плавающей запятой. При возмущении двойного корня можно ожидать погрешности порядка квадратных корней точности станка. То, что существует четыре корня, возможно, является следствием какой-то техники преобразования Лапласа. В символических вычислениях лучше работать только с символическими константами или рациональными числами.
Ответ №1:
Символьные вычисления и константы с плавающей запятой плохо сочетаются или вообще не сочетаются друг с другом. Одно простое исправление состоит в замене всех констант рациональными числами с использованием типа Rational данных sympy,
#Values of the oscilator
beta = sp.Rational(2,5)
omega_0 = sp.Rational(2,5)
f0 = 0
omega = 0
#Initial conditions
phi_0 = sp.Rational(314,400)
v_phi_0 = 0
Только с учетом этих изменений вывод будет корректным
Eq(phi(t), (C1 C2*t)*exp(-2*t/5))
Eq(phi(t), (157*t/500 157/200)*exp(-2*t/5))
и сюжет создается как
