2

I am setting up a rather complicated function that i wish to solve as an ODE. Preferably i want to use a solver which is accessible in most python environments so i am trying to use scipy.integrate.ode. The Problem is that my function produces a lot of values which i want to extract at each iteration.

Here is an example: I want to get someValues at each completed iteration of r.integrate() and save them somewhere:

from scipy.integrate import ode

def f(t, y):
    someValues = [10,5,4]
    return -1

y0 = 1

r = ode(f)
r.set_initial_value(y0, 0)

dt = 0.1
while r.successful() and r.y[0] >= 0:
    r.integrate(r.t + dt)
    # print(someValues)

However the solver forbids f() from returning additional values. Since it is not necessarily clear how the solver will acces the function, it is also not a good idea to just push it out to a global variable from inside f().

While it is possible to recalculate each step of f(), returning all values, after the solver is done, it would be preferrable to do it in one run for performance.

mathism
  • 63
  • 1
  • 6
  • I'm not sure if I understand the question completely, but are you basically trying to solve multiple (potentially coupled) first order odes? If so try odeint. Here is an example [here](https://stackoverflow.com/questions/66412602/is-this-correct-for-modeling-gravity-as-a-second-order-ode/66426700#66426700) – EntropicFroggy Mar 05 '21 at 06:02

1 Answers1

1

I have had this problem as well. With the Scipy ODE solvers, I do not think it is possible to do what you ask, and you must "recalculate each step of f()". Here is how I do it:

import numpy as np
from scipy.integrate import solve_ivp

def f_all(t, y):
    someValue = 8.0*y[0]
    dydt = y
    return dydt, someValue

def f(t, y):
    dydt, someValue = f_all(t, y)
    return dydt

def main():
    y0 = np.array([1.0])
    tspan = [0,1]
    t_eval = np.linspace(tspan[0],tspan[1], 100)
    sol = solve_ivp(f,tspan,y0,t_eval=t_eval,method='LSODA')
    ysol = sol.y.T

    someValue = np.empty((len(sol.t),),np.float64)
    for i in range(len(sol.t)):
        dydt, someValue[i] = f_all(sol.t[i], ysol[i])  
    return ysol, someValue

ysol, someValue = main()

Almost all of the time the ODE solve will take more time than doing a few hundred or thousand function calls after integration. But if you are concerned about performance, you can use the NumbaLSODA ode solver with numba to make all your code super fast:

import numpy as np
from NumbaLSODA import lsoda_sig, lsoda
import numba as nb

@nb.njit
def f_all(t, y):
    someValue = 8.0*y[0]
    dydt = y
    return dydt, someValue

@nb.cfunc(lsoda_sig)
def f(t, y_, dydt, p):
    y = nb.carray(y_, (1,))
    dydt_, someValue = f_all(t, y)
    dydt[0] = dydt_[0]
funcptr = f.address

@nb.njit
def main():
    y0 = np.array([1.0])
    tspan = [0,1]
    t_eval = np.linspace(tspan[0],tspan[1], 100)
    ysol, success = lsoda(funcptr, y0, t_eval)
    someValue = np.empty((len(t_eval),),np.float64)
    for i in range(len(t_eval)):
        dydt, someValue[i] = f_all(t_eval[i], ysol[i])  
    return ysol, someValue

ysol, someValue = main()

The NumbaLSODA version above is 100x faster than the Scipy version.

nicholaswogan
  • 631
  • 6
  • 13