Newmark方法和中心差分法求解单自由度时程曲线

 首先通过地震波绘制反应谱，随后根据周期为1s时的加速度反应谱值对地震加速度进行调幅，随后进行计算得到结构的位移时程曲线，速度时程曲线以及加速度反应谱

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import math
# ---------------------中心差分法----------------------- #
def Central_difference_method_displacenment( acceleration_g: list,time:list,step):
    length_t = len(time)
    displacement = [0.000] * length_t # 定义一个位移向量
    speed = [0] * length_t # 定义一个速度向量
    acceleration = [0] * length_t # 定义一个加速度向量
    for i in range(2, length_t): # 求解位移，单位为mm
        displacement[i] = (-acceleration_g[i] + (2 / step ** 2 - 4 * (math.pi) ** 2) * displacement[i - 1] + (
                    0.2 * math.pi / 2 / step - 1 / step ** 2) * displacement[i - 2]) / (
                                      1 / step ** 2 + 0.2 * math.pi / 2 / step)
    for i in range(length_t):
        displacement[i]=displacement[i] *1000
    displacement1 = displacement[:]
    for i in range(length_t):
        displacement1[i]=displacement1[i] /1000
    for i in range(1,len(displacement)-1): # 求解速度，单位为m/s^2
        speed[i]=(displacement1[i+1]-displacement1[i-1])/step/2
    for i in range(1, len(displacement) - 1): # 求解加速度，单位为g
        acceleration[i] = (displacement1[i + 1] -2*displacement1[i]+ displacement1[i - 1]) / step**2
        acceleration[i]=acceleration[i]/9.8
    return displacement,speed,acceleration
# ---------------------Newmark方法----------------------- #
def Newmark(acceleration_g:list,time:list,step):  # Newmark方法，β取1/6
    Newmark_length=len(time)*10
    Newmark_displacement = [0]*Newmark_length # 定义一个位移向量
    Nemark_speed = [0] * Newmark_length # 定义一个速度向量
    Nemark_acceleration = [0]*Newmark_length # 求解位移，单位为mm
    p1_m=[0]*(Newmark_length+1)
    acceleration_g_linear=[0]*Newmark_length
    for i in range(len(acceleration_g)-1):
        deta=acceleration_g[i+1]-acceleration_g[i]
        for j in range(i*10,i*10+10):
            if j ==i:
                acceleration_g_linear[j]=acceleration_g[i]
            else:
                acceleration_g_linear[j] = acceleration_g_linear[j-1]+0.1*deta
    a1_m = 6 / (step **2) + 3 / step * 0.2 * math.pi #系数a1与质量m的比值
    a2_m = 6 / step + 2 * 0.2 * math.pi #系数a2与质量m的比值
    a3_m = 2 + step / 2 * 0.2 * math.pi #系数a3与质量m的比值
    k1_m = 4 * math.pi ** 2 + a1_m
    for i in range(1,Newmark_length-1):
        p1_m[i]=-acceleration_g_linear[i]+a1_m*Newmark_displacement[i-1]+a2_m*Nemark_speed[i-1]+a3_m*Nemark_acceleration[i-1]
        Newmark_displacement[i]=p1_m[i]/k1_m # 求解位移
        Nemark_speed[i]=3/step*(Newmark_displacement[i]-Newmark_displacement[i-1])-2*Nemark_speed[i-1]-step/2\
                        *Nemark_acceleration[i-1] # 求解速度
        Nemark_acceleration[i]=6/step**2*(Newmark_displacement[i]-Newmark_displacement[i-1])-6/step*Nemark_speed[i-1]\