Python - Использование дельты Кронекера с ODEINT - PullRequest
Новая
       32

Python - Использование дельты Кронекера с ODEINT

1 голос
/ 10 октября 2019

Я пытаюсь построить вывод из ODE, используя дельта-функцию Кронекера, которая должна стать «активной» только в определенное время = t1. Это должно давать пилообразный отклик, когда начальное значение экспоненциально уменьшается до t = t1, где оно снова мгновенно возрастает, а затем снова снижается. Однако, когда я рисую это, похоже, что решатель видит дельта-функцию Кронекера как ноль за все время t. Есть ли способ сделать это в Python? 

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

def dy_dt(y,t):

    dy_dt = 500*KroneckerDelta(t,t1) - 2y

return dy_dt

t1 = 4
y0 = 500
t = np.arrange(0,10,0.1)

y = sp.odeint(dy_dt,y0,t)

plt.plot(t,y)

Ответы [ 2 ]

0 голосов
/ 11 октября 2019

В случае простой дельты Кронекера, использующей время, вы можете запустить оду по частям следующим образом: 

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

def dy_dt(y,t):
    return -2*y

t_delta = 4
tend = 10
y0 = [500]

t1 = np.linspace(0,t_delta,50)
y1 = odeint(dy_dt,y0,t1)

y0 = y1[-1] + 500  # execute Kronecker delta
t2 = np.linspace(t_delta,tend,50)
y2 = odeint(dy_dt,y0,t2)

t = np.append(t1, t2)
y = np.append(y1, y2)
plt.plot(t,y)

Другой вариант для сложных ситуаций - это events функциональность solve_ivp.

0 голосов
/ 10 октября 2019

Я думаю, что проблема может быть из-за внутренних ошибок округления, потому что 0.1 нельзя представить точно как питон float. Я бы попробовал 

import math

def dy_dt(y,t): 
    if math.isclose(t, t1):
        return 500 - 2*y
    else:
        return -2y

Также в документации odeint предлагается использовать параметр args вместо глобальных переменных, чтобы дать вашей производной функции доступ к дополнительным аргументам и заменить np.arange на np.linspace: 

import scipy.integrate as sp 
import matplotlib.pyplot as plt 
import numpy as np 
import math

def dy_dt(y, t, t1): 
    if math.isclose(t, t1):
        return 500 - 2*y
    else:
        return -2*y

t1 = 4 
y0 = 500 
t = np.linspace(0, 10, num=101) 
y = sp.odeint(dy_dt, y0, t, args=(t1,)) 
plt.plot(t, y)

Я не тестировал код, поэтому скажите, если с ним что-то не так.

РЕДАКТИРОВАТЬ:

При тестированиимой код Я посмотрел на t значения, для которых dy_dt оценивается. Я заметил, что odeint не только использует значения t, которые указаны, но и слегка их изменяет: 

...
3.6636447422787928
3.743098503914526
3.822552265550259
3.902006027185992
3.991829287543431
4.08165254790087
4.171475808258308
...

Теперь, используя мой метод, мы получаем 

math.isclose(3.991829287543431, 4) # False

, потому чтодопуск по умолчанию установлен на относительную ошибку не более 10 ^ (- 9), поэтому функция odeint «пропускает» удар производной на 4. К счастью, мы можем исправить это, указав более высокий порог ошибки: 

def dy_dt(y, t, t1): 
    if math.isclose(t, t1, abs_tol=0.01):
        return 500 - 2*y
    else:
        return -2*y

Сейчас dy_dt очень высоко для всех значений от 3,99 до 4,01. Можно уменьшить этот диапазон, если увеличить num аргумент linspace.

TL; DR

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

...