CODEWITHSUNDEEP
pythonscipydifferential-equations
kilojoules
kilojoules

Reputation: 10083

Solve a system of odes with scipy - how to reference different indexes?

I have a system of ODEs that depend on a matrix of data. Each ODE should reference a different column of data in its evaluation.

import numpy as np
n_eqns = 20
coeffs = np.random.normal(0, 1, (n_eqns, 20))
def dfdt(_, f, idx):
    return (f ** 2) * coeffs[idx, :].sum() - 2 * f * coeffs.sum()

from scipy.integrate import ode
f0 = np.random.uniform(-1, 1, n_eqns)
t0 = 0
tf = 1
dt = 0.001

r = ode(dfdt)
r.set_initial_value(f0, t0).set_f_params(range(n_eqns))

while r.successful() and r.t < tf:
    print(r.t+dt, r.integrate(r.t+dt))

How can I specify that each ODE should use the idx value associated with its index in the system of ODEs? The first equation should be passed idx=0, the second idx=1, and so on.

Upvotes: 2

Views: 263

Answers (1)

Wrzlprmft
Wrzlprmft

Reputation: 4415

The function dfdt takes and returns the state and derivative, respectively as arrays (or other iterables). Thus, all you have to do is to loop over all indices and apply your operations accordingly. For example:

def dfdt(t,f):
    output = np.empty_like(f)
    for i,entry in enumerate(f)
        output[i] = f[i]**2 * coeffs[i,:].sum() - 2*f[i]*coeffs.sum()
    return output

You can also write this using NumPy’s component-wise operations (which is quicker):

def dfdt(t,f):
    return f**2 * coeffs.sum(axis=1) - 2*f*coeffs.sum()

Finally note that using f for your state may be somewhat confusing since this is how ode denotes the derivative (which you call dfdt).

Upvotes: 2

Related Questions