Почему Dfun (градиент) не вызывается при использовании интегрировать.odeint в SciPy?

Может ли кто-нибудь привести пример предоставления якобиана функции integrate.odeint в SciPy? Я пытаюсь запустить этот код из учебника SciPy пример odeint, но кажется, что Dfun() (функция Якоби) никогда не вызывается.

from numpy import * # added
from scipy.integrate import odeint
from scipy.special import gamma, airy
y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
y0 = [y0_0, y1_0]


def func(y, t):
    return [t*y[1],y[0]]
    
    
def gradient(y,t):
    print 'jacobian'  # added
    return [[0,t],[1,0]]
    
    
x = arange(0,4.0, 0.01)
t = x
ychk = airy(x)[0]
y = odeint(func, y0, t)
y2 = odeint(func, y0, t, Dfun=gradient)
print y2 # added

python scipy odeint
person lumartor    schedule 08.07.2013    source источник

Ответы (1)


arrow_upward
17
arrow_downward

Под капотом scipy.integrate.odeint использует решатель LSODA из библиотеки ODEPACK FORTRAN. Чтобы иметь дело с ситуациями, когда функция, которую вы пытаетесь интегрировать, является жесткой, LSODA адаптивно переключается между двумя разными методами вычисления интеграла — методом Адамса, который быстрее, но не подходит. для жестких систем и BDF, который медленнее, но устойчив к жесткости.

Конкретная функция, которую вы пытаетесь интегрировать, не является жесткой, поэтому LSODA будет использовать Адамса на каждой итерации. Вы можете проверить это, вернув infodict (...,full_output=True) и проверив infodict['mused'].

Поскольку метод Адамса не использует якобиан, ваша функция градиента никогда не вызывается. Однако если вы зададите odeint жесткую функцию для интегрирования, такую ​​как уравнение Ван-дер-Поля:

def vanderpol(y, t, mu=1000.):
    return [y[1], mu*(1. - y[0]**2)*y[1] - y[0]]

def vanderpol_jac(y, t, mu=1000.):
    return [[0, 1], [-2*y[0]*y[1]*mu - 1, mu*(1 - y[0]**2)]]

y0 = [2, 0]
t = arange(0, 5000, 1)
y,info = odeint(vanderpol, y0, t, Dfun=vanderpol_jac, full_output=True)

print info['mused'] # method used (1=adams, 2=bdf)
print info['nje']   # cumulative number of jacobian evaluations
plot(t, y[:,0])

вы должны увидеть, что odeint переключается на использование BDF, и теперь вызывается функция Якоби.

Если вам нужен больший контроль над решателем, вам следует изучить scipy.integrate.ode, который представляет собой гораздо более гибкий объектно-ориентированный интерфейс для нескольких различных интеграторов.

person ali_m    schedule 09.07.2013
comment
Идеально. Более понятная, чем документация SciPy. Это именно то, что я хочу, спасибо ali_m! - person lumartor; 09.07.2013