Using a forloop to solve coupled differential equations in python

Question

I am trying to solve a set of differential equations, but I have been having difficulty making this work. My differential equations contain an "i" subscript that represents numbers from 1 to n. I tried implementing a forloop as follows, but I have been getting this index error (the error message is below). I have tried changing the initial conditions (y0) and other values, but nothing seems to work. In this code, I am using solve_ivp. The code is as follows:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import solve_ivp

def testmodel(t, y):
    X = y[0]
    Y = y[1]
    J = y[2]
    Q = y[3]
    a = 3
    S = 0.4
    K = 0.8
    L = 2.3
    n = 100 
    for i in range(1,n+1):
        dXdt[i] = K**a+(Q[i]**a) - S*X[i]
        dYdt[i] = (K*X[i])-(L*Y[i])
        dJdt[i] = S*Y[i]-(K*Q[i])
        dQdt[i] = K*X[i]/L+J[i]
        return dXdt, dYdt, dJdt, dQdt
t_span= np.array([0, 120])
times = np.linspace(t_span[0], t_span[1], 1000) 
y0 = 0,0,0,0
soln = solve_ivp(testmodel, t_span, y0, t_eval=times, 
vectorized=True)
t = soln.t
X = soln.y[0]
Y = soln.y[1]
J = soln.y[2]
Q = soln.y[3]

plt.plot(t, X,linewidth=2, color='red')
plt.show()

The error I get is

IndexError                Traceback (most recent call last)
<ipython-input-107-3a0cfa6e42ed> in testmodel(t, y)
     15     n = 100
     16     for i in range(1,n+1):
 --> 17     dXdt[i] = K**a+(Q[i]**a) - S*X[i]
   
 IndexError: index 1 is out of bounds for axis 0 with size 1

I have scattered the web for a solution to this, but I have been unable to apply any solution to this problem. I am not sure what I am doing wrong and what to actually change.

I have tried to remove the "vectorized=True" argument, but then I get an error that states I cannot index scalar variables. This is confusing because I do not think these values should be scalar. How do I resolve this problem, my ultimate goal is to plot these differential equations. Thank you in advance.

hello, dXdt,dYdt, etc are not defined... Q is not a list, just a single number, 0, you cant use indexing on a single number — Kevin Choi, Feb 25 '21 at 03:33
I also recommend you to use `class` not a single function when you do this kind of project = ] — Kevin Choi, Feb 25 '21 at 03:34
My answer is to what you actually code, but some half-sentences sound like you want to do something like in [this question](https://stackoverflow.com/questions/44783694/applying-a-set-of-ordinary-differential-equations-to-each-grid-cell). Meaning you have a structure of connected cells or compartments with a time dynamic inside each cell and some migration, transport or other coupling along the connections between cells. A non-biological example is sets of pendulums or metronoms on buffered platforms where the platforms are weakly coupled, giving surprising resonance patterns. — Lutz Lehmann, Feb 25 '21 at 07:33
@KevinChoi Thanks for your answer! What would be the best way to define those equations as a list? And with that, how could I integrate the "class" function--I am not sure how I would do that without severely changing the structure of the equations. — , Feb 25 '21 at 18:51
@py2912 hello, if you are not familar with `class`, forget about it. plz try that after you get familiar with coding stuffs. Also, numpy array would be a great option for this kind of matrix or array oriented research. numpy array would be a better way than defining those eqs as lists. — Kevin Choi, Feb 26 '21 at 05:09

Lutz Lehmann · Answer 1 · 2021-02-26T08:29:07.977

It is nice that you provide the standard solver with a vectorized ODE function for multi-point evalutions. But the default method is the explicit RK45, and explicit methods do not use Jacobi matrices. So there is no need for multi-point evaluations for difference quotients for the partial derivatives.

In essence, the coordinate arrays always have size 1, as the evaluation is at a single point, so for instance Q is an array of length 1, the only valid index is 0. Remember, in all "true" programming languages, array indices start at 0. It is only some CAS script languages that use the "more mathematical" 1 as index start. (Setting n=100 and ignoring the length of the arrays provided by the solver is wrong as well.)

You can avoid all that and shorten your routine by taking into account that the standard arithmetic operations are applied element-wise for numpy arrays, so

def testmodel(t, y):
    X,Y,J,Q = y
    a = 3; S = 0.4; K = 0.8; L = 2.3
    dXdt = K**a + Q**a - S*X
    dYdt = K*X - L*Y
    dJdt = S*Y - K*Q
    dQdt = K*X/L + J
    return dXdt, dYdt, dJdt, dQdt

Modifying your code for multiple compartments with the same dynamic

You need to pass the solver a flat vector of the state. The first design decision is how the compartments and their components are arranged in the flat vector. One variant that is most compatible with the existing code is to cluster the same components together. Then in the ODE function the first operation is to separate out these clusters.

    X,Y,J,Q = y.reshape([4,-1])

This splits the input vector into 4 pieces of equal length. At the end you need to reverse this split so that the derivatives are again in a flat vector.

    return np.concatenate([dXdt, dYdt, dJdt, dQdt])

Everything else remains the same. Apart from the initial vector, which needs to have 4 segments of length N containing the data for the compartments. Here that could just be

    y0 = np.zeros(4*N)

If the initial data is from any other source, and given in records per compartment, you might have to transpose the resulting array before flattening it.

Note that this construction is not vectorized, so leave that option unset in its default False.

For uniform interaction patterns like in a circle I recommend the use of numpy.roll to continue to avoid the use of explicit loops. For an interaction pattern that looks like a network one can use connectivity matrices and masks like in Using python built-in functions for coupled ODEs

@Thanks so much for your answer. In my case, I am trying to simulate a biological model where i = 1,2,3...n and n represents the number of oscillators. The example I gave above is just a modified version of the model that produces the same error I am getting. If I use your approach, then is there any way I can take this idea into account, or would it have to be a "single-oscillator" model? Since in your example, you do not use the "i" subscript at all. — , Feb 25 '21 at 18:50

Using a forloop to solve coupled differential equations in python

1 Answers1

Modifying your code for multiple compartments with the same dynamic