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

Я пытаюсь построить вывод 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)

python delta odeint
person Ponder    schedule 10.10.2019    source источник

Ответы (2)


arrow_upward
-1
arrow_downward

Я думаю, что проблема может быть во внутренних ошибках округления, потому что 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

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

person BurningKarl    schedule 10.10.2019
comment
Спасибо, что нашли время, чтобы изучить это! Очень подробный ответ :) - person Ponder; 11.10.2019
comment
Если мой ответ решил вашу проблему, вы можете принять его, чтобы повысить мою репутацию. Узнайте больше об этом здесь - person BurningKarl; 11.10.2019

arrow_upward
0
arrow_downward

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

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. .

person Matthew Flamm    schedule 11.10.2019