使用odeint在Python中求解Springmass系统

Question

我正在用Python使用odeint数值求解谐振子的二阶常微分方程。从我的图下面的代码来看，系统似乎根本没有像预期的那样阻尼。有人能帮我调试代码吗？下面是到目前为止我编写的代码

Mass Displacement and velocity plot

    import numpy as np
    from scipy.integrate import odeint
    import matplotlib.pyplot as plt
        
    
    ########################### Given Parameters ######################
    l = 125e-6 #m
    w = 20e-6 #m
    h = 5e-6 #m
    rho = 2300 #Kg/m3
    Q = 300
    A_0 = 20e-9 #nm
    gamma_sv = 30e-3 #J/m2
    H = 6e-20 #jouls
    R = 20e-9 #m
    E_t = 130e9 # Pa
    v_t = 0.3
    E_s = 70e9 #Pa
    v_s = 0.3
    
    ################## Calculated Parameters from Homework7 #############################
    I = (w*h**3)/12
    mu = rho*w*h
    C_1 = 1.8751
    K_c = (3*E_t*I)/(l**3)
    w_0 = (C_1**2/l**2)*np.sqrt(E_t*I/mu)
    w_d = w_0
    f_0 = w_0/(2*np.pi)
    m = K_c/w_0**2
    gamma = m*w_0/Q
    F_0 = (w_0**2*m*A_0)/Q
    
    # F_D = (F_0 * np.cos(t))/(A_0*m*w_0**2) 
    
    E_star = (E_s*E_t)/(((1-v_t**2)*E_s)+((1-v_s**2)*E_t)) # Estar
    a_0 = np.sqrt((H)/(24*np.pi*gamma_sv)) #molecular distance meters
    
    #NonDimensionalized Constants
    lambda = 1/(A_0*K_c)
    beta = 1
    alpha = 1/Q
    
    #Initialization
    tstart    = 0
    tstop     = 20
    increment = 0.01
    
    #intital Conditions
    x_init = np.array([0,1])
    
    t = np.arange(tstart, tstop+1, increment)
    
    #Function returning the dxdt
    
    def mydiff(x,t):
        #gamma = m*w_0/Q                                 #Damping Constant
        K_c = (3*E_t*I)/(l**3)                           #Stiffness of the mass
        m = K_c/w_0**2                                   #Mass
        F_D = (F_0 * np.cos(t))/(A_0*m*w_0**2)           #Driving force
            
        dx1dt = x[1]
        dx2dt = F_D - alpha*x[1] - beta*x[0]      #NonDimensional Equation being solved
                
        dxdt  = [dx1dt, dx2dt]
        return dxdt
    
    # Solving the ODE
    x = odeint(mydiff, x_init, t)
    
    x1 = x[:,0]
    x2 = x[:,1]
    
    #Plotting the Results
    # plt.plot(t,y)
    
    plt.plot(t,x1)
    plt.plot(t,x2)
    #plt.title("Damped and Driven Harmonic Oscillator with F_ts = 0")
    plt.xlabel("t/s")
    plt.ylabel("x(t)")
    plt.legend(["Displacement (x1)", "Velocity (x2)"], loc="upper right")
    plt.grid()
    plt.show()