如何使用 Scipy 的 solve_ivp 在积分时评估函数答案

【问题标题】：How to Evaluate Function at Integration Times Using Scipy's solve_ivp如何使用 Scipy 的 solve_ivp 在积分时评估函数
【发布时间】：2019-06-18 05:07:31
【问题描述】：

我正在使用 scipy.integrate 的 solve_ivp 方法来解决 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，因为对我来说这实际上需要大量的计算时间。

【问题讨论】：

标签： python scipy

【解决方案1】：

我认为这里的摩擦是避免使用 1D numpy ndarray 作为计算的基础对象。您可以在心理上将一维数组分配到您的 2 个单独的 foo 属性中，然后与 ODE 集成相比，foo_vectors_3 的计算将是微不足道的。您还可以添加辅助函数以从一维 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.]]

【讨论】：