Odeint python как работает
Перейти к содержимому

Odeint python как работает

  • автор:

Численное решение обыкновенных дифференциальных уравнений (ОДУ) в Python

Модуль scipy.integrate имеет две функции ode() и odeint(), которые предназначены для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (т.е. задача Коши).

Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

from scipy.integrate import odeint 

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат

odeint(func, y0, t[,args=(), . ]) 

Решение одного ОДУ

Допустим надо решить диф. уравнение 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 30.10.23 1:35:56, Your IP: 178.132.111.91; spyphy.zl3p.com/python/34_scipy_ode

Математика за компьютером

Если искомая функция является полиномом, для её дифференцирования и интегрирования можно использовать методы модуля numpy.polynomial.polynomial. Хочу обратить внимание, что самый левый коэффициент соответствует младшему члену полинома, правый — старшему. Функции polyder() и polyint() содержат дополнительные опциональные коэффициенты, которые позволяют использовать их более гибко.

Если данные заданы в виде массива, выполнить интегрирование можно с помощью метода trapz().

Наиболее универсальные способы интегрирования и дифференцирования предоставляют модули SymPy и MpMath. Последний устанавливается вместе с SymPy и служит для проведения вычислений произвольной точности. Он же выполняет численное дифференцирование функции. Результат сохраняется в виде объекта mpf, который хранит число с плавающей точкой.

Для вычисления интеграла можно воспользоваться методом integrate(), указав пределы интегрирования.

Как известно, не каждая функция может быть проинтегрирована аналитически, да и не всегда это нужно. Для численных расчётов пригодится модуль scipy.integrate, который содержит коллекцию разнообразных функций. В простейшем случае, для вычисления интеграла на заданном интервале подойдёт quad(). Результатом будут 2 числа, первое из которых — найденное значение, а второе — точность.

Для решения обыкновенных дифференциальных уравнений пригодится функция odeint() из этого же модуля. Она работает с уравнениями вида y’ = f(y,t) и принимает в качестве аргументов саму функцию f, а также начальное значение для y и интервал времени t (который должен быть массивом). Результатом является массив значений функции y в точках интервала t.

Моделирование движения тел в гравитационном поле на python (scipy.odeint, matplotlib, numpy). Как исправить?

Я моделирую движение вокруг Солнца с помощь python (scipy.odeint), строю графики с помощью matplotlib. Но мои решения не совпадают с учебником 6051052ae98f5582439581.pngФормулы, которые я использовал: 605105753e819095047053.png
Код:

from scipy.integrate import odeint import scipy.constants as constants import numpy as np import matplotlib.pyplot as plt M = 1.989 * (10 ** 30) G = constants.G alpha = 30 alpha0 = (alpha / 180) * np.pi v00 = 0.7 # km/s v0 = v00 * 1000 # m/s def dqdt(q, t): x = q[0] y = q[2] ax = - G * M * (x / ((x ** 2 + y ** 2) ** 1.5)) ay = - G * M * (y / ((x ** 2 + y ** 2) ** 1.5)) return [q[1], ax, q[3], ay] vx0 = v0 * np.cos(alpha0) vy0 = v0 * np.sin(alpha0) x0 = -1.5 * (10 ** 11) y0 = 0 * (10 ** 11) q0 = [x0, vx0, y0, vy0] N = 1000000 t = np.linspace(0.0, 100000000000.0, N) pos = odeint(dqdt, q0, t) x1 = pos[:, 0] y1 = pos[:, 2] plt.plot(x1, y1, 0, 0, 'ro') plt.ylabel('y') plt.xlabel('x') plt.grid(True) plt.show()

Результаты:
А: 6051062bda632524100789.png
Б: 605106660a511802129633.png
В: 605106922a4fc625762560.png
C: 605106cab31cf882297062.png
Как мне это исправить?
Может вы сможете подсказать другой метод решения, например, с помощью метода Эйлера или с использованием других библиотек. Буду очень благодарен за любую помощь)

  • Вопрос задан более двух лет назад
  • 490 просмотров

Digiratory

Лаборатория автоматизации и цифровой обработки сигналов

Оценка параметров ДУ в Python

Одной из распространенных задач, возникающих при исследовании различных объектов — построение математической модели. Нередко математическая модель представляется в виде системы дифференциальных уравнений, однако непосредственное измерение всех, входящих в них параметров, как правило, невозможно по различным причинам. В таком случае, одним из подходов является проведение идентификационных экспериментов и оценка параметров ДУ путем решения оптимизационной задачи.

В статье рассмотрен простой способ оценки параметров системы ДУ в форме Коши на языке Python.

Постановка задачи

Разработка простого модуля Python для решения задачи оценки параметров.
Исходные данные для задачи:

  • Представление модели в форме Коши, первые N величин — наблюдаемы (порядок уравнений не имеет значения с точки зрения математики, но упорядочивание сильно упростит жизнь при разработке);
  • Имеются данные эксперимента, на основе которого будет производится оценка параметров/

Исходные данные

Исходными данными, передаваемыми при создании экземпляра класса оценки параметров являются:

  1. Функция, реализующая вычисление ДУ.
  2. Массив экспериментальных данных (Два массива: время и значения)

Таким образом, получаем следующий класс с конструктором:

class parameter_estimator(): def __init__(self, x_data, y_data, f): self._x_data = x_data self._y_data = y_data self._f = f self._c = None self.n_observed = y_data.shape[1]

Тут просто сохраняются значения параметров и функция в поля класса.

Оценка параметров

Процедуру оценки параметров можно сделать на основе функции scipy.optimize.leastsq. Эта функция минимизирует сумму квадратов полученных величин и может принимать два интересных для начального понимания аргумента: func и x0.

Функция func должна принимать один аргумент (который м.б. вектором), являющийся параметрами, а x0 — начальные значения параметров.

Метод оценки параметров можно реализовать следующим образом:

def estimate(self, y0, guess): """ Произвести оценку параметров дифференциального уравнения с заданными начальными значениями параметров: y0 -- начальные условия ДУ guess -- параметры ДУ """ # Сохраняем число начальных условий self._y0_len = len(y0) # Создаем вектор оцениваемых параметров, # включающий в себя начальные условия self._est_values = np.concatenate((y0, guess)) # Решить оптимизационную задачу - решение в переменной c (c, kvg) = optimize.leastsq(self.f_resid, self._est_values) self._c = c # В возвращаемом значении разделяем начальные условия и параметры return c[self._y0_len:], c[0:self._y0_len]

Метод получает начальные значения начальных условий ДУ и начальные значения параметров системы. Функцию вычисления вектора ошибок между решением системы ДУ с текущими значениями параметров и экспериментальными значениями необходимо реализовывать самостоятельно, потому, следующим шагом будет ее написание:

def f_resid(self, p): """ Функция для передачи в optimize.leastsq При дальнейших вычислениях значения, возвращаемые этой функцией, будут возведены в квадрат и просуммированы """ delta = self._y_data - self.my_ls_func(self._x_data, p) return delta.flatten() # Преобразуем в одномерный массив

Функция получает вектор параметров системы (включая начальные условия ДУ).

Предложенная функция содержит вызов еще одной функции my_ls_func, выделенной для разделения функционала. Выглядит она следующим образом и вычисляет непосредственно решение системы ДУ при заданных параметрах:

def my_ls_func(self, x, teta): """ Определение функции, возвращающей значения решения ДУ в процессе оценки параметров x заданные (временные) точки, где известно решение (экспериментальные данные) teta -- массив с текущим значением оцениваемых параметров. Первые self._y0_len элементов -- начальные условия, остальные -- параметры ДУ """ # Для передачи функции создадим ламбда-выражение с подставленными # параметрами # Вычислим значения дифференциального уравнения в точках "x" r = integrate.odeint(lambda y, t: self._f(y, t, teta[self._y0_len:]), teta[0:self._y0_len], x) # Возвращаем только наблюдаемые переменные return r[:, 0:self.n_observed]

Таким образом, процедура оценки параметров работает следующим образом:

  1. Вычисляется решение системы ДУ при текущих (или начальных) значения параметров.
  2. Рассчитывается вектор ошибки между решением и экспериментальными данными.
  3. Производится корректировка значений параметров.
  4. Переход на шаг 1, если не достигнуто максимальное число итераций или минимальное рассогласование.

Визуализация

Для визуализации напишем следующую функцию:

def plot_result(self): """ Построить графическое представление результатов оценки параметров """ if self._c is None: print("Parameter is not estimated.") return sol, t = self.calcODE((self._c[self._y0_len:],), self._c[0:self._y0_len], min(self._x_data), max(self._x_data)) # Строим экспериментальные данные, как красные точки, # а результаты моделирования, как синюю линию pp.plot(self._x_data, self._y_data, '.r', t, sol, '-b') pp.xlabel('xlabel', ) pp.ylabel("ylabel", ) pp.legend(('data', 'fit'), loc=0) pp.show()

где calcODE — еще одна служебная функция для расчета системы ДУ:

def calcODE(self, args, y0, x0=0, xEnd=10, nt=101): """ Служебная функция для решения ДУ """ t = np.linspace(x0, xEnd, nt) sol = odeint(self._f, y0, t, args) return sol, t

Тестирование и пример использования

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

где \(\theta \) — угол, \(\omega\) — скорость.

Этот объект уже использовался в качестве иллюстрации к статье о фазовых портретах и частотному методу синтеза регуляторов.

Для получения «экспериментальных» данных воспользуемся решением системы ДУ с известными «истинными» параметрами:

def ode(y, t, k): """ Функция, реализующая систему ДУ маятника с трением """ theta, omega = y b = k[0] c = k[1] dydt = [omega, -b*omega - c*np.sin(theta)] return dydt def calcODE(args, y0, dy0, ts=10, nt=101): """ Вспомогательная функция для получения решения систему ДУ """ y0 = [y0, dy0] t = np.linspace(0, ts, nt) sol = odeint(ode, y0, t, args) return sol, t # Зададим истинные значения параметров системы b = 0.3 c = 5.0 args = ([b, c], ) y0 = 1 dy0 = 0 print("Real parameter: b = <>, c = <>".format(b, c)) print("Real initial condition: <> <>".format(y0, dy0)) sol, t = calcODE(args, y0, dy0)

Объявленная тут система будет использоваться и для оценки параметров, но уже без известных значений:

# Зададим истинные значения параметров системы b = 0.3 c = 5.0 args = ([b, c], ) y0 = 1 dy0 = 0 print("Real parameter: b = <>, c = <>".format(b, c)) print("Real initial condition: <> <>".format(y0, dy0)) sol, t = calcODE(args, y0, dy0) guess = [0.2, 0.3] # Начальные значения для параметров системы y0 = [0, 1] # Стартовые начальные значения для системы ДУ estimator = parameter_estimator(t, sol, ode) est_par = estimator.estimate(y0, guess) # Построим графики результатов оценки параметров estimator.plot_result() print("Estimated parameter: b=<>, c=<>".format(est_par[0][0], est_par[0][1])) print("Estimated initial condition: <>".format(est_par[1]))

Получившиеся результаты являются достаточно точно оценкой:

Real parameter: b = 0.3, c = 5.0
Real initial condition: 1 0
Estimated parameter: b = 0.299993833612628, c = 5.000066486504925
Estimated initial condition: [ 9.99976783e-01 1.58491793e-04]

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *