直观感受动力学数值解法的误差
▲图1图1所示的模型，弹簧刚度为，忽略其他因素，把小车从平衡位置往右移动，小车由初始位移开始自由振动，它经过平衡位置后，弹簧开始压缩，直到质量块达到，之后小车又开始往右运动。在该系统振动过程中没有摩擦力，也没有外力作用于小车。理论上，小车的位移区间是下面来看看数值解：importnumpyasnpimportmatplotlib.pyplotaspltdefmain(steps,dt):k=3#弹簧刚度x=0.5#初始位移m=1#质量v1=0#初始速度x1=xxxx=np.zeros(steps)foriinrange(steps):a=-k*x1/m#加速度v2=v1+a*dt#下一时刻的速度x2=x1+v1*dt+0.5*a*dt**2#下一时刻的位移#上述三个公式是中学物理的公式v1=v2x1=x2xxx[i]=x2fig,ax=plt.subplots(1,1,figsize=(10,6))ax.plot(range(steps),xxx,'r-')ax.set_xlabel('$t/s$',fontsize=14)ax.set_ylabel('$x2/m$',fontsize=14)fig.savefig('f46.png',dpi=400)plt.show()if__name__=="__main__":main(200,0.1)▲图2由图2可知，小车的振幅在不断增加，已经超过了理论上的这个范围了。▲图3考虑图3所示的一个时间步长。存在两个时间点和。粗线表示弹簧的内力。在此期间，小车正在移动，弹簧内力也随之变化。加速度与力直接相关，因此随着力的变化，加速度也会变化。然而，在同一时间段内，Python脚本使用的是时刻的弹簧内力计算了加速度的值，并在整个时间步长内使用该值。在代码中，加速度和力在整个时间步长内保持不变。因此，数值模拟中的力不是图3所示的曲线，而是一条直线（如虚线所示）。▲图4实际上，物理方程描述的是一个在时间上连续的系统，而Python脚本描述的是一个具有离散时间步长的系统，图4展示了实际力（曲线）和用于数值计算的力（阶梯线）。在附近，由于曲线趋于将阶梯线一分为二，因此数值模拟产生的误差不大。误差是两条曲线之间的面积，且曲线上下两侧的面积大致相等。在顶部和底部峰值处，计算力并不能模拟实际的力。因此，精度损失会逐渐累积。▲图5来看另一个例子。如图5所示的单自由度体系，简谐荷载作用在质点上,。根据图5b的弯矩图可求得柔度系数自振频率为动力系数则质点的动位移表达式令，用Newmark方法求解，并和解析解对比。importnumpyasnpimportmatplotlib.pyplotasplt#Newmark方法defnewmark1(steps,dt):M=1K=3F=1U_0=0#初始位移V_0=0#初始速度A_0=0#初始加速度#参数alpha=0.5;beta=0.25c0=1/(beta*dt**2)c1=alpha/(beta*dt)c2=1/(beta*dt)c3=1/(2*beta)-1c4=alpha/beta-1c5=dt*(alpha/(2*beta)-1)c6=dt*(1-alpha)c7=alpha*dt#有效刚度矩阵KK=K+c0*MinvKK=1/KKdata=np.zeros(steps)#存储各时间步长下的位移foriinrange(steps):H=c0*U_0+c2*V_0+c3*A_0F_1=np.sin(0.5*dt*i)*F+M*H#简谐荷载U_1=invKK*F_1data[i]=U_1##t+dt时刻的速度和加速度A_1=c0*(U_1-U_0)-c2*V_0-c3*A_0V_1=V_0+c6*A_0+c7*A_1A_0=A_1U_0=U_1V_0=V_1time=np.zeros(steps)data1=np.zeros(steps)foriinrange(steps):t=dt*itime[i]=tdata1[i]=np.sin(0.5*dt*i)*1.2*0.4plt.rcParams['font.sans-serif']=['NotoSansSC']#这个字体能够支持中英文字符plt.rcParams['axes.unicode_minus']=False#正常显示负号fig,ax=plt.subplots(1,1,figsize=(8,4))ax.plot(time,data,'r-',label='数值方法')ax.plot(time,data1,color='k',linestyle='--',label='解析法')ax.legend(frameon=True,#启用边框framealpha=0.8,#设置透明度fancybox=False,#边框圆角与否shadow=True,#是否启用阴影facecolor='lightgrey',#背景颜色edgecolor='dimgrey',#边框颜色loc='upperright')ax.set_xlabel('$t/s$',fontsize=14)ax.set_ylabel('$y/m$',fontsize=14)fig.savefig('f99.png',dpi=300)plt.show()if__name__=="__main__":newmark1(200,0.2)▲图6由图6可知，数值解和解析解在平衡位置相对误差不大，而在峰值位置相对误差较大。来源：数值分析与有限元编程
