0

I am working with a large system of ODEs that I automatically construct. I have a state array of length n, and an equal number of differential equations also stored in an array. Here there are 2, but in the real system, there would be over 100. The differential equations are constructed as a class and called via __call__, alternatively, ODEs are created and called with functools partial. I need to send the exact same state and time t to each differential function in the array and fill the derivative state vector dstate. I am wondering what the cleanest, most pythonic and, most importantly, fastest way of doing that is.

import numpy as np
from scipy import integrate
import pandas as pd

class F1:
    def __init__(self):
        pass
    def __call__(self, state):
        return 2 - state[0]**2 * state[1]

class F2:
    def __init__(self):
        pass
    def __call__(self, state):
        return 3 - state[1]*state[0]

fs = np.array([F1(), F2()])
state0 = np.ones(2)*4

def change(t, state):
    """
    This is the function I would like to change for something better!

    Ideally there would be something like: dstate = map(input_to_apply, list_of_functions)
    """
    dstate = np.zeros(2)
    for i,f in enumerate(fs):
        dstate[i] = f(state = state)
    return dstate

sol = integrate.solve_ivp(fun = change, t_span = (0, 15), y0 = state0)
res = pd.DataFrame(sol.y.T, columns = ['A', 'B'])
res.index = sol.t
ax = res.plot()
Patrickens
  • 321
  • 3
  • 14

1 Answers1

0

As for coding, there isn't (IMO) significant improvements on what you have done. You can replace single line loops:

def change(t, state):
    return np.array([f(state=state) for f in fs])

As for efficiency, check Calling functions by array index in Python. The best thing you can do is bundle the functions in one. In your example:

class F_list:
    def __init__(self):
        pass
    def __call__(self, state):
        return 2 - state[0]**2 * state[1], 3 - state[1]*state[0] 
    def change(self, t, state):
        return self(state)

F = F_list()
sol = integrate.solve_ivp(fun = F.change, t_span = (0, 15), y0 = state0)

Also note that using def or lambda is always faster than __call__, and if you want a class and unless you have a specific need for instances, you can use the classmethod decorator to avoid instantiating needlessly

Tarifazo
  • 4,118
  • 1
  • 9
  • 22