中心差分法解动力学方程
01 算法分析 将位移按照泰勒公式展开，得到前差分公式： 同样可得向后差分公式： 以上两式相减和相加分别得到： 以上两式忽略高阶小量，可得到 时刻速度和加速度表达式： 为了求解 时刻的位移，将 代入 时刻动力学方程 得到 其中 若已经求得 和 时刻的位移 和 ，则可以从 求得 时刻的位移。由 可知，只给定初值 和 不能求出 ，还必须确定 ，即该方法存在如何起步的问题。在向后差分公式 中取 得 其中 和 由初值条件给出。而 则由 求得。中心差分法解动力学方程的算法可归纳为(一)初始计算形成刚度矩阵 ，质量矩阵 和阻尼矩阵 由初值 和 求解 和 由时间步长 计算计算 计算有效质量矩阵 对 进行分解 (二) 对每一时间步计算 时刻的有效载荷解 时刻位移 如果需要，按照 计算 时刻的速度和加速度02 算例 用中心差分法解运动方程，时间步长 其中 初始条件 将初始条件代入方程，解得 每个时间步长计算 以及 03 编程实现 # 中心差分法# @表示矩阵乘法import numpy as npimport matplotlib.pyplot as pltstep = 10 #求解步数data = np.zeros( (2,step) ) # 存储各时间步长的数据K = np.mat('6,-2;-2,4')M = np.mat('2,0;0,1')Q = np.mat('0;10')a_0 = np.zeros((2,1))velocity_0 = np.zeros((2,1))# 求初始加速度acceleration_0 = np.linalg.inv(M) *(Q - K @ a_0)dt = 0.28c0 = 1/(dt**2)c1 = 0.5 * dtc2 = 2 * c0c3 = 1 / c2a_dt = a_0 - dt * velocity_0 + c3 * acceleration_0# 有效质量矩阵MM = c0 * MinvMM = np.linalg.inv(MM)TMP1 = K - c2*M A_t_sub_dt = a_dtA_t = a_0 A_t_plus_dt = np.zeros((2,1))for i in range(step): TMP2 = M @ A_t_sub_dt QQ = Q - TMP1 @ A_t - c0 * TMP2 A_t_plus_dt = invMM @ QQ data[0:,[i] ] = A_t_plus_dt A_t_sub_dt = A_t A_t = A_t_plus_dt#解线性方程组采用求逆的方法，计算规模大的，可以用LDLT分解.#time = np.zeros((step))for i in range(step): t = 0.28 * (i+1) time[i] = tlabels =['Δt','2Δt','3Δt','4Δt','5Δt','6Δt','7Δt','8Δt','9Δt','10Δt']fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(5,12) )ax1.plot(time, data[0, :], 'r-')ax1.set_xticks(time, labels)ax1.set_ylabel('$a_1$',fontsize = 14)ax2.plot(time, data[1, :], 'b-')ax2.set_xticks(time, labels)ax2.set_ylabel('$a_2$',fontsize = 14)fig.savefig('./f429.png', dpi = 400) plt.show()来源：数值分析与有限元编程
