I have an earthquake record and I’m trying to obtain the spectral acceleration plot by solving the following differential equation:

T = np.geomspace(0.1,10,1000) # periods of vibration
SA_undamped = [] # empty list of spectral accelerations
times = earthquake_record.index

for Ti in T: 
    ksi = 0
    wn = 2 * np.pi / Ti
    sol = scipy.integrate.solve_ivp(sdof, t_span=(0,17.5), y0=[0,0], t_eval=times)
    SAi = np.max(np.abs(sol['y'][0])) * wn**2
    SA_undamped.append(SAi)

where the function sdof gives the derivatives of the two ODEs resulting from the initial differential equation.

def sdof(t, y):
    u = y[0]
    v = y[1]
    return np.vstack((v, agfun(t) - wn**2 * u - 2 * ksi * wn * v)).reshape(2,)

And agfun(t) gives the ground acceleration at time t.
This works fine. Then, I employed the Newmark-beta method for numerical integration of the differential equation:

def newmark_SD(earthquake_record, T, damping=0, x0=0, v0=0, a0=0):
    times = earthquake_record.index
    dt = times[1] - times[0]
    shp = earthquake_record.shape
    agvals = earthquake_record.values.reshape(shp[0],)
    
    ksi = damping
    wn = 2 * np.pi / T
    beta = 1/6
    gama = 1/2    
    SD = []
    SV = [] 
    SA = []
    
    xo = x0
    vo = v0
    ao = a0
    
    ag0 = agvals[0]
    ago = ag0
    
    for n,i in enumerate(times):
        if n > 0:
            agn = agvals[n]
            dpi_ = ((agn - ago) + (1 / (beta * dt) + gama / beta * 2 * ksi * wn) * vo +
                    (1 / (2*beta) + dt * (gama / (2*beta) - 1) * 2 * ksi * wn) * ao)
            ki_ = wn**2 + gama / (beta * dt) * 2 * ksi * wn + 1 / (beta * dt**2)
            
            xn = dpi_ / ki_ + xo
            
            vn = gama / (beta * dt) * (xn - xo) - gama / beta * vo + dt * (1 - gama / (2*beta)) * ao + vo
            
            an = 1 / (beta * dt**2) * (xn - xo) - 1 / (beta * dt) * vo - 1 / (2*beta) * ao + ao
            
            SD.append(xn)
            SV.append(vn)
            SA.append(an)
            
            xo = xn
            vo = vn
            ao = an
            ago = agn
    
    return (SD,SV,SA) 

And used the following to get the whole response:

SA_newmark = []
for Ti in T: 
    spectra = newmark_SD(earthquake_record, Ti)
    wn = 2 * np.pi / Ti
    SAi = np.max(np.abs(spectra[0])) * wn**2
    SA_newmark.append(SAi)

The diagrams for the two solutions are shown in the following picture, which is set to display the interval from 0.1 to 1 to be able to see the differences, as in the larger periods they disappear and the two diagrams match perfectly.

So, there is a shift in the lower periods, which disappears in the larger periods. I’m not sure what to think of it. For any specific period, the integration is performed across the whole record, so I’m even more confused as to why there is a shift only in the lower periods.