Символьное решение линейных дифференциальных уравнений и систем методом преобразований Лапласа c применением SymPy13.10.2018 14:02

17x9uwkvvlnn2gxa0wc3hsonvnk.jpeg

Реализация алгоритмов на языке Python с использованием символьных вычислений очень удобна при решении задач математического моделирования объектов, заданных дифференциальными уравнениями. Для решения таких уравнений широко используются преобразования Лапласа, которые, говоря упрощенно, позволяют свести задачу к решению простейших алгебраических уравнений.
В данной публикации предлагаю рассмотреть функции прямого и обратного преобразования Лапласа из библиотеки SymPy, которые позволяют использовать метод Лапласа для решения дифференциальных уравнений и систем средствами Python.
Сам метод Лапласа и его преимущества при решении линейных дифференциальных уравнений и систем широко освещены в литературе, например в популярном издании [1]. В книге метод Лапласа приведен для реализации в лицензионных программных пакетах Mathematica, Maple и MATLAB (что подразумевает приобретение учебным заведением этого ПО) на выбранных автором отдельных примерах.
Попробуем сегодня рассмотреть не отдельный пример решения учебной задачи средствами Python, а общий метод решения линейных дифференциальных уравнений и систем с использованием функций прямого и обратного преобразования Лапласа. При этом сохраним обучающий момент: левая часть линейного дифференциального уравнения с условиями Коши будет формироваться самим студентом, а рутинная часть задачи, состоящая в прямом преобразовании Лапласа правой части уравнения, будет выполняться при помощи функции laplace_transform ().

История об авторстве преобразований Лапласа


Преобразования Лапласа (изображения по Лапласу) имеют интересную историю. Впервые интеграл в определении преобразования Лапласа появился в одной из работ Л. Эйлера. Однако в математике общепринято называть методику или теорему именем того математика, который открыл ее после Эйлера. В противном случае существовало бы несколько сотен различных теорем Эйлера.
В данном случае следующим после Эйлера был французский математик Пьер Симон де Лаплас (Pierre Simon de Laplace (1749–1827)). Именно он использовал такие интегралы в своей работе по теории вероятностей. Самим Лапласом не применялись так называемые «операционные методы» для нахождения решений дифференциальных уравнений, основанные на преобразованиях Лапласа (изображениях по Лапласу). Эти методы в действительности были обнаружены и популяризировались инженерами-практиками, особенно английским инженером-электриком Оливером Хевисайдом (1850–1925). Задолго до того, как была строго доказана справедливость этих методов, операционное исчисление успешно и широко применялось, хотя его законность ставилось в значительной мере под сомнение даже в начале XX столетия, и по этой теме велись весьма ожесточенные дебаты.

Функции прямого и обратного преобразования Лапласа


  1. Функция прямого преобразования Лапласа:
    sympy.integrals.transforms.laplace_transform(t, s, **hints).
    Функция laplace_transform () выполняет преобразования Лапласа функции f (t) вещественной переменной в функцию F (s) комплексной переменной, так что:

    $F(s)=\int_{0}^{∞}f(t)e^{-st}dt.$


    Эта функция возвращает (F, a, cond), где F (s) есть преобразование Лапласа функции f (t), a — полуплоскость, где определена F (s), и cond— вспомогательные условия сходимости интеграла.
    Если интеграл не может быть вычислен в закрытой форме, эта функция возвращает невычисленным LaplaceTransform объекта.
    Если установить опцию noconds=True, функция возвращает только F (s).
  2. Функция обратного преобразования Лапласа:
    sympy.integrals.transforms.inverse_laplace_transform(F, s, t, plane=None, **hints).
    Функция inverse_laplace_transform () осуществляет обратное преобразование Лапласа функции комплексного переменного F (s) в функцию f (t) вещественной переменной, так что:

    $f(t)=\frac{1}{2\pi i}\int_{c-∞\cdot i}^{c+∞\cdot i}e^{st}F(s)ds.$


    Если интеграл не может быть вычислен в закрытой форме, эта функция возвращает невычисленным InverseLaplaceTransform объекта.


Обратное преобразование Лапласа на примере определения переходной характеристики регуляторов


Передаточная функция ПИД-регулятора имеет вид [2]:

$W(s)=(1+\frac{K_{d}\cdot T_{d}\cdot s}{1+T_{d}\cdot s})\cdot K_{p}\cdot (1+\frac{1}{T_{i}\cdot s})\cdot \frac{1}{s}.$


Напишем программу для получения уравнений для переходных характеристик ПИД- и ПИ-регуляторов для приведенной передаточной функции, и дополнительно выведем время, затраченное на выполнение обратного визуального преобразования Лапласа.

Текст программы 
# Загрузка необходимых модулей
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
eimqdjovtna0wpsekvms9mig_xk.png
Обратное преобразование Лапласа часто используется при синтезе САУ, где Python может заменить дорогостоящих программных «монстров» типа MathCAD, поэтому приведенное использование обратного преобразования имеет практическое значение.

Преобразование Лапласа от производных высших порядков для решения задачи Коши


В продолжение нашего обсуждения будет приложение преобразований Лапласа (изображений по Лапласу) к поиску решений линейного дифференциального уравнения с постоянными коэффициентами вида:

$a\cdot x''(t)+b\cdot x'(t)+c\cdot x(t)=f(t). (1)$

Если a и b — константы, то

$L\left\{a\cdot f(t)+b\cdot q(t)\right\}=a\cdot L\left\{f(t)\right\}+b\cdot L\left\{q(t)\right\} (2)$


для всех 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)

Приведенный код также демонстрирует однозначность обратного преобразования Лапласа.

Если предположить, что $q(t)=f'(t)$ удовлетворяет условиям первой теоремы, то из этой теоремы будет следовать, что:

$L\left\{f''(t)\right\}=L\left\{q'(t)\right\}=sL\left\{q(t)\right\}-q(0)=sL\left\{f'(t)\right\}-f'(0)=s\left[sL\left\{f(t)-f(0)\right\}\right],$


и, таким образом,

$L\left\{f''(t)\right\}=s^{2}\cdot F(s)-s\cdot f(0)-f'(0).$


Повторение этого вычисления дает

$L\left\{f'''(t)\right\}=sL\left\{f''(t)\right\}-f''(0)=s^{3}F(s)-s^{2}f(0)-sf'(0)-f''(0).$


После конечного числа таких шагов мы получаем следующее обобщение первой теоремы:

$L\left\{f^{(n)}(t)\right\}=s^{n}L\left\{f(t)\right\}-s^{n-1}f(0)-s^{n-2}f'(0)-\cdot\cdot\cdot -f^{(n-1)}(0)=$

$=s^{n}F(s)-s^{n-1}f(0)-\cdot\cdot\cdot -sf^{(n-2)}(0)-f^{(n-1)}(0). (3)$

Применяя соотношение (3), содержащее преобразованные по Лапласу производные искомой функции с начальными условиями, к уравнению (1), можно получить его решение по методу, специально разработанному на нашей кафедре при активной поддержке Scorobey для библиотеки SymPy.

Метод решения линейных дифференциальных уравнений и систем уравнений, основанный на преобразованиях Лапласа, с использованием библиотеки SymPy


Для демонстрации метода используем простое дифференциальное уравнение, описывающее движение системы, состоящей из материальной точки заданной массы, закрепленной на пружине, к которой приложена внешняя сила. Дифференциальное уравнение и начальные условия для такой системы имеют вид:

$x''+4x=\sin(3t); x(0)=1.2; x'(0)=1, (4)$


где $x(0)$ — приведенное начальное положение массы, $x'(0)$ — приведенная начальная скорость массы.

Упрощённая физическая модель, заданная уравнением (4) при ненулевых начальных условиях [1]:
stf-v71tfn2unhr7393cpmfw0oc.jpeg

Система, состоящая из материальной точки заданной массы, закрепленной на пружине, удовлетворяет задаче Коши (задаче с начальными условиями). Материальная точка заданной массы первоначально находится в покое в положении ее равновесия.

Для решения этого и других линейных дифференциальных уравнений методом преобразований Лапласа удобно пользоваться следующей системой, полученной из соотношений (3):
$L\left\{f^{IV}(t)\right\}=s^{4}\cdot F(s)-s^{3}\cdot f(0)-s^{2}\cdot f^{I}(0)-s\cdot f^{II}(0)-f^{III}(0),$
$L\left\{f^{III}(t)\right\}=s^{3}\cdot F(s)-s^{2}\cdot f(0)-s\cdot f^{I}(0)-f^{II}(0),$
$L\left\{f^{II}(t)\right\}=s^{2}\cdot F(s)-s\cdot f(0)-f^{I}(0),$
$L\left\{f^{I}(t)\right\}=s\cdot F(s)-f(0), L\left\{f(t)\right\}=F(s).$
$L\left\{f(t)\right\}=F(s). (5) $

Последовательность решения средствами SymPy следующая:

  1. загружаем необходимые модули и явно определяем символьные переменные: 
    from sympy import *
import numpy as np
import matplotlib.pyplot as plt
var('s')
var('t', positive=True)
var('X', cls = Function)

  2. указываем версию библиотеки sympy, чтобы учесть ее особенности. Для этого нужно ввести такие строки: 
    import SymPy
print('Версия библиотеки sympy – %' %(sympy._version_))

  3. по физическому смыслу задачи переменная времени определяется для области, включающей ноль и положительные числа. Задаём начальные условия и функцию в правой части уравнения (4) с её последующим преобразование по Лапласу. Для начальных условий необходимо использовать функцию Rational, поскольку использование десятичного округления приводит к ошибке.
    x0=Rational(6,5) # Приведенное начальное положение массы
x01=Rational (1,1) # Приведенная начальная скорость
g=sin(3*t)
Lg=laplace_transform( g, t, s, noconds = True)

  4. пользуясь (5), переписываем преобразованные по Лапласу производные, входящие в левую часть уравнения (4), формируя из них левую часть этого уравнения, и сравниваем результат с правой его частью: 
    d2=s**2*X(s)-s*x0-x01
d0=X(s)
d=d2+4*d0
de=Eq(d,Lg)

  5. решаем полученное алгебраическое уравнение относительно преобразования X (s) и выполняем обратное преобразование Лапласа: 
    rez=solve(de,X(s))[0]
soln=inverse_laplace_transform(rez,s,t)

  6. осуществляем переход из работы в библиотеке SymPyв библиотеку NumPy: 
    f=lambdify(t,soln,'numpy')

  7. строим график обычным для Python методом: 
    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

njztkaopwkurg1k0cnh4nucbzf8.png
Получен график периодической функции, дающей положение материальной точки заданной массы. Метод преобразования Лапласа с использованием библиотеки SymPy дает решение не только без потребности сначала найти общее решение однородного уравнения и частное решение первоначального неоднородного дифференциального уравнения, но и без потребности использования метода элементарных дробей и таблиц Лапласа.
При этом учебное значение метода решения сохраняется за счёт необходимости использования системы (5) и перехода в NumPy для исследования решения более производительными методами.

Для дальнейшей демонстрации метода решим систему дифференциальных уравнений:
$\begin{cases}2x''+6x-2=0,\\y''-2x+2y=40\cdot \sin(3t),\end{cases} (6)$
с начальными условиями $x(0)=y(0)=y'(0)=0.$

Упрощённая физическая модель, заданная системой уравнений (6) при нулевых начальных условиях:
lqyfz7mwkd6eoesoou5ygcxngna.jpeg

Таким образом, сила 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)=%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()


Получим:
mv4lq6dnyzolk7_1hnuqut5e4hc.png

Для ненулевых начальных условий текст программы и график функций примет вид:

Текст программы 
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()


piltdbzyf4n58i4ek8wipsdruza.png

Рассмотрим решение линейного дифференциального уравнения четвёртого порядка с нулевыми начальными условиями:
$x^{(4)}+2\cdot x''+x=4\cdot t\cdot e^{t};$
$x(0)=x'(0)=x''(0)=x^{(3)}(0)=0.$

Текст программы: 
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)=%s\n' %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()


График решения:
vhaf82hrfe42xunqusxcwyn4k3w.png

Решим линейное дифференциальное уравнение четвёртого порядка:
$x^{(4)}+13x''+36x=0;$
с начальными условиями $x(0)=x''(0)=0$, $x'(0)=2$, $x^{(3)}(0)=-13$.

Текст программы: 
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)=%s\n' %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()


График решения:
lbuvrsjxqc_kvev-njvnlt8jt9a.png

Функции для решения ОДУ


Для имеющих аналитическое решение ОДУ и систем ОДУ применяется функция dsolve ():
sympy.solvers.ode.dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)

Давайте сравним производительность функции dsolve () с методом Лапласа. Для примера возьмём следующее дифференциальное уравнение четвёртой степени с нулевыми начальными условиями:
$x^{(IV)}(t)=3\cdot x'''(t)-x(t)=4\cdot t\cdot \exp(t);$
$x(0)=x'(0)=x''(0)=x'''(0)=0.$

Программа с использованием функции 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)=%s\n' %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
gcwmei24uj-q5cjqjnnzexy_ttq.png

Программа с использованием преобразований Лапласа: 
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)=%s\n' %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
deoxib8hfyfwpd12s0buzjnnwgk.png

Итак, функция 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. Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления

© Habrahabr.ru