I am solving an ODE using scipy.solve_ivp where I explicitly define the time output. I save the data and time array in a file and open it in another program for the analysis. Because I only save part of the data to save disk space, I recreate the entire time array in the analysis code. Doing so, using identical expressions, I end up with time steps slightly different. Is this due to floating point errors, some other artifact or can it be the way I reproduce the time array? ODE solve program In the program where I solve the ODE I define the time array as follows np.linspace(0, n_freq_sweep*dt, n_freq_sweep) with omega_z = 422500*2*np.pi n_dt = 1000 # Resolution (number of time-steps for one omega_z period) dt = 2*math.pi/(n_dt*omega_z) # Time-step duration i_init = 2000 # number of periods n_init = n_dt*i_init # Number of time-steps in total time_simu = np.linspace(0,n_init*dt,n_init) sol = solve_ivp(derivs, t_span = [time_simu[0], time_simu[-1]], y0 = [2e-6,0], method='RK45', t_eval=time_simu, rtol = 1e-6) np.savez(filename, time_simu = time_simu[limits[0]:limits[1]], r_z = sol.y[0][limits[0]:limits[1]], v_z = sol.y[1][limits[0]:limits[1]]) where derivs is a function computing the derivatives of r_z and v_z at the next time step, given initial conditions (y0). I then save sol using np.savez, only between two limits I define somewhere. Note that I only show you a part of the program, the next step is to solve_ivp for 41 other time periods with different conditions, this would make GB of data so I so solve_ivp() + savez() in a for loop. So the time array is not complete from 0 to end of simu time. Analysis program I open the data, retrieve r_z, v_z, and time_simu. I need to "fill in the gaps" of the time array to be able to track the phase of r_z and v_z with respect to some frequency reference. I recreate the time array as omega_z = 422500*2*math.pi n_dt = 1000 dt = 2*math.pi/(n_dt*omega_z) i_init = 2000 n_init = n_dt*i_init time_secular = np.linspace(0,n_init*dt,n_init) # time_simu is the time imported from the data array Time step comparison Then if I compare the time step I get from the data array (time_simu) with the time I get from the analysis code (time_secular) they are different by 2.8820544184511335e-19. print(time_simu[2] - time_simu[1]) print(time_secular[2] - time_secular[1]) outputs 2.366865088469783e-09 2.3668650887579883e-09 # ^^^^^^^ Discussion So the difference is 2.8820544184511335e-19, which is way too much to be a missing point in the linspace, a + or - 1 in i_init, in linspace, or whenever else. I did the operation in the same order, the same way, to avoid different floating point errors, but maybe I am missing something. solve_ivp cannot be the cause because it only outputs x and v, but time I save it directly from the array. Can it be np.savez? Continue reading...
