如何用MATLAB实现Newmark-β方法来计算单自由度体系的动力学响应？

Newmark-β方法是一种常用的数值积分方法，用于求解结构动力学方程。在MATLAB中，可以利用其内置的优化工具箱或者自定义函数来实现。以下是简化的步骤：

1. **设置系统模型**：
首先，你需要确定系统的刚度矩阵 `K`、质量矩阵 `M` 和外部荷载向量 `F`。对于单自由度体系，它们通常是简单的标量。

2. **定义时间步长**：
设定初始条件（如位移 `u0` 和速度 `v0`）、总时间跨度 `T`，以及时间步长 `dt`。

3. **选择Newmark参数**：
Newmark方法通常涉及两个参数：α和β。例如，你可以选择 α=0.5, β=0.25 的组合，这是经典Newmark法的标准值。

4. **创建积分函数**：
使用Newmark积分公式来估计加速度：

   a_n = (1 - β) * k * u_n + β * (k * u_n - k_n) / dt^2 + M * v_n;

其中，`u_n` 和 `v_n` 分别表示前一时刻的位移和速度，`k_n` 是前一时刻的内力。

5. **更新状态**：
根据加速度和当前时间步长计算新的位移和速度：

   u_(n+1) = u_n + v_n * dt + 0.5 * a_n * dt^2;
   v_(n+1) = v_n + 0.5 * (a_n + (1 - β) * a_(n+1)) * dt;

6. **循环迭代**：
对于每个时间步长，重复步骤4和5，直到达到总的模拟时间 `T`。

7. **结果存储**：
记录每一时刻的位移和速度，以便后续分析。

function [response] = newmark_beta(K, M, F, u0, v0, T, dt, alpha, beta)
    % 初始化
    response = [];
    u = u0; v = v0;
    
    % 时间迭代
    for t = 0:dt:T
        a = (1 - beta) * K * u + beta * (K * u - prev_k) / dt^2 + M * v; % 加速度
        u = u + v * dt + 0.5 * a * dt^2; % 更新位移
        v = v + 0.5 * (a + (1 - beta) * a) * dt; % 更新速度
        prev_k = K * u; % 存储内力
        
        % 结果存储
        response = [response; u; v];
    end
end