Что такое дифференциальные уравнения
Дифференциальные уравнения (ДУ) активно используются в моделировании динамических процессов, а также при анализе в широком спектре областей. Это могут быть как разные сферы науки, от популяционной экологии до гидродинамики, так и разные области экономики и менеджмента.
ДУ применяются для описания реальных процессов математически, преимущественно с целью прогнозирования и симулирования поведения, например:
- в механике ДУ описывают колебания струн, позволяют определять координаты и скорости объектов;
- в биологии — численность популяции биологических сообществ;
- в химии — кинетику химических реакций;
- в медицине — развитие эпидемий;
- в инвестировании используются для выбора оптимального портфеля;
- в банках и корпорациях — для расчета рисковых ситуаций.
ДУ применяются во многих областях, поэтому очень важно уметь их решать.
Современные языки программирования, такие как Python, имеют разные библиотеки, которые позволяют использовать многочисленные методы для решения ДУ.
Рассмотрим библиотеки и простые примеры их применения для решения разных ДУ. Но сначала кратко освежим теорию.
Дифференциальные уравнения — область математики, посвященная способам решения уравнений, в которых в качестве неизвестной выступает функция и ее производные различных порядков одного или нескольких аргументов.
Решить ДУ — найти такую функцию, которая при подстановке ее в уравнение обращает его в тождество.
Решение ДУ может быть:
- общим, когда решение ДУ описывает все множество решений;
- частным, когда решение имеет конкретные числовые значения.
Существует множество видов ДУ:
- обыкновенные,
- линейные,
- нелинейные,
- однородные,
- неоднородные,
- первого и высших порядков,
- в частных производных и пр.
Рассмотрим, как решать некоторые из видов ДУ с помощью Python.
Обыкновенные дифференциальные уравнения (ОДУ)
Они содержат функцию от одной независимой переменной. ОДУ первого порядка в общем виде можно записать так:
y'(x) = f(x)g(x),
где fx и g(x) — функции от переменной x.
Приведем простой пример такого уравнения:
y' = -2y
А также его решение с помощью библиотеки SciPy Python, разработанной для численного интегрирования и решения обыкновенных ДУ в частных производных.
Решить это уравнение можно следующим образом:
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
# Определение функции, представляющей ОДУ
def model(y, t):
dydt = -2 * y
return dydt
# Начальные условия
y0 = 1
# Время
t = np.linspace(0, 5, 100)
# Решение ОДУ методом odeint
solution_odeint = odeint(model, y0, t)
# визуализация решения
plt.plot(t, solution_odeint, label='y(t)')
plt.xlabel('Time')
plt.ylabel('y')
plt.legend()
plt.show()В этом примере использовались библиотеки:
- SciPy — для численного решения ДУ с помощью функции odeint, позволяющей интегрировать систему ОДУ;
- Numpy — для генерации массива элементов, соответствующих моментам времени, для которых вычисляется соответствующее решение;
- Matplotlib — для визуализации полученного решения.
Частные дифференциальные уравнения
Дифференциальным уравнением с частными производными (УЧП) называется уравнение относительно неизвестной функции нескольких переменных и ее частных производных, например:
dU/dt = -U
Решить это уравнение с помощью SciPy можно следующим образом:
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# Определение функции, представляющей УЧП
def model(t, u):
dudt = np.zeros_like(u)
dudt[0] = -u[0]
return dudt
# Начальные условия
u0 = np.array([1.0])
# Время
t = (0, 5)
# Решение УЧП методом solve_ivp
solution_ivp = solve_ivp(model, t, u0, method='RK45', dense_output=True)
# Визуализация решения
plt.plot(solution_ivp.t, solution_ivp.y[0], label='y(t)')
plt.xlabel('Time')
plt.ylabel('y')
plt.legend()
plt.show()В этом примере использовалась функция solve_ivp, в которую были переданы следующие параметры:
- model — уравнение, которое необходимо решить;
- t — интервал интегрирования;
- u0 — начальное условие;
- method — метод, с помощью которого велось решение (есть разные методы, «RK45» выбран как наиболее популярный: его подробнее рассмотрим ниже);
- dense_output — параметр, определяющий, вычислять или нет непрерывное решение.
Аналитическое решение дифференциальных уравнений с помощью SymPy
Библиотека SciPy используется для численного решения ДУ. А если необходимо аналитическое? Для этого можно воспользоваться библиотекой SymPy.
SymPy работает с символьными математическими выражениями, поэтому функции библиотеки позволяют осуществлять вычисления с большей точностью.
Программа, решающая ДУ с помощью SymPy, в ответе дает не набор чисел, а выражение, записанное математическими символами.
Для демонстрации возможностей библиотеки решим ДУ из первого примера.
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
# Определяем переменные
t = sp.symbols('t')
y = sp.Function('y')
# Определяем уравнения
equation = sp.Eq(y(t).diff(t), -2 * y(t))
# Его решение
solution = sp.dsolve(equation)
print(solution)На этом можно закончить работу программы, так как в выводе появится ответ
Eq(y(t), C1*exp(-2*t)),
точно описывающий решение ДУ.
Теперь определим численные значения полученной функции на том же интервале и построим график, чтобы сравнить с ответами в предыдущих примерах.
# Зададим те же значения времени, что и в первом примере
t1=np.linspace(0, 5, 100)
# Вычислим значения
y1=np.exp(-2*t1)
plt.plot(t1, y1, label='y(t)=C1*exp(-2*t)')
plt.xlabel('Time')
plt.ylabel('y')
plt.legend()
plt.show()Сравнив полученные графики в этом примере и в первом, легко убедиться, что результаты совпадают, но при использовании SymPy выводится искомая формула.
Численное решение дифференциальных уравнений методом Рунге-Кутты
Методы Рунге-Кутты — семейство алгоритмов, которые по шагам численно дают математическое приближение в качестве решения ОДУ.
Для демонстрации метода решим то же уравнение, что и в первом примере.
Для каждого из методов подготовим свою функцию, реализующую соответствующий названию алгоритм:
- метод Рунге-Кутты 1-го порядка (или метод Эйлера) реализуем с помощью функции rungekutta1;
- 2-го порядка — rungekutta2;
- 4-го или «RK45», который широко распространен, — rungekutta4.
Соответствующий код:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
def model(y, t):
dydt = -2 * y
return dydt
def rungekutta1(f, y0, t):
n = len(t)
y = np.zeros((n, 1))
y[0] = y0
for i in range(n - 1):
y[i+1] = y[i] + (t[i+1] - t[i]) * f(y[i], t[i])
return y
def rungekutta2(f, y0, t):
n = len(t)
y = np.zeros((n, 1))
y[0] = y0
for i in range(n - 1):
h = t[i+1] - t[i]
y[i+1] = y[i] + h * f(y[i] + f(y[i], t[i]) * h / 2., t[i] + h / 2.)
return y
def rungekutta4(f, y0, t):
n = len(t)
y = np.zeros((n, 1))
y[0] = y0
for i in range(n - 1):
h = t[i+1] - t[i]
k1 = f(y[i], t[i])
k2 = f(y[i] + k1 * h / 2., t[i] + h / 2.)
k3 = f(y[i] + k2 * h / 2., t[i] + h / 2.)
k4 = f(y[i] + k3 * h, t[i] + h)
y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)
return y
# Начальные условия
y0 = 1
# Время
t = np.linspace(0, 5, 100)
sol = odeint(model, y0, t)
sol1 = rungekutta1(model, y0, t)
sol2 = rungekutta2(model, y0, t)
sol4 = rungekutta4(model, y0, t)
plt.scatter(t, sol, label='odient',c='k')
plt.plot(t, sol1, label='rungekutta1')
plt.plot(t, sol2, label='rungekutta2')
plt.plot(t, sol4, label='rungekutta4')
plt.xlabel('Time')
plt.ylabel('y')
plt.legend()
plt.show()Реализовав этот код, можно увидеть, как отличаются решения разными методами и насколько они дают близкое решение к odeint библиотеки SciPy.
Решение систем дифференциальных уравнений
Рассмотрим следующую задачу:
Решить систему ОДУ с помощью SciPy:
y'1=2y2+x,
y'2=2y1+3sin(x)
y1(0)=0
y2(0)=1.
Решение этой задачи приведем в несколько этапов:
- Вызовем необходимые библиотеки.
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt- Создадим функцию, выразив производную из уравнения. Подобный метод применялся в примере ранее. Отличие в том, что переменная Y содержит оба значения — y1 и y2. Также созданная функция вернет сразу два соответствующих уравнениям значения.
def model(Y, x):
y1, y2 = Y
return [2 * y2 + x,
2 * y1 + 3 * np.sin(x)]- Определим начальные условия, исходя из задачи.
y1_0 = 0
y2_0 = 1- Зададим интервал, на котором будет вычисляться решение.
x = np.linspace(0, 1, 50)- Определим переменную sol, которая будет соответствовать решению с помощью функции odeint. В эту функцию необходимо передать следующие параметры:
- функцию, соответствующую уравнениям системы;
- начальные условия;
- интервал интегрирования, определенный ранее.
sol = odeint(model, (y1_0, y2_0), x)- Передадим решения в соответствующие переменные. Переменная sol в этом случае — матрица, где первый столбец соответствует решению первого уравнения, второй — второго.
y1 = sol[:,0]
y2 = sol[:,1]- Визуализируем полученное решение.
plt.plot(x, y1, label='y1')
plt.plot(x,y2, label='y2')
plt.xlabel("Time")
plt.ylabel("Y")
plt.legend()
plt.show()Преимущества использования Python и рассмотренных библиотек
- Простота и читаемость кода на Python позволяет даже неспециалисту в языке подготовить небольшую программу для решения ДУ.
- Также написание такого скрипта сильно упрощает использование библиотеки SciPy, которая имеет специальный инструментарий для решения ДУ.
- Инструменты SciPy содержат широкий набор методов решения ОДУ и УЧП.
- Применение функций SciPy существенно снижает объем кода, который нужно подготовить, и помогает получить решение в удобном для дальнейшего использования или импорта виде.
- Достоинства SciPy сделали библиотеку популярной среди ученых и инженеров, а значит, подготовленный с ее использованием код будет понятен и привычен для таких специалистов.
Преимущества использования Python и SymPy:
- Библиотека SymPy разработана полностью на Python, не требует сложной установки и легко интегрируется в разработанные ранее программы.
- Как и SciPy, библиотека доступна бесплатно, что существенно упрощает ее использование во многих проектах.
- SymPy не занимает много памяти компьютера, так как зависит только от библиотеки mpmath, а не большого числа разных пакетов.