直观感受动力学数值解法的误差

▲图1

图1所示的模型，弹簧刚度为    ，忽略其他因素，把小车从平衡位置往右移动    ，小车由初始位移开始自由振动，它经过平衡位置后，弹簧开始压缩，直到质量块达到    ，之后小车又开始往右运动。在该系统振动过程中没有摩擦力，也没有外力作用于小车。理论上，小车的位移区间是    

下面来看看数值解：

import numpy as np
import matplotlib.pyplot as plt

def main(steps, dt):
    k = 3    # 弹簧刚度
    x = 0.5   # 初始位移
    m = 1     #质量
    v1 = 0   # 初始速度
    x1 = x

    xxx = np.zeros(steps)

    for i in range(steps):
        a  = -k *x1 /m   # 加速度
        v2 = v1 + a* dt   # 下一时刻的速度
        x2 = x1 + v1*dt + 0.5 * a * dt **2# 下一时刻的位移
        # 上述三个公式是中学物理的公式
        v1 = v2
        x1 = x2
        xxx[i] = x2

    fig, 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方法求解，并和解析解对比。

import numpy as np
import matplotlib.pyplot as plt

# Newmark方法

def newmark1(steps, dt):

    M = 1
    K = 3
    F = 1

    U_0 = 0   #初始位移
    V_0 = 0   #初始速度
    A_0 = 0   #初始加速度

    # 参数
    alpha = 0.5; beta = 0.25
    c0 = 1/(beta * dt**2)
    c1 = alpha /(beta * dt)
    c2 = 1/(beta * dt)
    c3 = 1/(2*beta) - 1
    c4 = alpha /beta - 1
    c5 = dt * (alpha /(2*beta ) - 1)
    c6 = dt * (1- alpha )
    c7 = alpha * dt


    # 有效刚度矩阵
    KK = K + c0 * M 
    invKK = 1/KK


    data = np.zeros( steps ) # 存储各时间步长下的位移

    for i in range(steps):
        H = c0 * U_0 + c2 * V_0 + c3 * A_0
        F_1 = np.sin(0.5 * dt*i) * F + M * H   #简谐荷载
        U_1 = invKK * F_1
        data[i ] = U_1

        ## t+dt时刻的速度和加速度
        A_1 = c0 *(U_1 - U_0) - c2*V_0 -c3*A_0
        V_1 = V_0 + c6* A_0 + c7* A_1

        A_0 = A_1
        U_0 = U_1
        V_0 = V_1


    time = np.zeros(steps)
    data1 = np.zeros(steps)
    for i in range(steps):
        t = dt * i
        time[i] = t
        data1[i] = np.sin(0.5 * dt*i) * 1.2 * 0.4


    plt.rcParams['font.sans-serif'] = ['Noto Sans SC']    # 这个字体能够支持中英文字符
    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= 'upper right') 
    
    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)