Scipy odeint оценивает при отрицательных значениях времени - PullRequest
Новая
       10

Scipy odeint оценивает при отрицательных значениях времени

2 голосов
/ 18 июня 2019

Я пытаюсь найти систему ODE с scipy.integrate.odeint из граничного условия в последний раз и работаю до первоначального времени (как описано здесь: Обратная интеграция во времени с использованием scipy odeint ).

Однако odeint переходит к отрицательным значениям времени - вне фактического диапазона решений, которые я ищу - и это вызывает ошибку, потому что мои ODE зависят от квадратного корня времени, а моя функция вместо этого возвращает комплексное значение действительного числа.

Вот пример, который воспроизводит проблему: 

    import numpy as np 
    from scipy.integrate import odeint
    tmax = 4e4
    tmin = 1
    t = np.linspace(tmax,tmin,1e3)
    param0 = [1] #value of x at tmax

    def func(param,t):
        x = param[0]
        dxdt = 1e-10/np.sqrt(t)
        print(t) #show what values of t are tried by odeint
        return dxdt

    res = odeint(func,param0,t)

Очень быстро print(t) показывает отрицательные значения, а res заполняется nan.

Есть ли способ предотвратить переход odeint к отрицательным значениям? Почему он пытается значения вне моего входного массива t?

Для моего реального кода я нашел несколько способов избежать результатов nan, таких как добавление if t<0: t=0 в функцию (что не сработало бы здесь) или наложение очень маленьких максимальных временных шагов (hmax<tmin) но это делает вычисления намного длиннее.

Обратите внимание, что Scipy odeint Неотрицательное решение связано, но моя проблема немного другая: меня не касается отрицательное решение x, а отрицательный аргумент t.

1 Ответ

0 голосов
/ 18 июня 2019

Я бы рекомендовал использовать solve_ivp вместо старого odeint. Обратите внимание, что solve_ivp ожидает f (t, y) вместо f (y, t) в правой части вашей оды: 

import numpy as np 
from scipy.integrate import solve_ivp
tmax = int(4e4)
tmin = 1
t = np.linspace(tmax, tmin, int(1e3))
param0 = [1] #value of x at tmax

def func(t, param):
    dxdt = 1e-10/np.sqrt(t)
    return dxdt

res = solve_ivp(func, y0=param0, t_span=[tmax, tmin], t_eval=t)

, где t_eval=t гарантирует, что решение будет сохранено в моменты времени t. В противном случае решатель выбирает моменты времени.

...