用MATLAB得到一维变带宽存储整体刚度矩阵

假设有 $n$ 个节点，每个节点有 $d$ 个自由度，即总自由度数为 $nd$。设每个节点的位移为 $u_i$，则整体位移向量为 $\mathbf{u}=[u_1,u_2,\cdots,u_n]^T$。根据有限元方法，整体刚度矩阵可以表示为：

$$
\mathbf{K}=\sum_{e=1}^{ne} \mathbf{k}_e
$$

其中，$ne$ 为单元数，$\mathbf{k}_e$ 为第 $e$ 个单元的局部刚度矩阵。对于一维线性单元，其局部刚度矩阵为：

$$
\mathbf{k}_e=\frac{EA}{L_e}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}
$$

其中，$E$ 为模量，$A$ 为截面积，$L_e$ 为单元长度。为了方便计算，我们可以先将局部刚度矩阵变形为 $2d\times 2d$ 的形式，再将其按节点自由度排列为 $nd\times nd$ 的整体刚度矩阵。具体实现过程如下：

% 假设有 5 个节点，每个节点有 2 个自由度
n = 5;
d = 2;
nd = n * d;

% 定义单元长度、截面积和模量
L = 1;
A = 1;
E = 1;

% 定义节点坐标，这里假设节点等间距分布在 [0, L] 区间内
x = linspace(0, L, n);

% 定义单元节点编号
ele = [1, 2;
       2, 3;
       3, 4;
       4, 5];

% 定义局部刚度矩阵
ke = E * A / L * [1, -1; -1, 1];

% 初始化整体刚度矩阵
K = zeros(nd);

% 循环遍历每个单元，计算对应的局部刚度矩阵
for e = 1:size(ele, 1)
    % 获取该单元的两个节点编号
    n1 = ele(e, 1);
    n2 = ele(e, 2);
    
    % 计算该单元的长度
    Le = x(n2) - x(n1);
    
    % 将局部刚度矩阵变形为 2d x 2d 的形式
    ke2 = zeros(2*d);
    ke2(1:2, 1:2) = ke;
    ke2(3:4, 3:4) = ke;
    
    % 将局部刚度矩阵按节点自由度排列为 nd x nd 的整体刚度矩阵
    idx = [(n1-1)*d+1:n1*d, (n2-1)*d+1:n2*d];
    K(idx, idx) = K(idx, idx) + ke2 * E * A / Le;
end

最终得到的 $K$ 即为一维变带宽存储的整体刚度矩阵。