Newmark-beta法求解运动微分方程低周期结果滞后于真实解

我有一条地震记录，尝试通过求解运动微分方程获取谱加速度曲线。先用scipy.integrate.solve_ivp得到了正常结果，随后自己实现了Newmark-beta法数值积分，发现低振动周期（0.1-1s）下结果滞后于solve_ivp的解，但高周期时两者完全吻合。所有周期的积分都覆盖完整记录，困惑为何仅低周期出现滞后。

1. 参考解法（solve_ivp）代码 #

import numpy as np
import scipy.integrate

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)

对应的sdof函数：

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

注：agfun(t)用于获取对应时刻的地面加速度

2. 自定义Newmark-beta法代码 #

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) 

调用代码：

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)

核心原因：时间步长与结构周期的匹配性问题 #

Newmark-beta法的精度高度依赖时间步长dt与结构自振周期T的比值：

• 对于低周期结构（T=0.1-1s），自振频率wn=2π/T很高，对应的振动周期远小于地震记录的时间步长dt（假设地震记录是常见的0.01s或0.02s采样率，T=0.1s时，dt/T=0.1或0.2，已接近临界值）。
• 你使用的是beta=1/6, gama=1/2的线性加速度法（显式Newmark的一种），这种方法的稳定条件是dt ≤ T/π（无阻尼情况）。当dt接近或超过这个阈值时，数值解会出现相位滞后、振幅衰减的问题，低周期结构更容易触发这个条件。
• 高周期结构（T=1-10s）的dt/T比值很小（比如T=10s，dt=0.01s时dt/T=0.001），远低于稳定阈值，因此数值解精度足够，和solve_ivp的结果吻合。

代码中的具体问题 #

你的Newmark-beta实现存在两个关键错误，进一步放大了低周期下的误差：

1. 初始加速度计算错误：
   初始时刻（n=0）的结构加速度a0应该由运动微分方程推导得出：a0 = -wn²*x0 - 2*ksi*wn*v0 - agfun(t0)，但你直接将a0设为0，这会引入初始误差，低周期结构对初始条件更敏感，误差被放大。
2. 地面加速度项的处理错误：
   Newmark-beta法中，地面加速度的激励应该是绝对加速度的修正项，你的dpi_计算中使用agn - ago是错误的，正确的激励项应该是-agn（运动方程为ü + 2ksiwnu̇ + wn²u = -ag(t)），而非地面加速度的差值。

修正后的Newmark-beta代码 #

def newmark_SD_corrected(earthquake_record, T, damping=0, x0=0, v0=0):
    times = earthquake_record.index
    dt = times[1] - times[0]
    agvals = earthquake_record.values.flatten()
    
    ksi = damping
    wn = 2 * np.pi / T
    beta = 1/6
    gama = 1/2    
    SD = [x0]  # 保留初始位移
    SV = [v0]  # 保留初始速度
    SA = []
    
    # 计算初始加速度
    ao = -wn**2 * x0 - 2 * ksi * wn * v0 - agvals[0]
    SA.append(ao)
    
    xo = x0
    vo = v0
    
    for n in range(1, len(times)):
        agn = agvals[n]
        
        # 修正后的等效荷载与刚度
        dpi = -agn + wn**2 * xo + 2 * ksi * wn * vo + (1/(beta*dt**2))*xo + (gama/(beta*dt))*vo + (1/(2*beta))*ao
        ki = wn**2 + (gama/(beta*dt))*2*ksi*wn + 1/(beta*dt**2)
        
        xn = dpi / ki
        
        # 计算速度和加速度
        vn = vo + (1 - gama)*dt*ao + gama*dt*((xn - xo)/(beta*dt**2) - vo/(beta*dt) - ao/(2*beta))
        an = (xn - xo)/(beta*dt**2) - vo/(beta*dt) - ao/(2*beta)
        
        SD.append(xn)
        SV.append(vn)
        SA.append(an)
        
        xo = xn
        vo = vn
        ao = an
    
    return (SD, SV, SA)

额外建议 #

• 对于低周期结构，建议采用隐式Newmark-beta法（比如beta=1/4, gama=1/2的平均加速度法），这种方法是无条件稳定的，不会受dt/T比值的限制，精度更可靠。
• 验证时可以先测试无阻尼单自由度系统在简谐激励下的响应，对比解析解，确认代码正确性后再用于地震记录分析。

内容的提问来源于stack exchange，提问作者n6r5