Returning odd result when solving differential equations with solve_ivp

Question

I am trying to solve a series of differential equations as shown in the code below. But when I do simulate these differential equations, I am getting an inaccurate result. Here, I have 4 variables and I first reshaped my vector and then concatenated the values in hopes of having equal vector lengths for each of my variables. With this, I set my initial condition as the number of variables * N (N represents the number of cells I want to simulate). The for loop illustrated in the code is for simulating populations of cells as opposed to single cells and thus I set my for loop from values 1 to N. The trouble I am having, however, is that when I return my results and plot the variables (X, Y, Z, V), my X, Z, and V values are at 0, when in reality I should be seeing oscillations. Additionally, my F value, which represents the mean-field, is supposed to also show oscillations but has negative values when I use a print statement (I am not supposed to get any negative values for any of these equations). I have checked my code many times over and I am not sure why I am not getting any oscillations when all of the equations were listed properly. I think I am not passing the variables correctly into the solver or if I need to use a different way of having vectors of equal length. Does anyone have any idea why I am not getting oscillations for X, Z, V, but am getting oscillations for Y? And how I can correct this error? Any help is much appreciated (code and plot are below).

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

def testing(t,y):
    X,Y,Z,V = y.reshape([4,-1])
    N = 1000; v1 = 6.8355; K1 = 2.7266; n = 5.6645; v2 = 8.4297; K2 = 0.2910 
    k3 = 0.1177; v4 = 1.0841; K4 = 8.1343; k5 = 0.3352; v6 = 4.6645 
    K6 = 9.9849; k7 = 0.2282; v8 = 3.5216; K8 = 7.4519; vc = 6.7924
    Kc = 4.8283; K = 1; L = 0 
    dXdt = np.zeros(N)
    dYdt = np.zeros(N)
    dZdt = np.zeros(N)
    dVdt = np.zeros(N)
    for i in range(1,N+1):
        F = 1/N * sum([V[i]])
        dXdt[i] = (v1*(K1**n))/((K1**n)+(Z[i]**n))-((v2*X[i])/(K2+X[i])) + ((vc*K*F)/(Kc+(K*F))) + L
        dYdt[i] = (k3*X[i]) - ((v4*Y[i])/(K4+Y[i]))
        dZdt[i] = (k5*Y[i]) - ((v6*Z[i])/(K6+Z[i]))
        dVdt[i] = (k7*X[i]) - ((v8*V[i])/(K8+V[i]))
        return np.concatenate([dXdt,dYdt,dZdt,dVdt])
t_span = (0, 96)
t = np.linspace(t_span[0], t_span[1], 3000) 
N = 1000
y0 = np.zeros(4*N)
soln = solve_ivp(testing, t_span, y0, t_eval=t)
t = soln.t; X = soln.y[0]; Y = soln.y[1]; Z = soln.y[2]; V = soln.y[3]

resulting plot: CLICK TO SEE PLOT

user_3088138 · Answer

There are several problems in the ODE function testing

The indentation of the return statement is wrong, it is in the loop, so triggered in the first run through it, for i=0.
F = 1/N * sum([V[i]]) is visibly wrong, you want to put F = 1/N * sum(V) outside the loop, so it is computed once for all i.
N should be the actual length of the received arrays, so that the mean calculation actually computes the mean, N=len(V).
The solver tries to evaluate the fractional power of Z at some negative value, probably in some intermediary stage in the RK45 method. Make it safe via absolute value, to keep some sign use Z*abs(Z)**(n-1).

In your interpretation of the solution, the indices are also not correct.

X = soln.y[0]; Y = soln.y[1]; Z = soln.y[2]; V = soln.y[3] does not correspond to the state component layout. The components for the first cell are X = soln.y[0]; Y = soln.y[N]; Z = soln.y[2*N]; V = soln.y[3*N], or in short, X,Y,Z,V = soln.y[0::N].

A corrected version using numpy array operations is

def testing(t,y):
    X,Y,Z,V = y.reshape([4,-1])
    N = len(X); v1 = 6.8355; K1 = 2.7266; n = 5.6645; v2 = 8.4297; K2 = 0.2910 
    k3 = 0.1177; v4 = 1.0841; K4 = 8.1343; k5 = 0.3352; v6 = 4.6645 
    K6 = 9.9849; k7 = 0.2282; v8 = 3.5216; K8 = 7.4519; vc = 6.7924
    Kc = 4.8283; K = 1; L = 0 
    F = 1/N * sum(V)
    dXdt = (v1*(K1**n))/((K1**n)+(Z*abs(Z)**(n-1)))-((v2*X)/(K2+X)) + ((vc*K*F)/(Kc+(K*F))) + L
    dYdt = (k3*X) - ((v4*Y)/(K4+Y))
    dZdt = (k5*Y) - ((v6*Z)/(K6+Z))
    dVdt = (k7*X) - ((v8*V)/(K8+V))
    return np.concatenate([dXdt,dYdt,dZdt,dVdt])

and gives, even for N=10, the nice oscillations

Returning odd result when solving differential equations with solve_ivp

1 Answers