Метод Odeint() в Scipy Python

Одна или несколько переменных (неизвестных) и некоторые их производные составляют дифференциальное уравнение. Другими словами, производные переменных определяются дифференциальным уравнением.

Дифференциальные уравнения классифицируются в соответствии с их порядком. Старшая производная, иногда называемая дифференциальным коэффициентом, которая присутствует в дифференциальном уравнении, определяет его порядок.

В этом уроке мы будем решать систему обыкновенных дифференциальных уравнений, используя метод Python Scipy.

В Python Scipy есть метод odeint() в модуле scipy.integrate, который решает набор обыкновенных дифференциальных уравнений путем их интегрирования.

Синтаксис:

scipy.integrate.odeint(func full_output=0, ml=None, mu=None, , y0, t, args=(), Dfun=None, col_deriv=0, rtol=None, atol=None, tcrit=None, h0=0.0, mxhnil=0, mxordn=12, mxords=5, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, printmessg=0, tfirst=False)

Где параметры:

• func (callable): вычисление производной y в момент времени t. Аргумент tfirst должен иметь значение True, если подпись является вызываемой(t, y,…)
• y0 (массив): начальное состояние y(может быть вектором).
• t (массив): список временных интервалов, для которых необходимо решить y. Первым компонентом этой последовательности должна быть начальная точка значения. Повторяющиеся значения допускаются; эта последовательность должна быть монотонно возрастающей или убывающей.
• args (tuple): дополнительные параметры для передачи функции.
• dfun: градиент функции(якобиан). Параметр tfirst должен быть установлен в значение True, если подпись является вызываемой(t, y,…)
• col_deriv (boolean): True, если Dfun определяет производные по строкам, а не по столбцам(что занимает больше времени).
• full_output (boolean): True, если второй вывод должен быть словарем необязательных выходных данных.
• printmessg (boolean): следует ли печатать сообщение о конвергенции.
• tfirst (boolean): В этом случае первые два аргумента func(и Dfun, если он указан) должны быть t, y, а не обычные y, t.

Метод odeint() возвращает y (массив с начальным значением y0 в первой строке и значением y для каждого выбранного момента времени в t) типа array.

Давайте рассмотрим пример, выполнив следующие шаги:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

• Вектор [тета, омега] должен быть y. Мы используем Python для реализации этой системы.

def pend_fun(y_, t_, b_, c_):
    theta_, omega_ = y_
    dydt_ = [omega_, -b_*omega_ - c_*np.sin(theta_)]
    return dydt_

• Мы предполагаем, что b = 0,20 и c = 4,5 являются постоянными.

b_ = 0.20
c_ = 4.5

• Первоначально маятник покоится с омега(0) = 0, а начальные обстоятельства таковы, что он примерно вертикальен с тета(0) = пи – 0,1. Вектор начальных условий тогда.

y0_ = [np.pi - 0.1, 0.0]

• В диапазоне 0 = t = 10 мы получим решение на 101 равномерно расположенной выборке. Итак, вот наш диапазон времени.

t_ = np.linspace(0, 10, 101)

• Чтобы сгенерировать решение, вызовите odeint. Мы используем аргумент args для передачи аргументов b и c в odeint, чтобы их можно было передать в функцию ожидания.

sol_ =integrate.odeint(pend_fun, y0_, t_, args=(b_, c_))

Ответ — массив, сохраняющий форму(101, 2). Тета(t) — это первый столбец, а Омега — второй(t). Код ниже отображает оба элемента.

plt.plot(t_, sol_[:, 0], 'b', label='theta(t)')
plt.plot(t_, sol_[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

• Сравнение Odeint и Solve_ivp
• Odeint 2-го порядка
• Метод complex_ode()
• Rk45

Сравнение Odeint и Solve_ivp

Класс scipy.integrate.ode и функция scipy.integrate.solve_ivp используют функцию определения системы, которая по умолчанию требует, чтобы первые два параметра func были в порядке, противоположном этим аргументам. Необходимо установить для аргумента tfirst значение True, чтобы использовать функцию с сигнатурой func(t, y,…).

Методsolve_ivp() устраняет проблему с начальным значением в системе ODE. Учитывая начальное значение, эта функция численно интегрирует набор обыкновенных дифференциальных уравнений.

dy / dt = f(t, y)
y(t0) = y0

Здесь y(t) — векторная функция ND(состояние), t — одномерная независимая переменная(время), а f(t, y) — векторная функция ND, которая устанавливает дифференциальные уравнения. Учитывая начальное значение y(t0)=y0, цель состоит в том, чтобы найти y(t), примерно удовлетворяющую дифференциальным уравнениям.

Правая часть жестких решателей ОДУ должна быть комплексно-дифференцируемой, однако некоторые решатели предлагают интегрирование в комплексной области(удовлетворяют уравнениям Коши-Римана [11]).

Передайте y0 со сложным типом данных, чтобы решить сложную проблему. Переписать задачу индивидуально для реального и воображаемого разделов — еще одна альтернатива, которая всегда доступна.

Синтаксис:

scipy.integrate.solve_ivp(fun, dense_output=False, events=None,  t_span, y0, method='RK45', t_eval=None,vectorized=False, args=None, **options)

Где параметры:

• fun (callable): Система находится на правой стороне. Забава с вызывающей сигнатурой(т,у). Здесь t — скаляр, а ndarray y имеет две возможности: альтернативно, если он имеет форму(n), fun должна вернуть объект, подобный массиву, с формой(n). Он также может иметь форму(n, k), и в этом случае fun должна вернуть массив, форма которого(n, k) такова, что каждый столбец эквивалентен одному столбцу в y. Векторизованный аргумент определяет, какая из двух альтернатив выбрана. Векторизованная реализация позволяет быстрее получить якобиан посредством аппроксимации конечными разностями(требуется для жестких решателей).
• _span (два кортежа с плавающей запятой): представляет интервал интегрирования от t0 до t1.
• y0 (array_data, например shape(n)): используется для указания начального состояния.
• method (строка или odsolver): используется для указания того, какой метод интеграции использовать, например RK45, Radau, RK23, DOP853, BDF и LSODA.
• t_eval (array_data): используется для указания времени сохранения вычисленного решения.
• Density_output (boolean): используется для указания, хотим ли мы вычислить непрерывное решение или нет.
• event (callable): используется для отслеживания событий.
• Vectorized (boolean): используется для указания, хотим ли мы реализовать функцию векторизованным способом или нет.
• args (tuple): для передачи дополнительных аргументов.
• **options: параметры, предоставленные выбранному решателю. Ниже приведен список всех опций, доступных для ранее реализованных решателей.
• first_step (None, float): размер первого шага. Поскольку значение по умолчанию — «Нет», решение должен принять алгоритм.
• max_step (float): Максимально допустимый размер шага. Размер шага не ограничен и полностью зависит от решателя, если np.inf используется по умолчанию.
• rtol, atol ( array_data, float): допуски, как относительные, так и абсолютные. Решение поддерживает локальные оценки ошибок, которые меньше, чем atol плюс rtol, умноженный на абс(y). Здесь абсолютная точность контролируется atol, а относительная точность (количество точных цифр) контролируется rtol (количество правильных десятичных знаков). Установите atol меньшим, чем наименьшее значение, которое можно предсказать из rtol * abs(y), чтобы rtol доминировал над допустимой ошибкой и чтобы получить желаемый rtol. Нет никакой гарантии, что будет столько же точных чисел, если atol больше, чем rtol * abs(y). С другой стороны, чтобы получить требуемый атол, отрегулируйте rtol так, чтобы оно никогда не превышало atol и rtol * abs(y). Может быть выгодно передать массив с формой(n) для atol, если масштабы компонентов y отличаются друг от друга.
• jac (array_data, разреженная матрица): правая часть матрицы Якоби системы относительно y — это то, что требуется для методов Радау, BDF и LSODA.(n, n)-мерная матрица Якобиана имеет элемент I j), равный d f_i/d y_j. Якобиан можно описать тремя разными способами.
  1. Якобиан считается постоянным, если тип матрицы или массива подобен массиву. Не поддерживается LSODA.
  2. Якобиан при необходимости будет называться jac(t, y), если он доступен для вызова, и в этом случае предполагается, что он зависит как от t, так и от y. Возвращаемое значение для алгоритмов Радау и BDF может быть разреженной матрицей.
  3. Если «Нет»(по умолчанию), для аппроксимации якобиана будут использоваться конечные разности.
• jac_sparsity (array_data, разреженная матрица): описывает структуру разреженности матрицы Якоби для аппроксимации конечной разностью. Его форма должна быть(n, n). Если jac не имеет значения None, этот аргумент не учитывается. Вычисления будут значительно ускорены, если обеспечить структуру разреженности и якобиан будет иметь лишь несколько ненулевых элементов в каждой строке [10]. Нулевая запись указывает, что связанный с якобианом элемент всегда равен нулю. Если None(по умолчанию), предполагается, что якобиан плотный. Не поддерживается «LSODA», вместо этого используйте lband и uband.
• lband, uband(int): Параметры метода LSODA, т.е. jac[i, j]!= 0 только для i – lband <= j <= i + uband, определяют полосу пропускания якобиана. По умолчанию не используется. Чтобы их установить, ваш метод jac должен возвращать якобиан в упакованном формате, что требует, чтобы возвращаемый массив имел n столбцов и строки uband + lband + 1 с записанными в них якобианскими диагоналями. Jac упакованный[uband + I – j, j] в частности равен Jac[i, j]. Scipy.linalg.solve Banded использует тот же формат(см. иллюстрацию). Чтобы уменьшить количество якобианских элементов, оцениваемых с помощью конечных разностей, эти параметры можно альтернативно использовать с jac=None.
• min_step (float): наименьший размер шага, разрешенный алгоритмом «LSODA». Значение по умолчанию для минимального шага — 0.

Метод solve_ivp() возвращает группу объектов.

Давайте возьмем пример и поймем, как работаетsolve_ivp(), выполнив следующие шаги:

• Импортируйте необходимые библиотеки или методы, используя приведенный ниже код Python.

from scipy import integrate

• Простое экспоненциальное затухание со случайно выбранными моментами времени.

def exponential_decay_fun(t_, y_): return -0.4 * y_
sol = integrate.solve_ivp(exponential_decay_fun, [0, 10], [1, 3, 7])
print(sol.t)

Odeint 2-го порядка

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

Мы собираемся использовать метод odeint() Python Scipy, описанный в приведенном выше подразделе.

Импортируйте необходимые библиотеки или методы, используя приведенный ниже код Python.

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

Определите функцию второго порядка и найдите ее производную, используя приведенный ниже код.

def fun(u_,x_):
  return(u_[1],-2*u_[1]-2*u_[0]+np.cos(2*x_))
y0_ =[0,0]
xs_ = np.linspace(1,10,300)
us_ = integrate.odeint(fun,y0_,xs_)
ys_ = us_[:,0]

Постройте вышеуказанную производную, используя приведенный ниже код.

plt.plot(xs_,ys_,'-')
plt.plot(xs_,ys_,'r*')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.show()

Метод complex_ode()

В Python Scipy есть метод complex_ode() в модуле scipy.integrate, который представляет собой обертку для сложных систем.

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

Синтаксис:

scipy.integrate.complex_ode(f, jac=None)

Где параметры:

• f (вызываемый): Rhs уравнения. Y.shape ==(n,) и t — скаляр. Установите f_args, используя set f params(*args).
• jac (callable): Якобиан правой руки «df[i] / dy[j] равен jac[i,j]. Вызов «set f params(*args)» устанавливает переменную «jac args». «.

Rk45

Метод Рунге-Кутты-Фельберга, иногда известный как метод Фельберга, представляет собой численный аналитический подход, используемый для решения обыкновенных дифференциальных уравнений. На основе широкого класса процедур Рунге-Кутты она была создана немецким математиком Эрвином Фельбергом.

Метод Фельберга уникален, поскольку это встроенный метод семейства Рунге-Кутты, что означает, что он объединяет множество оценок одной и той же функции для создания методов с разными порядками и одинаковыми константами ошибок.

Подход, представленный в публикации Фельберга 1969 года, известный как метод RKF45, представляет собой метод 4-го порядка с оценкой ошибки 5-го порядка. Встроенная методика высшего порядка, позволяющая автоматически определять размер адаптивного шага, может быть использована для оценки и регулирования погрешности решения путем проведения дополнительного расчета.

В этом разделе мы будем использовать описанный выше метод Rk45 с методомsolve_ivp() для решения дифференциальных уравнений, поскольку метод odeint() не поддерживает какой-либо метод параметров для Rk45.

Мы будем использовать тот же пример, который мы использовали в приведенном выше подразделе, выполнив следующие шаги:

Импортируйте необходимые библиотеки или методы, используя приведенный ниже код Python.

from scipy import integrate

Простое экспоненциальное затухание со случайно выбранными моментами времени и с помощью метода RK45 с использованием приведенного ниже кода.

def exponential_decay_fun(t_, y_): return -0.4 * y_
sol = integrate.solve_ivp(exponential_decay_fun, [0, 10], [1, 3, 7], method = 'RK45')
print(sol.t)