运行时间沃宁：在Double_scalars和Runtime Warnning中遇到的无效值：在Double_scalars中遇到的溢出

Question

我正在编码NEWMARK BETA数值动态方程求解器，但我 有一个错误：

d:\university\payan nameh\cssr python codes\newmarkbeta_method.py:46: RuntimeWarning: overflow encountered in double_scalars
  p[I]=-1.556*u2[I]+a1*y[I-1]+a2*dy[I-1]+a3*dy2[I-1]
d:\university\payan nameh\cssr python codes\newmarkbeta_method.py:48: RuntimeWarning: invalid value encountered in double_scalars
  dy[I]=2*(y[I]-y[I-1])/dt-dy[I-1]
d:\university\payan nameh\cssr python codes\newmarkbeta_method.py:49: RuntimeWarning: invalid value encountered in double_scalars
  dy2[I]=4*(y[I]-y[I-1])/(dt**2)-4*dy[I-1]/dt-dy2[I-1]

我尝试使用小数点，但没有用。

from mpmath import *
import numpy as np
from numpy.core.numeric import ones_like, zeros_like
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint

chichi = pd.read_excel(r'D:\university\payan nameh\excel datas\chi-chi_w_scaled4g.xlsx')
u2 = chichi.iloc[:,1:].values
u2 = u2.transpose()
u2 = u2*9.81
u2 = u2.flatten()
TS = 1.5 #second
omega = (2*np.pi)/(TS)
kesai = 0.15
dt = 0.005
Ti = 89.995
L = 10
ky = 0.05
H = 3
Ts = np.linspace(0, Ti, int(Ti/dt+1.0))
Dt = zeros_like(Ts)
kmod = ky*9.81+11*Dt/L
d1 = 1-0.36*H
Y0init = zeros_like(Ts)
Y1init = zeros_like(Ts)
dy2 = zeros_like(Ts)
dy2_p = zeros_like(Ts)
dd2 = zeros_like(Ts)    
dd1 = zeros_like(dd2)
p = zeros_like(dd2)
p_prime = zeros_like(dd2)
y =  zeros_like(Ts)
y_p =  zeros_like(Ts)
dy =  zeros_like(Ts)
dy_p =  zeros_like(Ts)
a1 = 4/(dt**2)+(4*kesai*omega)/dt
a2 = 4/dt+2*kesai*omega
a3 = 1
k_hat = omega**2+a1
a1_prime = 4/(dt**2)+(4*kesai*omega)/(dt*d1)
a2_prime = 4/dt+2*kesai*omega/d1
a3_prime = 1
k_hat_prime = (omega**2)/d1+a1_prime
for I in range(1,len(Ts)):
   p[I]=-1.556*u2[I]+a1*y[I-1]+a2*dy[I-1]+a3*dy2[I-1]
   y[I]=p[I]/k_hat
   dy[I]=2*(y[I]-y[I-1])/dt-dy[I-1]
   dy2[I]=4*(y[I]-y[I-1])/(dt**2)-4*dy[I-1]/dt-dy2[I-1]
   if -kmod[I-1]-0.2313*H*dy2[I]-u2[I]>0:
       p_prime[I]=1.556*kmod[I-1]/d1+a1_prime*y[I-1]+a2_prime*dy[I-1]+a3_prime*dy2[I-1]
       y[I]=p_prime[I]/k_hat_prime
       dy[I]=2*(y[I]-y[I-1])/dt-dy[I-1]
       dy2[I]=4*(y[I]-y[I-1])/(dt**2)-4*dy[I-1]/dt-dy2[I-1]
       dd2[I]=-kmod[I-1]-0.2313*H*dy2[I]-u2[I]
       dd1[I]=dd2[I]*dt+dd1[I-1]
       if dd1[I]>0:
           Dt[I]=dd1[I]*dt+Dt[I-1]
           kmod[I]=ky*9.81+11*Dt[I]/L
       else:
           Dt[I]=Dt[I-1]
           kmod[I]=ky*9.81+11*Dt[I]/L
   else:
    Dt[I]=Dt[I-1]
    kmod[I]=ky*9.81+11*Dt[I]/L
plt.plot(Ts,y)
plt.show()