【进阶篇】有限元法：理论与MATLAB实现

# 2.1 有限元方程的求解

有限元方程的求解是有限元法中的关键步骤，它决定了有限元模型的精度和效率。常用的求解方法有直接法和迭代法。

### 2.1.1 直接法

直接法是通过求解一个大型线性方程组来直接得到未知变量的值。常用的直接法有高斯消元法和LU分解法。直接法的优点是求解精度高，收敛速度快。但是，对于大型有限元模型，直接法需要大量的计算资源，并且计算复杂度随模型规模的增加而急剧增加。

### 2.1.2 迭代法

迭代法是通过反复迭代的方式逼近未知变量的值。常用的迭代法有雅可比迭代法、高斯-赛德尔迭代法和共轭梯度法。迭代法的优点是计算资源需求较少，并且收敛速度对于大型模型也不受影响。但是，迭代法的求解精度和收敛速度取决于迭代次数和迭代算法的选择。

# 2. MATLAB中有限元法的实现

### 2.1 有限元方程的求解

有限元法的核心在于求解描述物理问题的偏微分方程组。在MATLAB中，可以使用直接法或迭代法来求解有限元方程。

#### 2.1.1 直接法

直接法是一种一次性求解所有未知数的方法。MATLAB中常用的直接法有高斯消元法和Cholesky分解法。

**高斯消元法**

% 刚度矩阵
K = [2 -1 0; -1 2 -1; 0 -1 1];

% 载荷向量
F = [1; 0; 0];

% 求解未知位移
U = K \ F;

**逻辑分析：**

* 高斯消元法将刚度矩阵K转换为上三角矩阵，然后通过回代求解未知位移U。
* MATLAB中使用反斜杠运算符(\)来求解线性方程组，其中K\F表示求解方程组K*U = F。

**Cholesky分解法**

% 刚度矩阵
K = [2 -1 0; -1 2 -1; 0 -1 1];

% Cholesky分解
[L, U] = chol(K);

% 求解未知位移
Y = L \ F;
X = U \ Y;

**逻辑分析：**

* Cholesky分解将刚度矩阵K分解为下三角矩阵L和上三角矩阵U的乘积。
* 求解未知位移时，先将载荷向量F分解为Y，然后求解X。
* Cholesky分解比高斯消元法更稳定，但计算量更大。

#### 2.1.2 迭代法

迭代法是一种逐次逼近解的方法。MATLAB中常用的迭代法有Jacobi迭代法和共轭梯度法。

**Jacobi迭代法**

% 刚度矩阵
K = [2 -1 0; -1 2 -1; 0 -1 1];

% 载荷向量
F = [1; 0; 0];

% 初始猜测
U0 = zeros(size(F));

% 迭代次数
maxIter = 100;

% 迭代求解
for i = 1:maxIter
    for j = 1:size(K, 1)
        U0(j) = (F(j) - K(j, :) * U0) / K(j, j);
    end
end

**逻辑分析：**

* Jacobi迭代法逐行更新未知位移，直到满足收敛条件。
* MATLAB中使用for循环来实现迭代过程。

**共轭梯度法**

% 刚度矩阵
K = [2 -1 0; -1 2 -1; 0 -1 1];

% 载荷向量
F = [1; 0; 0];

% 初始猜测
U0 = zeros(size(F));

% 迭代次数
maxIter = 100;

% 共轭梯度法求解
[U, flag] = pcg(K, F, 1e-6, maxIter);

**逻辑分析：**

* 共轭梯度法是一种更有效的迭代法，它利用共轭梯度方向来加速收敛。
* MATLAB中使用pcg函数来求解共轭梯度方程组，其中K为刚度矩阵，F为载荷向量，1e-6为收敛容差，maxIter为最大迭代次数。

# 3.1 结构力学分析

有限元法在结构力学分析中得到了广泛的应用，可以解决各种复杂的结构问题，如梁的弯曲、板的振动等。

#### 3.1.1 梁的弯曲

梁的弯曲分析是结构力学中常见的问题，有限元法可以准确地计算梁的变形和内力。

% 梁的几何参数
L = 1;  % 梁长
b = 0.1; % 梁宽
h = 0.2; % 梁高
E = 200e9; % 杨氏模量
I = b*h^3/12; % 截面惯性矩

% 施加边界条件
fixed_end = 1; % 固定端
free_end = 2; % 自由端
u_fixed = 0; % 固定端的位移
theta_fixed = 0; % 固定端的转角
P = 1000; % 集中力大小

% 网格划分
n_elements = 10; % 单元数
x = linspace(0, L, n_elements+1); % 节点坐标

% 组装刚度矩阵和载荷向量
K = zeros(2*(n_elements+1)); % 刚度矩阵
F = zeros(2*(n_elements+1), 1); % 载荷向量

for i = 1:n_elements
    % 单元刚度矩阵
    k = [E*I/L^3, -E*I/L^2, -E*I/L^3, E*I/L^2;
         -E*I/L^2, 4*E*I/L^3, E*I/L^2, -2*E*I/L^3;
         -E*I/L^3, E*I/L^2, E*I/L^3, -E*I/L^2;
         E*I/L^2, -2*E*I/L^3, -E*I/L^2, 4*E*I/L^3];
    % 单元载荷向量
    f = [0; -P/2; 0; P/2];
    % 组装刚度矩阵和载荷向量
    K(2*i-1:2*i+2, 2*i-1:2*i+2) = K(2*i-1:2*i+2,