Scipy solve_ivp sensitivity to random phase shifts

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


• 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. Nov 12, 2022 at 1:13
• 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 ? Nov 12, 2022 at 17:51
• Also how do you get a random $\phi_j$? How many have you compared, and how have you sampled them? Nov 12, 2022 at 18:09

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.