- PVSM.RU - https://www.pvsm.ru -
Реализация алгоритмов на языке Python с использованием символьных вычислений очень удобна при решении задач математического моделирования объектов, заданных дифференциальными уравнениями. Для решения таких уравнений широко используются преобразования Лапласа, которые, говоря упрощенно, позволяют свести задачу к решению простейших алгебраических уравнений.
В данной публикации предлагаю рассмотреть функции прямого и обратного преобразования Лапласа из библиотеки SymPy, которые позволяют использовать метод Лапласа для решения дифференциальных уравнений и систем средствами Python.
Сам метод Лапласа и его преимущества при решении линейных дифференциальных уравнений и систем широко освещены в литературе, например в популярном издании [1]. В книге метод Лапласа приведен для реализации в лицензионных программных пакетах Mathematica, Maple и MATLAB (что подразумевает приобретение учебным заведением этого ПО) на выбранных автором отдельных примерах.
Попробуем сегодня рассмотреть не отдельный пример решения учебной задачи средствами Python, а общий метод решения линейных дифференциальных уравнений и систем с использованием функций прямого и обратного преобразования Лапласа. При этом сохраним обучающий момент: левая часть линейного дифференциального уравнения с условиями Коши будет формироваться самим студентом, а рутинная часть задачи, состоящая в прямом преобразовании Лапласа правой части уравнения, будет выполняться при помощи функции laplace_transform().
Преобразования Лапласа (изображения по Лапласу) имеют интересную историю. Впервые интеграл в определении преобразования Лапласа появился в одной из работ Л. Эйлера. Однако в математике общепринято называть методику или теорему именем того математика, который открыл ее после Эйлера. В противном случае существовало бы несколько сотен различных теорем Эйлера.
В данном случае следующим после Эйлера был французский математик Пьер Симон де Лаплас (Pierre Simon de Laplace (1749-1827)). Именно он использовал такие интегралы в своей работе по теории вероятностей. Самим Лапласом не применялись так называемые «операционные методы» для нахождения решений дифференциальных уравнений, основанные на преобразованиях Лапласа (изображениях по Лапласу). Эти методы в действительности были обнаружены и популяризировались инженерами-практиками, особенно английским инженером-электриком Оливером Хевисайдом (1850-1925). Задолго до того, как была строго доказана справедливость этих методов, операционное исчисление успешно и широко применялось, хотя его законность ставилось в значительной мере под сомнение даже в начале XX столетия, и по этой теме велись весьма ожесточенные дебаты.
Эта функция возвращает (F, a, cond), где F(s) есть преобразование Лапласа функции f(t), a<Re(s) – полуплоскость, где определена F(s), и cond– вспомогательные условия сходимости интеграла.
Если интеграл не может быть вычислен в закрытой форме, эта функция возвращает невычисленным LaplaceTransform объекта.
Если установить опцию noconds=True, функция возвращает только F(s).
Если интеграл не может быть вычислен в закрытой форме, эта функция возвращает невычисленным InverseLaplaceTransform объекта.
Передаточная функция ПИД-регулятора имеет вид [2]:
Напишем программу для получения уравнений для переходных характеристик ПИД- и ПИ-регуляторов для приведенной передаточной функции, и дополнительно выведем время, затраченное на выполнение обратного визуального преобразования Лапласа.
# Загрузка необходимых модулей
import time
start = time.time()
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
var('s Kp Ti Kd Td ') # Объявляем символьные переменные
var('t ', positive=True) # Ограничение на символьную переменную времени
Kp=2;Ti=2;Kd=4;Td=Rational(1,2) # Представление дробной символьной переменной
fp=(1+(Kd*Td*s)/(1+Td*s))*Kp*(1+1/(Ti*s))*(1/s) # Передаточная функция ПИД-регулятора с оператором s
ht=inverse_laplace_transform(fp,s,t) # Переходная характеристика ПИД-регулятора, получаемая методом обратного преобразования Лапласа
Kd=0
fpp=(1+(Kd*Td*s)/(1+Td*s))*Kp*(1+1/(Ti*s))*(1/s) # Передаточная функция ПИ-регулятора (Kd=0) с оператором s
htt=inverse_laplace_transform(fpp,s,t) # Переходная характеристика ПИ-регулятора, получаемая методом обратного преобразования Лапласа
stop = time.time()
print('Время на обратное визуальное преобразование Лапласа: %s s ' %N((stop-start),3))
# Переходим из символьной области в численную
f=lambdify(t, ht, 'numpy')
F=lambdify(t, htt, 'numpy')
tt = np.arange(0.01,20,0.05)
# Построение графика
plt.title('Переходные характеристики регуляторов n с передаточными функциями:n ПИД - W(s)=%s n ПИ - W(s)=%s' %(fp,fpp))
plt.plot(tt,f(tt), color='r', linewidth=2, label='ПИД-регулятор: h(t)=%s' %ht)
plt.plot(tt,F(tt), color='b', linewidth=2, label='ПИ-регулятор: h(t)=%s' %htt)
plt.grid(True)
plt.legend(loc='best')
plt.show()
Получим:
Время на обратное визуальное преобразование Лапласа: 2.68 s
Обратное преобразование Лапласа часто используется при синтезе САУ, где Python может заменить дорогостоящих программных “монстров” типа MathCAD, поэтому приведенное использование обратного преобразования имеет практическое значение.
В продолжение нашего обсуждения будет приложение преобразований Лапласа (изображений по Лапласу) к поиску решений линейного дифференциального уравнения с постоянными коэффициентами вида:
Если a и b — константы, то
для всех s, таких, что существуют оба преобразования Лапласа (изображения по Лапласу) функций f(t) и q(t).
Проверим линейность прямого и обратного преобразований Лапласа с помощью ранее рассмотренных функций laplace_transform() и inverse_laplace_transform(). Для этого в качестве примера примем f(t)=sin(3t), q(t)=cos(7t), a=5, b=7 и используем следующую программу.
from sympy import*
var('s a b')
var('t ', positive=True)
a=5
f=sin(3*t)
b=7
q=cos(7*t)
L1=a*laplace_transform(f, t, s, noconds = True) # Прямое преобразование a*L{f(t)}
L2=b*laplace_transform(q, t, s, noconds = True) # Прямое преобразование b*L{q(t)}
L=factor(L1+L2) # Сумма a*L{f(t)}+b*L{q(t)}
print(L)
LS=factor(laplace_transform(a*f+b*q, t, s, noconds = True)) # Прямое преобразование L{a*f(t)+b*q(t)}
print(LS)
print(LS==L)
L_1=a*inverse_laplace_transform(L1/a, s, t) # Обратное преобразование a* L^-1{f(t)}
L_2=b*inverse_laplace_transform(L2/b, s,t) # Обратное преобразование b* L^-1{q(t)}
L_S=L_1+L_2 # a* L^-1{f(t)}+b* L^-1{q(t)}
print(L_S)
L_1_2=inverse_laplace_transform(L1+L2,s,t) #Обратное преобразование L^-1{a*f(t)+b*q(t)}
print(L_1_2)
print(L_1_2==L_S)
Получим:
(7*s**3 + 15*s**2 + 63*s + 735)/((s**2 + 9)*(s**2 + 49))
(7*s**3 + 15*s**2 + 63*s + 735)/((s**2 + 9)*(s**2 + 49))
True
5*sin(3*t) + 7*cos(7*t)
5*sin(3*t) + 7*cos(7*t)
Приведенный код также демонстрирует однозначность обратного преобразования Лапласа.
Если предположить, что удовлетворяет условиям первой теоремы, то из этой теоремы будет следовать, что:
и, таким образом,
Повторение этого вычисления дает
После конечного числа таких шагов мы получаем следующее обобщение первой теоремы:
Применяя соотношение (3), содержащее преобразованные по Лапласу производные искомой функции с начальными условиями, к уравнению (1), можно получить его решение по методу, специально разработанному на нашей кафедре при активной поддержке Scorobey [1] для библиотеки SymPy.
Для демонстрации метода используем простое дифференциальное уравнение, описывающее движение системы, состоящей из материальной точки заданной массы, закрепленной на пружине, к которой приложена внешняя сила. Дифференциальное уравнение и начальные условия для такой системы имеют вид:
где — приведенное начальное положение массы, — приведенная начальная скорость массы.
Упрощённая физическая модель, заданная уравнением (4) при ненулевых начальных условиях [1]:
Система, состоящая из материальной точки заданной массы, закрепленной на пружине, удовлетворяет задаче Коши (задаче с начальными условиями). Материальная точка заданной массы первоначально находится в покое в положении ее равновесия.
Для решения этого и других линейных дифференциальных уравнений методом преобразований Лапласа удобно пользоваться следующей системой, полученной из соотношений (3):
Последовательность решения средствами SymPy следующая:
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t', positive=True)
var('X', cls = Function)
import SymPy
print('Версия библиотеки sympy – %' %(sympy._version_))
x0=Rational(6,5) # Приведенное начальное положение массы
x01=Rational (1,1) # Приведенная начальная скорость
g=sin(3*t)
Lg=laplace_transform( g, t, s, noconds = True)
d2=s**2*X(s)-s*x0-x01
d0=X(s)
d=d2+4*d0
de=Eq(d,Lg)
rez=solve(de,X(s))[0]
soln=inverse_laplace_transform(rez,s,t)
f=lambdify(t,soln,'numpy')
x=np.linspace(0,6*np.pi,100)
plt.title('Функция, дающая положение материальной точки n заданной массы:n х (t)=%s' %soln)
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t)',fontsize=12)
plt.plot(x,f(x),'g',linewidth=2)
plt.show()
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X', cls = Function)
import sympy
print("Версия библиотеки sympy – %s" %(sympy.__version__))
x0=Rational(6, 5) # Приведенное начальное положение материальной точки заданной массы
x01=Rational(1,1) # Приведенная начальная скорость материальной точки заданной массы
g=sin(3*t)
Lg=laplace_transform(g, t, s, noconds=True) # Прямое преобразование Лапласа
d2=s**2*X(s)-s*x0-x01
d0=X(s)
d=d2+4*d0
de=Eq(d, Lg)
rez=solve(de,X(s))[0] # Решение алгебраического уравнения
soln=inverse_laplace_transform(rez,s,t) # Обратное преобразование Лапласа
f=lambdify(t, soln ,"numpy")
x=np.linspace(0,6*np.pi,100)
plt.title('Функция, дающая положение материальной точки n заданной массы:n х (t)=%s' %soln)
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2)
plt.show()
Получаем:
Версия библиотеки sympy – 1.3
Получен график периодической функции, дающей положение материальной точки заданной массы. Метод преобразования Лапласа с использованием библиотеки SymPy дает решение не только без потребности сначала найти общее решение однородного уравнения и частное решение первоначального неоднородного дифференциального уравнения, но и без потребности использования метода элементарных дробей и таблиц Лапласа.
При этом учебное значение метода решения сохраняется за счёт необходимости использования системы (5) и перехода в NumPy для исследования решения более производительными методами.
Для дальнейшей демонстрации метода решим систему дифференциальных уравнений:
с начальными условиями
Упрощённая физическая модель, заданная системой уравнений (6) при нулевых начальных условиях:
Таким образом, сила f(t) внезапно прилагается ко второй материальной точке заданной массы в момент времени t = 0, когда система находится в покое в ее положении равновесия.
Решение системы уравнений идентично ранее рассмотренному решению дифференциального уравнения (4), поэтому привожу текст программы без пояснений.
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X Y', cls = Function)
x0 = 0
x01 =0
y0 = 0
y01 =0
g=40*sin(3*t)
Lg=laplace_transform( g, t, s, noconds = True)
de1=Eq(2*(s**2*X(s)-s*x0-x01)+6*X(s)-2*Y(s))
de2=Eq(s**2*Y(s)-s*y0-y01-2*X(s)+2*Y(s)-Lg)
rez=solve([de1,de2],X(s),Y(s))
rezX=expand(rez[X(s)])
solnX = inverse_laplace_transform(rezX,s,t)
rezY=expand(rez[Y(s)])
solnY = inverse_laplace_transform(rezY,s,t)
f=lambdify(t,solnX ,"numpy")
F=lambdify(t,solnY ,"numpy")
x = np.linspace(0,4*np.pi,100)
plt.title('Функции положения материальных точек заданных масс:n x(t)=%sn y(t)=%s' %(solnX,solnY))
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t), y(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2, label='x(t)')
plt.plot(x,F(x),'b', linewidth=2, label='y(t)')
plt.legend(loc='best')
plt.show()
Получим:
Для ненулевых начальных условий текст программы и график функций примет вид:
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X Y', cls = Function)
x0 =0
x01 =-1
y0 =0
y01 =-1
g=40*sin(t)
Lg=laplace_transform( g, t, s, noconds = True)
de1=Eq(2*(s**2*X(s)-s*x0-x01)+6*X(s)-2*Y(s))
de2=Eq(s**2*Y(s)-s*y0-y01-2*X(s)+2*Y(s)-Lg)
rez=solve([de1,de2],X(s),Y(s))
rezX=expand(rez[X(s)])
solnX = (inverse_laplace_transform(rezX,s,t)).evalf().n(3)
rezY=expand(rez[Y(s)])
solnY =( inverse_laplace_transform(rezY,s,t)).evalf().n(3)
f=lambdify(t,solnX ,"numpy")
F=lambdify(t,solnY ,"numpy")
x = np.linspace(0,4*np.pi,100)
plt.title( 'Функции положения материальных точек заданных масс:n x(t)= %s n y(t)=%s ' %(solnX,solnY))
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t), y(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2, label='x(t)')
plt.plot(x,F(x),'b', linewidth=2, label='y(t)')
plt.legend(loc='best')
plt.show()
Рассмотрим решение линейного дифференциального уравнения четвёртого порядка с нулевыми начальными условиями:
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X', cls = Function)
# Начальные условия
x0=0
x01=0
x02=0
x03=0
g = 4*t*exp(t)
Lg=laplace_transform( g, t, s, noconds = True) # Прямое преобразование Лапласа
d4=s**4*X(s)-s**3*x0-s**2*x01-s*x02-x03
d2=s**2*X(s)-s*x0-x01
d0=X(s)
d=factor(d4+2*d2+d0)
de= Eq(d, Lg)
rez = solve(de,X(s))[0] # Решение алгебраического уравнения
soln = inverse_laplace_transform(rez,s,t) # Обратное преобразование Лапласа
f=lambdify(t,soln ,"numpy")
x = np.linspace(0,6*np.pi,100)
plt.title( 'Решение:n х (t)=%sn' %soln,fontsize=11)
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2)
plt.show()
График решения:
Решим линейное дифференциальное уравнение четвёртого порядка:
с начальными условиями , , .
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X', cls = Function)
# Начальные условия
x0 =0
x01 =2
x02=0
x03=-13
d4=s**4*X(s)-s**3*x0-s**2*x01-s*x02-x03
d2=s**2*X(s)-s*x0-x01
d0=X(s)
d=factor(d4+13*d2+36*d0)
de = Eq(d, 0)
rez = solve(de,X(s))[0] # Решение алгебраического уравнения
soln = inverse_laplace_transform(rez,s,t) # Обратное преобразование Лапласа
f=lambdify(t,soln ,"numpy")
x = np.linspace(0,6*np.pi,100)
plt.title('Решение:n х (t)=%sn' %soln,fontsize=11)
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2)
plt.show()
График решения:
Для имеющих аналитическое решение ОДУ и систем ОДУ применяется функция dsolve():
sympy.solvers.ode.dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)
Давайте сравним производительность функции dsolve() с методом Лапласа. Для примера возьмём следующее дифференциальное уравнение четвёртой степени с нулевыми начальными условиями:
import time
start = time.time()
import numpy as np
from sympy import *
import matplotlib.pyplot as plt
var('t C1 C2 C3 C4')
u = Function("u")(t)
de = Eq(u.diff(t,t,t,t) -3*u.diff(t,t,t)+3*u.diff(t,t)-u.diff(t),4*t*exp(t)) # Запись дифференциального уравнения
des = dsolve(de,u) # Решение дифференциального уравнения
#Начальные условия
eq1=des.rhs.subs(t,0)
eq2=des.rhs.diff(t).subs(t,0)
eq3=des.rhs.diff(t,t).subs(t,0)
eq4=des.rhs.diff(t,t,t).subs(t,0)
# Решение системы алгебраических уравнений для начальных условий
seq=solve([eq1,eq2-1,eq3-2,eq4-3],C1,C2,C3,C4)
rez=des.rhs.subs([(C1,seq[C1]),(C2,seq[C2]),(C3,seq[C3]),(C4,seq[C4])])
F = Lambda(t,rez)
f=lambdify(t, rez, 'numpy')
x = np.linspace(0,6*np.pi,100)
stop = time.time()
print ('Время решения уравнения с использованием функции dsolve(): %s s' %round((stop-start),3))
plt.title('Решение с использованием функции dsolve():n х (t)=%sn' %rez,fontsize=11)
plt.grid(True)
plt.xlabel('Time t seconds',fontsize=12)
plt.ylabel('f(t)',fontsize=16)
plt.plot(x,f(x),color='#008000', linewidth=3)
plt.show()
Получим:
Время решения уравнения с использованием функции dsolve(): 1.437 s
import time
start = time.time()
from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t ', positive=True)
var('X', cls = Function)
#Начальные условия
x0 =0
x01 =0
x02=0
x03=0
g = 4*t*exp(t) # Запись левой части дифференциального уравнения
Lg=laplace_transform( g, t, s, noconds = True) # Прямое преобразование Лапласа
d4=s**4*X(s)-s**3*x0-s**2*x01-s*x02-x03
d3=s**3*X(s)-s**2*x0-s*x01-x02
d2=s**2*X(s)-s*x0-x01
d1=s*X(s)-x0
d0=X(s)
d=factor(d4-3*d3+3*d2- d1) # Запись правой части дифференциального уравнения
de= Eq(d, Lg)
rez = solve(de,X(s))[0] # Решение алгебраического уравнения
soln =collect(inverse_laplace_transform(rez,s,t),t) # Обратное преобразование Лапласа
f=lambdify(t,soln , 'numpy')
x = np.linspace(0,6*np.pi,100)
stop = time.time()
print ('Время решения уравнения с использованием преобразования Лапласа: %s s' %round((stop-start),3))
plt.title ('Решение с использованием преобразования Лапласа:n х (t)=%sn' %soln,fontsize=11)
plt.grid(True)
plt.xlabel(' t ',fontsize=12)
plt.ylabel('x(t)',fontsize=12)
plt.plot(x,f(x),'g', linewidth=2)
plt.show()
Получим:
Время решения уравнения с использованием преобразования Лапласа: 3.274 s
Итак, функция dsolve() (1.437 s) решает уравнение четвёртого порядка быстрее, чем выполняется решение по методу преобразований Лапласа (3.274 s) более чем в два раза. Однако при этом следует отметить, что функция dsolve() не решает системы дифференциальных уравнений второго порядка, например, при решении системы (6) с использованием функция dsolve() возникает ошибка:
from sympy import*
t = symbols('t')
x, y = symbols('x, y', Function=True)
eq = (Eq(Derivative(x(t),t,2),-3*x(t)+y(t)), Eq(Derivative(y(t),t,2),2*x(t)-2*y(t)+40*sin(3*t)))
rez=dsolve(eq)
print(list(rez))
Получим:
raiseNotImplementedError
NotImplementedError
Данная ошибка означает, что решение системы дифференциальных уравнений с помощью функции dsolve() не может быть представлено символически. Тогда как при помощи преобразований Лапласа мы получили символическое представление решения, и это доказывает эффективность предложенного метода.
Примечание.
Для того чтобы найти необходимый метод решения дифференциальных уравнений с помощью функции dsolve(), нужно использовать classify_ode(eq, f(x)), например:
from sympy import *
from IPython.display import *
import matplotlib.pyplot as plt
init_printing(use_latex=True)
x=Symbol('x')
f=Function('f')
eq=Eq(f(x).diff(x,x)+f(x),0)
print(dsolve(eq,f(x)))
print(classify_ode(eq, f(x)))
eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
print(classify_ode(eq, f(x)))
rez=dsolve(eq, hint='almost_linear_Integral')
print(rez)
Получим:
Eq(f(x), C1*sin(x) + C2*cos(x))
('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')
('separable', '1st_exact', 'almost_linear', '1st_power_series', 'lie_group', 'separable_Integral', '1st_exact_Integral', 'almost_linear_Integral')
[Eq(f(x), -acos((C1 + Integral(0, x))*exp(-Integral(-tan(x), x))) + 2*pi), Eq(f(x), acos((C1 + Integral(0,x))*exp(-Integral(-tan(x), x))))]
Таким образом, для уравнения eq=Eq(f(x).diff(x,x)+f(x),0) работает любой метод из первого списка:
nth_linear_constant_coeff_homogeneous,
2nd_power_series_ordinary
Для уравнения eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x) работает любой метод из второго списка:
separable, 1st_exact, almost_linear,
1st_power_series, lie_group, separable_Integral,
1st_exact_Integral, almost_linear_Integral
Чтобы использовать выбранный метод, запись функции dsolve() примет вид, к примеру:
rez=dsolve(eq, hint='almost_linear_Integral')
Данная статья ставила своей целью показать, как использовать средства библиотек SciPy и NumPy на примере решения систем линейных ОДУ операторным методом. Таким образом, были рассмотрены методы символьного решения линейных дифференциальных уравнений и систем уравнений методом Лапласа. Проведен анализ производительности этого метода и методов, реализованных в функции dsolve().
Ссылки:
1. Дифференциальные уравнения и краевые задачи: моделирование и вычисление с помощью Mathematica, Maple и MATLAB. 3-е издание.: Пер. с англ. — М.: ООО «И.Д. Вильяме», 2008. — 1104 с.: ил. — Парал. тит. англ.
2. Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления [2]
Автор: Елена Чёрная
Источник [3]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/python/295597
Ссылки в тексте:
[1] Scorobey: https://habr.com/users/scorobey/
[2] Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления: https://habr.com/post/346338/
[3] Источник: https://habr.com/post/426041/?utm_source=habrahabr&utm_medium=rss&utm_campaign=426041
Нажмите здесь для печати.