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")
solve_ivp
? $\endgroup$