[PYTHON] scipy.integrate.ode에 적응 단계 크기 사용

scipy.integrate.ode에 적응 단계 크기 사용

scipy.integrate.ode에 대한 (간략한) 문서는 두 가지 방법 (dopri5와 dop853)이 stepize 제어와 dense 출력을 가지고 있다고 말합니다. 예제와 코드 자체를 살펴보면 통합 자로부터 출력을 얻을 수있는 매우 간단한 방법 만 볼 수 있습니다. 즉, 여러분은 적분기를 일부 고정 된 dt만큼 앞으로 이동시키고 그 시간에 함수 값을 얻은 다음 반복하는 것처럼 보입니다.

내 문제는 상당히 가변적 인 시간 계를 가지고 있으므로 필요한 공차를 달성하기 위해 필요한 시간 단계마다 값을 얻고 싶습니다. 즉, 초기에 상황이 천천히 변하기 때문에 출력 시간 단계가 커질 수 있습니다. 그러나 일이 흥미로워 짐에 따라, 출력 시간 간격은 더 작아야합니다. 저는 평등 한 출력에서 ​​밀도가 높은 출력을 실제로 원하지 않습니다. 단지 적응 형 함수가 사용하는 시간 간격을 원합니다.

관련 개념 (거의 정반대)은 "밀도가 높은 출력"이며, 취해진 단계는 스테퍼가 신경 쓰는만큼 크지 만 함수의 값은 (일반적으로 스테퍼의 정확도와 비교할 수있는 정확도로) 보간됩니다 네가 원해. Fortran의 기본 scipy.integrate.ode는 분명히이 기능을 사용할 수 있지만 ode에는 인터페이스가 없습니다. 반면에 odeint는 다른 코드를 기반으로하며 밀도가 높은 출력을 나타내는 것으로 보입니다. 오른쪽 사이드가 호출 될 때마다 출력 할 수 있으며, 출력 시간과 관련이 없음을 알 수 있습니다.

그래서 나는 미리 원하는 출력 시간 단계를 결정할 수 있다면 적응 능력을 이용할 수 있습니다. 불행히도, 내가 좋아하는 시스템의 경우, 통합을 실행할 때까지 대략적인 시간 척도가 시간의 함수인지 알지조차 모릅니다. 그래서 나는 하나의 통합 단계를 취하는 개념을 고밀도 출력 개념과 결합해야 할 것입니다.

분명 scipy 1.0.0은 새로운 인터페이스를 통해 고밀도 출력을 지원했습니다. 특히 scipy.integrate.odeint에서 scipy.integrate.solve_ivp로 이동하는 것이 좋습니다. scipy.integrate.solve_ivp는 키워드 dense_output으로 사용됩니다. True로 설정하면 리턴 된 오브젝트는 시간의 배열로 호출 할 수있는 속성 sol을 가지며, 그 시간에 통합 함수 값을 리턴합니다. 그래도이 질문에 대한 문제는 해결되지 않지만 많은 경우에 유용합니다.

해결법

1. ==============================

   1.SciPy 0.13.0부터,

   SciPy 0.13.0부터,

   import numpy as np
   from scipy.integrate import ode
   import matplotlib.pyplot as plt
   
   
   def logistic(t, y, r):
       return r * y * (1.0 - y)
   
   r = .01
   t0 = 0
   y0 = 1e-5
   t1 = 5000.0
   
   backend = 'dopri5'
   # backend = 'dop853'
   solver = ode(logistic).set_integrator(backend)
   
   sol = []
   def solout(t, y):
       sol.append([t, *y])
   solver.set_solout(solout)
   solver.set_initial_value(y0, t0).set_f_params(r)
   solver.integrate(t1)
   
   sol = np.array(sol)
   
   plt.plot(sol[:,0], sol[:,1], 'b.-')
   plt.show()
   

   결과:

   두 팀 모두 동일한 백엔드를 사용하지만 팀 D와 약간 다른 결과가 나타납니다. 나는 이것이 dopri5의 FSAL 속성과 관련이 있다고 의심한다. 팀의 접근 방식에서, 나는 7 번째 단계의 결과 k7이 파기 된 것으로 생각하여 k1이 새로 계산되었다.

   참고 : 초기 값을 설정 한 후에 set_solout을 설정하면 작동하지 않는 알려진 버그가 있습니다. SciPy 0.17.0부터 수정되었습니다.

2. ==============================

   2.나는 동일한 결과를 얻으려고 이것을보고 있었다. ode 인스턴스화에서 nsteps = 1을 설정하여 단계별 결과를 얻으려면 해킹을 사용할 수 있습니다. 모든 단계에서 UserWarning을 생성합니다 (이 경우에는 catch되고 억제 될 수 있음).

   나는 동일한 결과를 얻으려고 이것을보고 있었다. ode 인스턴스화에서 nsteps = 1을 설정하여 단계별 결과를 얻으려면 해킹을 사용할 수 있습니다. 모든 단계에서 UserWarning을 생성합니다 (이 경우에는 catch되고 억제 될 수 있음).

   import numpy as np
   from scipy.integrate import ode
   import matplotlib.pyplot as plt
   import warnings
   
   
   def logistic(t, y, r):
       return r * y * (1.0 - y)
   
   r = .01
   t0 = 0
   y0 = 1e-5
   t1 = 5000.0
   
   #backend = 'vode'
   backend = 'dopri5'
   #backend = 'dop853'
   
   solver = ode(logistic).set_integrator(backend, nsteps=1)
   solver.set_initial_value(y0, t0).set_f_params(r)
   # suppress Fortran-printed warning
   solver._integrator.iwork[2] = -1
   
   sol = []
   warnings.filterwarnings("ignore", category=UserWarning)
   while solver.t < t1:
       solver.integrate(t1, step=True)
       sol.append([solver.t, solver.y])
   warnings.resetwarnings()
   sol = np.array(sol)
   
   plt.plot(sol[:,0], sol[:,1], 'b.-')
   plt.show()
   

   결과:

3. ==============================

   3.integrate 메소드는 단일 내부 단계를 리턴하도록 메소드에 알리는 부울 인수 단계를 채택합니다. 그러나 'dopri5'및 'dop853'솔버는이를 지원하지 않습니다.

   integrate 메소드는 단일 내부 단계를 리턴하도록 메소드에 알리는 부울 인수 단계를 채택합니다. 그러나 'dopri5'및 'dop853'솔버는이를 지원하지 않습니다.

   다음 코드는 'vode'솔버가 사용될 때 솔버가 취한 내부 단계를 얻는 방법을 보여줍니다.

   import numpy as np
   from scipy.integrate import ode
   import matplotlib.pyplot as plt
   
   
   def logistic(t, y, r):
       return r * y * (1.0 - y)
   
   r = .01
   
   t0 = 0
   y0 = 1e-5
   t1 = 5000.0
   
   backend = 'vode'
   #backend = 'dopri5'
   #backend = 'dop853'
   solver = ode(logistic).set_integrator(backend)
   solver.set_initial_value(y0, t0).set_f_params(r)
   
   sol = []
   while solver.successful() and solver.t < t1:
       solver.integrate(t1, step=True)
       sol.append([solver.t, solver.y])
   
   sol = np.array(sol)
   
   plt.plot(sol[:,0], sol[:,1], 'b.-')
   plt.show()
   

   결과:

4. ==============================

   4.FYI는 이미 답변을 받아 들였지만, 계산 된 궤도를 따라 어디에서나 밀도가 높은 출력 및 임의 샘플링이 PyDSTool에서 기본적으로 지원된다는 역사적인 기록을 지적해야합니다. 여기에는 해석기가 내부적으로 사용하는 모든 적응 적으로 결정된 시간 단계에 대한 기록도 포함됩니다. 이것은 dopri853과 radau5와의 인터페이스이며, 오른쪽 정의를위한 (훨씬 느린) 파이썬 함수 콜백에 의존하기보다는 그것들과 인터페이스하는 데 필요한 C 코드를 자동 생성합니다. 이 모든 기능은 내 지식에 비추어 볼 때 파이썬 중심의 다른 솔버에서는 기본적으로 또는 효율적으로 제공되지 않습니다.

   FYI는 이미 답변을 받아 들였지만, 계산 된 궤도를 따라 어디에서나 밀도가 높은 출력 및 임의 샘플링이 PyDSTool에서 기본적으로 지원된다는 역사적인 기록을 지적해야합니다. 여기에는 해석기가 내부적으로 사용하는 모든 적응 적으로 결정된 시간 단계에 대한 기록도 포함됩니다. 이것은 dopri853과 radau5와의 인터페이스이며, 오른쪽 정의를위한 (훨씬 느린) 파이썬 함수 콜백에 의존하기보다는 그것들과 인터페이스하는 데 필요한 C 코드를 자동 생성합니다. 이 모든 기능은 내 지식에 비추어 볼 때 파이썬 중심의 다른 솔버에서는 기본적으로 또는 효율적으로 제공되지 않습니다.

5. ==============================

   5.dopri5와 dop853에서도 사용할 수있는 또 다른 옵션이 있습니다. 기본적으로 해석기는 중간 값을 계산하는 데 필요한만큼 logistic () 함수를 자주 호출하므로 결과가 저장됩니다.

   dopri5와 dop853에서도 사용할 수있는 또 다른 옵션이 있습니다. 기본적으로 해석기는 중간 값을 계산하는 데 필요한만큼 logistic () 함수를 자주 호출하므로 결과가 저장됩니다.

   import numpy as np
   from scipy.integrate import ode
   import matplotlib.pyplot as plt
   
   sol = []
   def logistic(t, y, r):
       sol.append([t, y])
       return r * y * (1.0 - y)
   
   r = .01
   
   t0 = 0
   y0 = 1e-5
   t1 = 5000.0
   # Maximum number of steps that the integrator is allowed 
   # to do along the whole interval [t0, t1].
   N = 10000
   
   #backend = 'vode'
   backend = 'dopri5'
   #backend = 'dop853'
   solver = ode(logistic).set_integrator(backend, nsteps=N)
   solver.set_initial_value(y0, t0).set_f_params(r)
   
   # Single call to solver.integrate()
   solver.integrate(t1)
   sol = np.array(sol)
   
   plt.plot(sol[:,0], sol[:,1], 'b.-')
   plt.show()
   

from https://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode by cc-by-sa and MIT license