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?

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?

• python
• arrays
• numpy
• floating-point

• 1 Are both operations done on the same machine using same Python installation? – matszwecja Commented 9 hours ago
• I am unable to get the first number, I am always getting 2.3668650887579883e-09 so I'd check if the code for time_simu is the same as you'd posted in question. – matszwecja Commented 8 hours ago
• The order of assignments is not chronological, ` n_dt` gets used before it is assigned a value (in the local context, it could still have a value from far away). To define a subdivision with ´n` intervals, use linspace(0, n*dt, n+1). Use the arange function where it is more appropriate. – Lutz Lehmann Commented 8 hours ago
• 1 I suppose that you might just want to check that math.pi and np.pi are exactly the same? In principle, those libraries might choose to store different precision. – lastchance Commented 8 hours ago
• 1 @MartinBrown : RK45 nor solve_ivp do not change the time steps in question. The sample times are externally given in t_eval and do not change. – Lutz Lehmann Commented 7 hours ago

1 Answer 1

I guess that you are not comparing the difference at the exact same position. Notice how you store

np.savez(filename, time_simu = time_simu[limits[0]:limits[1]], ...)

That is: You are not storing complete time_simu, but only in the range limits[0]:limits[1] (What exactly are the values of limits[0] and limits[1]?). When you later reload your stored values, you only reload the stored range, so what is now at index 2 (namely at index 2 of the reloaded time_simu[limits[0]:limits[1]]) is not necessarily the same as what was at index 2 before storing (namely at index 2 of time_simu, without slicing it).

In other words, when you compare …

print(time_simu[2] - time_simu[1])
print(time_secular[2] - time_secular[1])

… you are not necessarily comparing the difference between the same two indices, because in time_secular, you do not seem to apply the slice limits[0]:limits[1] before comparison.

The problem here is that when you create time_secular = np.linspace(0,n_init*dt,n_init) and, likewise, time_simu = np.linspace(0,n_init*dt,n_init), the differences between different neighboring indices will not all be exactly the same. So when you look at the difference of successors in different places, you may find slightly different results. You can check this, for example, as follows:

import math
import numpy as np

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)
reference_diff = time_secular[2] - time_secular[1]

print("reference difference:", reference_diff)
all_diffs = set()
for t_up, t_lo in zip(time_secular[1:], time_secular[:-1]):
    all_diffs.add(t_up - t_lo)
print("\nall occurring diffs in `time_secular`:")
print(", ".join(str(d) for d in sorted(all_diffs)))

And you will (hopefully) get the following result:

reference difference: 2.3668650887579883e-09

all occurring diffs in `time_secular`:
2.366865088469783e-09, 2.3668650886866233e-09, 2.3668650887408334e-09,
2.366865088754386e-09, 2.366865088757774e-09, 2.3668650887579858e-09,
2.3668650887579874e-09, 2.3668650887579883e-09, 2.366865088757989e-09,
2.3668650887579924e-09, 2.366865088757999e-09, 2.3668650887580123e-09,
2.3668650887580387e-09, 2.3668650887580917e-09, 2.3668650887581975e-09,
2.366865088758621e-09, 2.366865088759468e-09, 2.366865088761162e-09,
2.3668650887679384e-09, 2.3668650887950435e-09, 2.3668650889034637e-09,
2.3668650893371446e-09

Note that the second value that you found, 2.366865088469783e-09, is also part of all the occurring differences in time_secular.

As you do not show when or how you "fill in the gaps" in the reloaded time array, I cannot further pin down the problem. It might be as simple as an off-by-one error in indexing, or really due to the fact that you compare the difference at index 1,2 of the sliced, reloaded time_simu with the difference at index 1,2 of the unsliced, reconstructed time_secular.