编制矩阵位移法程序

矩阵位移法（Matrix Displacement Method）是一种常用的结构分析方法，可以用于求解各种复杂结构的内力、位移等问题。下面是一个简单的 MATLAB 程序示例，用于求解一个简单梁的内力和位移。

假设有一个简单支承的梁，长度为 L，截面积为 A，弹性模量为 E，载荷为 P。将梁分成 n 个小段，每个小段的长度为 h = L/n。在每个节点处，假设有两个未知位移 u 和 v，每段小梁的内力可以表示为以下矩阵形式：

$$
\begin{bmatrix}
N_i \\
Q_i
\end{bmatrix} =
\begin{bmatrix}
\frac{A_i}{h_i} & -\frac{A_i}{h_i} \\
-\frac{A_i}{h_i} & \frac{A_i}{h_i}
\end{bmatrix}
\begin{bmatrix}
u_{i-1} \\
v_{i-1}
\end{bmatrix} +
\begin{bmatrix}
0 \\
P_i
\end{bmatrix}
$$

其中，$N_i$ 和 $Q_i$ 分别表示第 i 个节点处的轴力和剪力，$A_i$ 表示该节点处的截面积，$P_i$ 表示该节点处的载荷。

根据上述公式，可以编写如下 MATLAB 程序：

% 输入参数
L = 5;       % 梁的长度
A = 1e-3;    % 梁的截面积
E = 2e11;    % 梁的弹性模量
P = 1e3;     % 梁的载荷
n = 10;      % 将梁分成 n 个小段

% 计算每段小梁的长度和节点编号
h = L/n;
x = linspace(0, L, n+1);

% 初始化矩阵和向量
K = zeros(2*n, 2*n);  % 刚度矩阵
f = zeros(2*n, 1);    % 载荷向量

% 计算每段小梁的刚度矩阵和载荷向量，并组装成总的刚度矩阵和载荷向量
for i = 1:n
    A_i = A;
    h_i = h;
    P_i = P/n;
    
    k_i = [A_i/h_i, -A_i/h_i; -A_i/h_i, A_i/h_i] * E/h_i;
    f_i = [0; P_i];
    
    K(2*i-1:2*i, 2*i-1:2*i) = K(2*i-1:2*i, 2*i-1:2*i) + k_i;
    f(2*i-1:2*i) = f(2*i-1:2*i) + f_i;
end

% 设置边界条件，第一个节点 u=0，最后一个节点 v=0
K(1, :) = [1, zeros(1, 2*n-1)];
K(2*n, :) = [zeros(1, 2*n-1), 1];
f(1) = 0;
f(2*n) = 0;

% 解方程，得到节点位移和内力
u_v = K\f;
u = u_v(1:2:end);
v = u_v(2:2:end);

% 计算每个节点处的轴力和剪力
N = A * E * diff(u) / h;
Q = P/n + A * E * diff(v) / h;

% 绘制位移和内力图
subplot(2, 1, 1);
plot(x, [0; u; 0]);  % 位移图
xlabel('Position (m)');
ylabel('Displacement (m)');
title('Displacement Diagram');
grid on;

subplot(2, 1, 2);
plot(x(1:end-1), N); hold on;  % 轴力图
plot(x(1:end-1), Q);          % 剪力图
xlabel('Position (m)');
ylabel('Internal Force (N)');
title('Internal Force Diagram');
legend('Axial Force', 'Shear Force');
grid on;

上述程序将梁分成了 10 个小段，计算了每个节点处的位移、轴力和剪力，并绘制了位移图和内力图。你可以根据具体的问题和要求进行修改和扩展。