中心差分法解动力学方程

同样可得向后差分公式：

以上两式相减和相加分别得到：

以上两式忽略高阶小量，可得到    时刻速度和加速度表达式：

为了求解    时刻的位移，将    代入    时刻动力学方程

得到

其中

若已经求得    和    时刻的位移    和    ，则可以从    求得    时刻的位移。由    可知，只给定初值    和    不能求出    ，还必须确定    ，即该方法存在如何起步的问题。

在向后差分公式    中取    得

其中    和    由初值条件给出。而    则由    求得。

中心差分法解动力学方程的算法可归纳为

(一)初始计算

1. 形成刚度矩阵        ，质量矩阵        和阻尼矩阵        
2. 由初值        和        求解        和        
3. 由时间步长        计算
4. 计算        
5. 计算有效质量矩阵        
6. 对        进行分解        

(二) 对每一时间步

1. 计算        时刻的有效载荷
2. 解        时刻位移        
3. 如果需要，按照        计算        时刻的速度和加速度

其中

初始条件

将初始条件代入方程，解得

每个时间步长计算

以及

# 中心差分法
# @表示矩阵乘法

import numpy as np
import matplotlib.pyplot as plt
step = 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.28
c0 = 1/(dt**2)
c1 = 0.5 * dt
c2 = 2 * c0
c3 = 1 / c2
a_dt = a_0 - dt * velocity_0 + c3 * acceleration_0

# 有效质量矩阵
MM = c0 * M
invMM = np.linalg.inv(MM)
TMP1 = K - c2*M 
A_t_sub_dt = a_dt
A_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] = t

labels =['Δ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()