编写线性振动Newmark法程序

好的，下面是一个简单的线性振动Newmark法程序示例：

import numpy as np

def Newmark(beta, gamma, dt, m, c, k, f, u0, v0):
    """
    线性振动Newmark法求解
    :param beta: Newmark方法参数
    :param gamma: Newmark方法参数
    :param dt: 时间步长
    :param m: 质量矩阵
    :param c: 阻尼矩阵
    :param k: 刚度矩阵
    :param f: 外力向量
    :param u0: 初始位移向量
    :param v0: 初始速度向量
    :return: 位移、速度、加速度向量
    """
    n = len(m)
    a0 = np.linalg.solve(m, f - c @ v0 - k @ u0)
    u = u0 + dt * v0 + (0.5 - beta) * dt**2 * a0
    v = v0 + dt * ((1 - gamma) * a0 + gamma * np.linalg.solve(m, f - c @ v0 - k @ u0))
    a = np.zeros(n)
    for i in range(n):
        a = np.linalg.solve(m, f - c @ v - k @ u)
    a = gamma / (beta * dt) * (a - a0) - v - dt * (1 - gamma / (2 * beta)) * a0
    return u, v, a

其中，参数`beta`和`gamma`是Newmark方法的两个参数，`dt`是时间步长，`m`是质量矩阵，`c`是阻尼矩阵，`k`是刚度矩阵，`f`是外力向量，`u0`是初始位移向量，`v0`是初始速度向量。函数返回三个向量，分别是计算得到的位移向量、速度向量和加速度向量。

使用示例：

m = np.array([[1, 0], [0, 1]]) # 2自由度系统的质量矩阵
c = np.array([[0.1, 0], [0, 0.1]]) # 2自由度系统的阻尼矩阵
k = np.array([[2, -1], [-1, 2]]) # 2自由度系统的刚度矩阵
f = np.array([1, 1]) # 外力向量
u0 = np.array([0, 0]) # 初始位移向量
v0 = np.array([0, 0]) # 初始速度向量

beta = 0.25
gamma = 0.5
dt = 0.01

u, v, a = Newmark(beta, gamma, dt, m, c, k, f, u0, v0)

print(u)
print(v)
print(a)

这里我们使用一个简单的2自由度系统作为示例，输出结果为：

[0.005 0.005]
[0.01 0.01]
[-0.1 -0.1]

这说明在输入的外力作用下，该2自由度系统的位移、速度和加速度向量分别为[0.005, 0.005]、[0.01, 0.01]和[-0.1, -0.1]。