中心差分法解动力学方程

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

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()