I am trying to solve a coupled system of ODE's using the solve_ivp function from scipy. The general form of the equation is given via $$\dot{y}(t) = M(t)y(t).$$ The time dependence of matrix is basically sourced by a sum of weighted cosines, and each term (should) be shifted by some random phase factor. The input signal is given by the general form $$\phi(t) = \sum_{j}A_j \cos (\omega_jt+\phi_j),$$ where $\phi_j \in [0,2\pi)$. Now I am looking at the power spectrum of the signal, and upon inclusion of the random phases, the only change I see is a shift of the onset of the spectrum from $1$ onto $0$. Now my question is, could the source of this change be some numerical artifact in the sense of numerical ode solvers being sensitive to such phase shifts?


The specific form of the matrix is $$ M(t)= \begin{pmatrix}0&-\gamma&0\\ \gamma & 0 & \phi(t) \\ 0 & -\phi(t) & 0 \end{pmatrix},$$ and the initial condition is that $y(0) = \hat{e}_z$. The output I am looking at is the $x$-component of the vector returned by solve_ivp, in particular the power spectrum of this signal.

Edit 2:

I am sampling the random phases as phi = 2*np.pi*np.random.sample(N), such that it goes into the input signal as

def input(w,t):
    input = 0
    wlist = np.linspace(0,1,N)
    phi = 2*np.pi*np.random.rand(N)
    for i in range(N):
        input += np.sqrt(f(wlist[i])*1/N)*np.cos(2*np.pi*w*(1+wlist[i]**2/2)*t+phi[i])
    return input

# Full input as the average of U sampled inputs

def input_full(w,t,U):
    input_full = 0
    for i in range(U):
        input_full += input(w,t)
    return 1/U*input_full

The function $f(w)$ going into the definition is a PDF which gives the weight for every mode entering the sum of cosines defined above. Then this is given into solve_ivp via

#defining inital state, span of times, number of samples, and coupling
sample = 100
dt = 1/4
U = 50
wC = .01
MINI  = np.array([0,0,1])
t_span = np.array([0,sample*dt])
tlist = np.linspace(0,sample*dt,sample,endpoint=False)

#defining the rhs of the ODE
def rhs_full(t,M,x,w):
    return np.array([-(1-x)*w*M[1],(1-x)*w*M[0] - wC*input_full(w,t,U)*M[2],wC*input_full(w,t,U)*M[1]])

#the call to solve_ivp
sol = scipy.integrate.solve_ivp(rhs_full,t_span,MINI,args=p,t_eval=tlist,method="RK45")

  • 3
    $\begingroup$ You'll have to provide a more complete description of the problem, along with a more complete description of the output you get. For example, explain the exact structure of $M(t)$ and tell us how the "input signal" enters the problem. $\endgroup$ Commented Nov 12, 2022 at 1:13
  • 1
    $\begingroup$ Can you add more details (e.g graphs) on your power Spektrum and it's calculation. Also, can you write down the call to solve_ivp ? $\endgroup$
    – Laurent90
    Commented Nov 12, 2022 at 17:51
  • 1
    $\begingroup$ Also how do you get a random $\phi_j$? How many have you compared, and how have you sampled them? $\endgroup$
    – user9794
    Commented Nov 12, 2022 at 18:09

1 Answer 1


I first thought you were integrating an ODE of the type $y'=f(y,p)$ with $p$ some random phase vector drawn at the start of the simulation. Actually, looking at your code, it seems you are not integrating a traditional ODE, but a stochastic one, since you draw $p$ at each call of $f$. It may be useful to look into more advanced stochastic integrators that are constructed for this kind of problems.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.