python - solve_ivp - Usar tamaños de paso adaptables con scipy.integrate.ode

scipy.integrate.rk45 example

Question

La documentación (breve) para scipy.integrate.ode dice que dos métodos ( dopri5 y dop853 ) tienen control de paso de página y salida densa. Al mirar los ejemplos y el código en sí, solo puedo ver una forma muy simple de obtener resultados de un integrador. A saber, parece que acaba de avanzar el integrador hacia adelante mediante un dt fijo, obtiene los valores de la función en ese momento y repite.

Mi problema tiene escalas de tiempo bastante variables, por lo que me gustaría obtener los valores en cualquier momento que necesite evaluar para alcanzar las tolerancias requeridas. Es decir, desde el principio, las cosas cambian lentamente, por lo que los pasos del tiempo de salida pueden ser grandes. Pero a medida que las cosas se ponen interesantes, los pasos de tiempo de salida tienen que ser más pequeños. En realidad, no quiero una salida densa a intervalos iguales, solo quiero los pasos de tiempo que usa la función de adaptación.

EDITAR: salida densa

Una noción relacionada (casi lo opuesto) es "salida densa", donde los pasos tomados son tan grandes como los pasos a seguir, pero los valores de la función se interpolan (usualmente con precisión comparable a la precisión del paso a paso) a lo que sea usted quiere. El scipy.integrate.ode subyacente scipy.integrate.ode es aparentemente capaz de esto, pero ode no tiene la interfaz. odeint , por otro lado, se basa en un código diferente, y evidentemente hace una producción densa. (Puede emitir cada vez que se llame a su lado derecho para ver cuándo sucede eso, y ver que no tiene nada que ver con los tiempos de salida).

Así que aún podría aprovechar la adaptabilidad, siempre que pudiera decidir los pasos de tiempo de salida que quiero con anticipación. Desafortunadamente, para mi sistema favorito, ni siquiera sé cuáles son las escalas de tiempo aproximadas como funciones de tiempo, hasta que ejecuto la integración. Así que tendré que combinar la idea de tomar un paso integrador con esta noción de producción densa.

Answer 1

El método de integrate acepta un step argumento booleano que le dice al método que devuelva un solo paso interno. Sin embargo, parece que los solucionadores ''dopri5'' y ''dop853'' no lo admiten.

El siguiente código muestra cómo puede obtener los pasos internos que realiza el solucionador cuando se usa el solucionador ''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()

Resultado:

Answer 2

He estado mirando esto para tratar de obtener el mismo resultado. Resulta que puede usar un truco para obtener los resultados paso a paso al establecer nsteps = 1 en la instanciación de oda. Generará un UserWarning en cada paso (esto puede capturarse y suprimirse).

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()

resultado:

Answer 3

Aquí hay otra opción que también debería funcionar con dopri5 y dop853 . Básicamente, el solucionador llamará a la función logistic() tantas veces como sea necesario para calcular los valores intermedios, de modo que almacenaremos los resultados:

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()

Answer 4

Para su información, aunque ya se ha aceptado una respuesta, debo señalar el registro histórico de que la salida densa y el muestreo arbitrario desde cualquier lugar a lo largo de la trayectoria calculada se admite de forma nativa en PyDSTool . Esto también incluye un registro de todos los pasos de tiempo determinados por adaptación utilizados internamente por el solucionador. Esto interactúa tanto con dopri853 como con radau5 y genera automáticamente el código C necesario para interactuar con ellos en lugar de basarse en devoluciones de llamadas de función python (mucho más lentas) para la definición del lado derecho. Ninguna de estas características se proporciona de forma nativa o eficiente en ningún otro solucionador enfocado en python, que yo sepa.

Answer 5

Desde SciPy 0.13.0 ,

Ahora se puede acceder a los resultados intermedios de la familia dopri de solucionadores de ODE mediante una función de devolución de llamada de solout .

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()

Resultado:

El resultado parece ser ligeramente diferente de Tim D''s, aunque ambos usan el mismo back-end. Sospecho que esto tiene que ver con la propiedad dopri5 de dopri5 . En el enfoque de Tim, creo que el resultado k7 de la séptima etapa se descarta, por lo que k1 se calcula de nuevo.

Nota: Hay un error conocido con set_solout que no funciona si lo configura después de establecer los valores iniciales . Se corrigió a partir de SciPy 0.17.0 .