Численное решение обыкновенных дифференциальных уравнений (ОДУ) в Python
Модуль scipy.integrate имеет две функции ode() и odeint(), которые предназначены для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (т.е. задача Коши).
Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.
from scipy.integrate import odeint Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат
Решение одного ОДУ
Допустим надо решить диф. уравнение 1-го порядка
import numpy as np from scipy. integrate import odeint import matplotlib.pyplot as plt # create function def dydt(y, t): return -y*t t = np.linspace( -2, 2, 51) # vector of time y0 = 1 # start value y = odeint (dydt, y0, t) # solve eq. y = np.array(y).flatten() plt.plot( t, y,'-sr', linewidth=3) # graphic plt.show() # display Решение системы ОДУ
Пусть теперь мы хотим решить (автономную) систему диф. уравнений 1-го порядка
import numpy as np from scipy. integrate import odeint import matplotlib.pyplot as plt # create function def f(y, t): y1, y2 = y return [y2, - y2 - y1] t = np.linspace( 0, 10, 41) # vector of time y0 = [0, 1] # start value w = odeint(f, y0, t) # solve eq. y1 = w[:,0] y2 = w[:,1] fig = plt.figure(facecolor='white') plt.plot(t, y1, '-o', t, y2, '-o', linewidth=2) plt.ylabel("z") plt.xlabel("t") plt.grid(True) plt.show() # display Выходной массив w состоит из двух столбцов — y1(t) и y2(t).
Также без труда можно построить фазовые траектории:
fig2 = plt.figure(facecolor='white') plt.plot(y1, y2, linewidth=2) plt.grid(True) plt.show() Платы ARDUINO
Now 24.07.23 17:11:55, Your IP: 95.143.190.109; spyphy.zl3p.com/python/34_scipy_ode
Решение ОДУ#
Подмодуль scipy.integrate содержит в себе множество методов решения систем обыкновенных дифференциальных уравнений (ОДУ, ordinary differential equation , ODE ).
Т.к. везде ниже будут строиться графики полученного и точного решений уравнения, определим функцию для построения графиков этих решений.
from matplotlib import pyplot as plt def plot_solution(ax, x, exact_sol, sol_x, sol_y, xlabel="$t$"): ax.plot(x, exact_sol(x), color="green", label="Точное решение") ax.scatter(sol_x, sol_y, color="red", marker=".", label="Приближенное решение") ax.set_xlabel(xlabel) ax.set_ylabel("$y$") ax.legend()
Задача Коши#
Функция scipy.integrate.solve_ivp (solve initial value problem) позволяет численно решать задачу Коши вида
Уравнение первого порядка#
В качестве самого простого примера решим уравнение
точное решение которого \(y=e^t\) .
Чтобы численно решить это уравнение средствами SciPy , необходимо определить функцию правой части ( f(t, y) в ячейке ниже), задать начальное значение \(y(0)\) в виде списка из одного элемента ( y_0 в ячейке ниже), а также задать интервал независимой переменной, на которой необходимо решить дифференциальное уравнение.
import numpy as np from scipy import integrate from matplotlib import pyplot as plt def f(t, y): return y def exact_solution(t): return np.exp(t) y_0 = [1] t_0 = 0 t_final = 3 solution = integrate.solve_ivp(f, (t_0, t_final), y_0) print(solution)
message: 'The solver successfully reached the end of the integration interval.' nfev: 26 njev: 0 nlu: 0 sol: None status: 0 success: True t: array([0. , 0.10001999, 1.06609106, 2.30431769, 3. ]) t_events: None y: array([[ 1. , 1.10519301, 2.9040598 , 10.01740317, 20.08580546]]) y_events: None
Метод возвращает структуру со сравнительно большим количеством полей, из которых можно получить отладочную информацию о сходимости алгоритма, а также получить приближенную оценку \(\tilde(t)\) в наборе точек \(\
fig, ax = plt.subplots(figsize=(8, 8), layout="tight") x = np.linspace(t_0, t_final, 100) plot_solution(ax, x, exact_solution, solution.t, solution.y[0])
Хоть правая часть уравнения и не зависит явным образом от \(t\) , функция f(t, x) все равно объявляется с первым параметром t .
Для данного уравнения функция scipy.integrate.solve_ivp выдала решение, содержащее совсем небольшое количество точек, так как метод быстро сошелся. Параметром t_eval можно в явно виде задать сетку, в узлах которой требуется получить оценку точного решения.
solution = integrate.solve_ivp(f, (t_0, t_final), y_0, t_eval=np.linspace(t_0, t_final, 50)) fig, ax = plt.subplots(figsize=(8, 8), layout="tight") plot_solution(ax, x, exact_solution, solution.t, solution.y[0])
Уравнения высшего порядка#
Предполагается, что в общем случае \(y=y(t)\) является векторной и решается система обыкновенных дифференциальных уравнений.
\[\begin
Если имеется дифференциальное уравнение высшего порядка разрешенное относительно старшей переменной вида
то необходимо свести его к системе из \(n\) уравнение первого порядка.
В качестве примера решим уравнение
точное решение которого — \(z = \sin t\) . Чтобы его решить, необходимо свести уравнение к системе двух дифференциальных уравнений первого порядка. Введём обозначения \(y_1(t) = z(t)\) и \(y_2(t) = z'(t)\) , тогда систему выше можно переписать в виде
В таком случае функция правой части уравнения должна принимать на вход вектор \(y\) и возвращать вектор производны \(y’=f(t, y)\) в виде списка или массива NumPy .
def exact_solution(t): return np.sin(t) def f(t, y): return [ y[1], -y[0] ] y_0 = [0, 1] t_0 = 0 t_final = np.pi t_eval = np.linspace(0, np.pi, 10) solution = integrate.solve_ivp(f, (t_0, t_final), y_0, t_eval=t_eval) print(solution)
message: 'The solver successfully reached the end of the integration interval.' nfev: 44 njev: 0 nlu: 0 sol: None status: 0 success: True t: array([0. , 0.34906585, 0.6981317 , 1.04719755, 1.3962634 , 1.74532925, 2.0943951 , 2.44346095, 2.7925268 , 3.14159265]) t_events: None y: array([[ 0.00000000e+00, 3.42054448e-01, 6.42796558e-01, 8.66269432e-01, 9.85009559e-01, 9.84744759e-01, 8.66190242e-01, 6.42630108e-01, 3.41454338e-01, -5.87748431e-04], [ 1.00000000e+00, 9.39714990e-01, 7.66028577e-01, 5.00009935e-01, 1.73634556e-01, -1.73906214e-01, -5.00472507e-01, -7.66382176e-01, -9.39862784e-01, -9.99961909e-01]]) y_events: None
В решении \(y\) тоже возвращается в виде массива из \(n\) строк и \(m\) столбцов, где \(n\) — порядок системы уравнений, а \(m\) — количество точек, в которых построена оценка решения. В данном примере первая строка массива соответствует искомой зависимости \(y_1(t)=z(t)\) , а вторая строка — \(y_2(t)\) , которая соответствует производной решения \(z'(t)\) .
x = np.linspace(0, np.pi, 50) fig, ax = plt.subplots(figsize=(8, 8), layout="tight") plot_solution(ax, x, exact_solution, solution.t, solution.y[0])
Выбор метода решения ОДУ#
Решение таких систем ОДУ отнюдь не тривиально. Разработано множество методов их решения и ряд из них “зашит” в подмодуле scipy.integrate , среди которых:
- Явные ме́тоды Ру́нге — Ку́тты 2-го ( RK23 ), 4-го ( RK45 ) и 8-го ( DOP853 ) порядков;
- Неявный метод Руне — Кутты 5-го порядка Radau ;
- Неявные методы BDF и LSODA .
Основным критерием выбора является жесткость системы ОДУ: явные методы плохо проявляют себя на жестких системах. SciPy рекомендует использовать по умолчанию метод RK45 и если он плохо/долго сходится переключаться на метод Radau иди BDF .
Краевая задача#
Функция scipy.integrate.solve_bvp (solve boundary value problem`) предназначена для решения системы ОДУ с краевыми условиями
В качестве аргументов функция solve_bvp принимает:
- функцию f , задающую правую часть уравнения,
- функцию bc (boundary condition), задающую невязку для граничных условий,
- массив x , определяющий сетку значений независимой переменной \(x\) ,
- массив y , задающий “догадку” об итоговом решении \(y = y(x)\) .
Функция bc должна принимать на вход \(y(a)\) и \(y(b)\) (векторные величины в общем случае) и возвращать вектор невязки на граничных условиях.
В качестве примера решим последнюю систему,
но с граничными условиями на отрезке \(x \in [0, \frac<\pi>]\) . Как и в случае задачи Коши здесь потребуется свести уравнение второго порядка к системе из двух уравнений первого порядка.
import numpy as np from scipy import integrate from matplotlib import pyplot as plt def boundary_residual(ya, yb): return np.array([ ya[0] - 0, yb[0] - 1 ]) a, b = 0, np.pi/2 N = 30 x = np.linspace(a, b, N) y_guess = np.zeros((2, N), dtype=float) sol = integrate.solve_bvp(f, boundary_residual, x, y_guess) fig, ax = plt.subplots(figsize=(8, 8), layout="tight") plot_solution(ax, x, exact_solution, sol.x, sol.y[0])
Интерполяция и аппроксимация