Scipy solve_ivp sensitivity to random phase shifts

Question

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?

Edit:

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")

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.