用matlab代码求解板的质量矩阵和刚度矩阵，从而求解1-10阶模态振型

在MATLAB中，对于板的质量矩阵和刚度矩阵的求解通常涉及到结构力学中的有限元素法(FEM)。以下是一个简单的例子，假设我们有一个二维矩形薄板，使用一阶四节点四边形单元（Quad4），并忽略了剪切效应。这只是一个基础的步骤，实际代码会更复杂，包括网格划分、边界条件设定等。

首先，我们需要创建一些函数来计算节点位移、面积、单元刚度矩阵（Ke）和质量矩阵（Me）：

function [Ke, Me] = element_stiffness_mass(nodeCoord, E, nu, h, A)
    % nodeCoord 是每个节点的坐标数组，E, nu 是材料属性，h 是厚度
    % A 是单个单元面积
    
    Ke = zeros(4);  % 初始化单元刚度矩阵
    Me = zeros(4);  % 初始化单元质量矩阵
    
    % 计算节点位置相对于单元中心的位置
    d1 = nodeCoord(2,:) - nodeCoord(1,:);
    d2 = nodeCoord(4,:) - nodeCoord(3,:);
    
    % 计算相关的长度和角度
    l1 = sqrt(d1.^2 + d2.^2);
    alpha = atan2(d2, d1);
    
    % 一阶四节点单元的刚度矩阵和质量矩阵计算
    Ke = (1/(12*E*h*(1-nu^2))) * [
        (1-nu) * (d1.*d1 + d2.*d2), -(1+nu) * d1.*d2, -(1+nu) * d1.*d2, (1-nu) * (d1.*d1 + d2.*d2);
        -(1+nu) * d1.*d2, (1-nu) * (d1.*d1 + d2.*d2), (1-nu) * (d1.*d1 + d2.*d2), -(1+nu) * d1.*d2;
        -(1+nu) * d1.*d2, (1-nu) * (d1.*d1 + d2.*d2), (1-nu) * (d1.*d1 + d2.*d2), -(1+nu) * d1.*d2;
        (1-nu) * (d1.*d1 + d2.*d2), -(1+nu) * d1.*d2, -(1+nu) * d1.*d2, (1-nu) * (d1.*d1 + d2.*d2)];
    Me = rho*A*[1 0; 0 1];  % 质量矩阵，忽略剪切效应
    
end

然后，你需要创建一个网格，并在循环中对每个单元调用这个函数：

% 网格设置和边界条件
nodes = ...;  % 节点坐标列表
elements = ...;  % 元素连接节点的索引
E, nu, h = ...;  % 材料属性
rho = ...;  % 密度
A = ...;  % 单元面积

K = zeros(numNodes, numNodes);
M = zeros(numNodes, numNodes);

for e = 1:length(elements)
    nodesInElement = elements(e,:);  % 当前元素的四个节点
    KeEl, MeEl = element_stiffness_mass(nodes(nodesInElement,:), E, nu, h, A);
    K(nodesInElement(1),:) = K(nodesInElement(1),:) + KeEl;
    K(:,nodesInElement(1)) = K(:,nodesInElement(1)) + KeEl';
    M(nodesInElement(1),:) = M(nodesInElement(1),:) + MeEl;
    M(:,nodesInElement(1)) = M(:,nodesInElement(1)) + MeEl';
end

% 消除自由度，添加边界条件
% ...

% 求解前10阶模态
[V, D] = eig(K, M);
D = diag(sqrt(D));  % 取模态频率的平方根
V = normalize(V, 'two_norm', [], 1);  % 归一化振型