Можете ли вы привести простой пример адаптивной пошаговой функции scipy.integrate.LSODA? - PullRequest
Новая
       85

Можете ли вы привести простой пример адаптивной пошаговой функции scipy.integrate.LSODA?

1 голос
/ 08 июля 2019

Мне нужно понять механизм функции scipy.integrate.LSODA.

Я написал тестовый скрипт, который интегрирует простую функцию.Согласно веб-странице LSODA входами функций могут быть функции rhs, t_min, initial y и t_max.С другой стороны, когда я запускаю код, я ничего не получаю.Что мне делать? 

import scipy.integrate as integ
import numpy as np

def func(t,y):
    return t

y0=np.array([1])
t_min=1
t_max=10
N_max=100
t_min2=np.linspace(t_min,t_max,N_max)
first_step=0.01

solution=integ.LSODA(func,t_min,y0,t_max)
solution2=integ.odeint(func,y0,t_min2)

print(solution.t,solution.y,solution.nfev,'\n')
print(solution2)

Решение дать 

1 [ 1.] 0

[[  1.00000000e+00]
 [  9.48773662e+00]
 [  9.00171421e+01]
 [  8.54058901e+02]
 [  8.10308559e+03]]

1 Ответ

1 голос
/ 08 июля 2019

1.) Вы только инициируете экземпляр класса решателя LSODA, никаких вычислений не происходит, только инициализация массивов с начальными данными. Чтобы получить odeint -подобный интерфейс, используйте solve_ivp с опцией method='LSODA'.

2.) Без опции tfirst=True решатель LSODA будет решать y'(t)=t, а odeint будет решать y'(t)=y(t)

Чтобы получить сопоставимые результаты, необходимо также выровнять допуски, поскольку допуски по умолчанию могут быть разными. Таким образом, можно вызывать такие методы, как 

print "LSODA"
solution=integ.solve_ivp(func,[t_min, t_max],y0,method='LSODA', atol=1e-4, rtol=1e-6)
print "odeint"
solution2=integ.odeint(func,y0,t_min2, tfirst=True, atol=1e-4, rtol=1e-6)

Даже тогда вы не получите никакой информации о внутренних шагах odeint, даже если для кода FORTRAN есть опция для этого, оболочка Python не копирует ее. Вы можете добавить оператор печати в функцию ODE func, чтобы увидеть, в каких местах эта функция на самом деле вызывается, это должно составлять в среднем 2 вызова с аргументами close-by на каждый внутренний шаг.

Это сообщает 

LSODA
1.0 [1.]
1.00995134265 [1.00995134]
1.00995134265 [1.01005037]
1.01990268529 [1.02010074]
1.01990268529 [1.02019977]
10.0 [50.50009903]
10.0 [50.50009903]
odeint
1.0 [1.]
1.00109084546 [1.00109085]
1.00109084546 [1.00109204]
1.00218169092 [1.00218407]
1.00218169092 [1.00218526]
11.9106363102 [71.43162985]

где сообщаемые шаги в выводе LSODA равны 

[ 1.          1.00995134  1.01990269 10.        ] [[ 1.          1.01005037  1.02019977 50.50009903]] 7

Конечно, метод высокого порядка объединит линейный полином y'=t с квадратичным полиномом y(t)=0.5*(t^2+1), практически без ошибок, независимо от размера шага.

...