Как оценить функцию во времена интеграции, используя Scipy solve_ivp

Question

Я использую метод execute_ivp от scipy.integrate для решения ivp, и я хочу иметь возможность оценивать функцию по временным шагам, которые я даю для интеграции, но я не знаю, как это сделать.

Я мог бы вернуться к каждому из элементов в интеграции, но это заняло бы смешное количество времени в дополнение к времени, которое уже требуется для решения IVP, так что я бы предпочел вычислитьих в то же время, что фактический метод вычисляет значения во время интеграции.

import scipy.integrate
import numpy

class Foo:
    def __init__(self):
        self.foo_vector_1 = numpy.zeros(3)
        self.foo_vector_2 = numpy.zeros(3)
        self.foo_vector_3 = numpy.zeros(3)

foo = Foo()

d_vector_1 = lambda foo: # gets the derivative of foo_vector_1
d_vector_2 = lambda foo: # gets the derivative of foo_vector_2

def get_foo_vector_3_value(foo):
    return # returns the ACTUAL VALUE of foo_vector_3, NOT its derivative

def dy(t, y):
    foo.foo_vector_1 = numpy.array((y[0],y[1],y[2]))
    foo.foo_vector_2 = numpy.array((y[3],y[4],y[5]))
    return numpy.array((d_vector_1(foo),d_vector_2(foo))).flatten().tolist()

foo.foo_vector_1 = numpy.array((1,2,3))
foo.foo_vector_2 = numpy.array((4,5,6))

y0 = numpy.array((foo.foo_vector_1, foo.foo_vector_2)).flatten().tolist()

sol = scipy.integrate.solve_ivp(dy, (0,10), y0, t_eval=numpy.arange(0,1000,1))

foo_vectors_1 = numpy.column_stack((sol.y[0], sol.y[1], sol.y[2]))
foo_vectors_2 = numpy.column_stack((sol.y[3], sol.y[4], sol.y[5]))
foo_vectors_3 = ????????

В идеале, я смог бы получить значение foo_vectors_3 без необходимости сбрасывать foo в цикле по всем спискамвекторов foo, потому что для меня это потребовало бы значительного времени вычислений.

Answer 1

Я думаю, что трение здесь позволяет избежать использования 1D numpy ndarray в качестве базового объекта для вычислений.Вы можете мысленно распределить массив 1D в ваши 2 отдельных атрибута foo. Тогда вычисление foo_vectors_3 будет тривиальным по сравнению с интеграцией ODE.Вы также можете добавить вспомогательные функции на карту из 1D ndarray для solve_ivp и вашего foo_vectors и обратно.

In [65]: import scipy.integrate 
    ...: import numpy as np 
    ...:  
    ...: def d_vec1(t, y): 
    ...:     # put in your function here instead of just returning 1 
    ...:     return 1 * np.ones_like(y) 
    ...:  
    ...: def d_vec2(t, y): 
    ...:     # put in your function here instead of just returning 2 
    ...:     return 2 * np.ones_like(y) 
    ...:  
    ...: def eval_foo3(t, y): 
    ...:     return y[0:3,:] + y[3:,:]  # use your own function instead 
    ...:  
    ...: def dy(t, y): 
    ...:     return numpy.array((d_vec1(t, y[0:3]), d_vec2(t, y[3:]))).flatten() 
    ...:  
    ...: v1 = np.array([1, 2, 3]) 
    ...: v2 = np.array([4, 5, 6]) 
    ...: y0 = numpy.array((v1, v2)).flatten() 
    ...: t_eval = np.linspace(0, 10, 11) 
    ...: sol = scipy.integrate.solve_ivp(dy, (0, 10), y0, t_eval=t_eval) 
    ...:  
    ...: foo3 = eval_foo3(sol.t, sol.y) 
    ...: print(sol.y[0:3]) 
    ...: print(sol.y[3:]) 
    ...: print(foo3)

[[ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10. 11.]
 [ 2.  3.  4.  5.  6.  7.  8.  9. 10. 11. 12.]
 [ 3.  4.  5.  6.  7.  8.  9. 10. 11. 12. 13.]]
[[ 4.  6.  8. 10. 12. 14. 16. 18. 20. 22. 24.]
 [ 5.  7.  9. 11. 13. 15. 17. 19. 21. 23. 25.]
 [ 6.  8. 10. 12. 14. 16. 18. 20. 22. 24. 26.]]
[[ 5.  8. 11. 14. 17. 20. 23. 26. 29. 32. 35.]
 [ 7. 10. 13. 16. 19. 22. 25. 28. 31. 34. 37.]
 [ 9. 12. 15. 18. 21. 24. 27. 30. 33. 36. 39.]]